Logo Search packages:      
Sourcecode: semidef-oct version File versions  Download package

lin_prog.m

function [x,z,info] = lin_prog(A,b,c);

% [x,z,info] = lin_prog(A,b,c);
% 
% solves the linear program
%
% (primal) minimize    c'*x
%          subject to  A*x + b >= 0
%
% (dual)   maximize    -b'*z
%          subject to  A'*z = c
%                      z >= 0
%
% Input arguments:
% A:     nxm-matrix.
% b:     n-vector.
% c:     m-vector.
%
% Output arguments
% x:     primal solution.
% z:     dual solution.
% info:  string.  info = 'optimal', or 'may be infeasible', or 
%        'may be unbounded below', or 'maxiters exceeded'.
%
%
% The code adds an upper bound to the primal constraints,
% and it really solves  the problem
%
% (primal) minimize    c'*x
%          subject to  A*x + b >= 0
%                      e'*(A*x + b) <= M
%
% (dual)   maximize    -b'*(z-w*e) - M*w
%          subject to  A'*(z-w*e) = c
%                      z >= 0, w >= 0
%
% with M = 1e3 * || [A b] ||_1
% and e is the vector with all components one.
% 
% If no optimal solution can be found with strict inequality
% e'*(b-A*x) < M, then the original problem is infeasible, or
% unbounded below, or the default value of M is too small.  
%

[n,m] = size(A);

if (size(b,1) < size(b,2)),  b = b';  end;
if ( size(b,1) ~= n ) | (size(b,2) ~= 1),
   error('b must be an n-vector.');
end;
if (size(c,1) < size(c,2)),  c = c';  end;
if ( size(c,1) ~= m ) | (size(c,2) ~= 1),
   error('c must be an m-vector.');
end;

M = 1e3 * max(sum(abs([A b])));

%
% phase 1
% 

disp(' '); disp(' PHASE 1.');
abstol = 1e-8;  nu = 10.0;  maxiters = 100;
[x0,Z0,z0,ul,iters] = ...
   phase1([b A], ones(n,1), M, nu, abstol, maxiters);

if iters == maxiters,
   info = 'maxiters exceeded';
   x = ; z = ;
   return;
end;

if (ul(1) > 0.0) 
   info = 'may be infeasible';
   x=; z=;
   return;
end;


%
% phase 2
% 

disp(' '); disp(' PHASE 2.');
M = max(M, 1e3*sum(b+A*x0));    % M must be greater than e'(b-A*x0) 
                                % for bigM.m 
abstol = 1e-8;  reltol = 1e-3;  nu = 10.0;  maxiters = 100;
[x,z,w,ul,iters] = ...
   bigM([b A], ones(n,1), c, x0, M, nu, abstol, reltol, 0.0, maxiters);

if iters == maxiters
   info = 'maxiters exceeded';
   x=[]; z=[];
   return;
end;

if sum(b+A*x) > 0.9*M 
   info = 'may be unbounded below';
   x=[]; z=[];
   return;
end;

info = 'optimal';

Generated by  Doxygen 1.6.0   Back to index