lipsol.m

线性规划算法

      %    removed since the indices still reflect their existence 
      iAeqindex = Aindex(find (Aindex >= Aeqstart & Aindex <= Aeqend) ) ... 
         - extraAineqvars - numRowsAineq; 
       
      nnzAineqindex = nnz(iAineqindex); 
      nnzAeqindex = nnz(iAeqindex); 
       
      lambda.ineqlin(iAineqindex) = full(-y(1:nnzAineqindex)); 
      yAeq = y(nnzAineqindex+extraAineqvars+1:end); 
      lambda.eqlin(iAeqindex) = full(-yAeq(1:nnzAeqindex)); 
       
      nnzubindex = nnz(ubindexfinal); 
      lambda.upper(ubindexfinal) = full(w(windexfinal(1:nnzubindexfinal))); 
      % if we recast some upper bound constraints as xpos/xneg constraints 
      extraindex = Aindex((Aineqend+1):(Aeqstart-1));  
      lambda.upper(ilbinfubfin) = full(-y(extraindex)); 
      lambda.lower(lbindexfinal) = full(z(zindexfinal(1:nnzlbindexfinal))); 
       
      if Fixed_exist & ~Sgtons_exist 
         Alambda = Aineq_orig'*lambda.ineqlin + Aeq_orig'*lambda.eqlin; 
         lambda.lower(fixedlower) = f_orig(fixedlower) + Alambda(fixedlower); 
         lambda.upper(fixedupper) = - f_orig(fixedupper) - Alambda(fixedupper); 
          
      elseif Sgtons_exist & ~Fixed_exist 
         % sgrows: what eqlin lambdas we are solving for (row in A) 
         % sgcols: what column in A that singleton is in 
         result = -f_orig - Aineq_orig'*lambda.ineqlin - Aeq_orig'*lambda.eqlin ... 
            + lambda.lower - lambda.upper; 
         lambda.eqlin(sgrows) = result(sgcols); 
      elseif Fixed_exist & Sgtons_exist %  
         % solve for singletons first 
         Alambda = Aineq_orig'*lambda.ineqlin + Aeq_orig'*lambda.eqlin; 
         result = -f_orig - Aineq_orig'*lambda.ineqlin - Aeq_orig'*lambda.eqlin ... 
            + lambda.lower - lambda.upper; 
         lambda.eqlin(sgrows) = result(sgcols); 
         % now the fixed variables multipliers 
         Alambda = Aineq_orig'*lambda.ineqlin + Aeq_orig'*lambda.eqlin; 
         lambda.lower(fixedlower) = f_orig(fixedlower) + Alambda(fixedlower); 
         lambda.upper(fixedupper) = - f_orig(fixedupper) - Alambda(fixedupper); 
      end 
   end 
end  % ComputeLambda 
 
if (converged) 
   exitflag = 1; 
elseif (iter == maxiter+1) 
   message = sprintf('Maximum number of iterations exceeded.\n'); 
   exitflag = 0; 
else 
   exitflag = -1; 
end 
 
if (diagnostic_level >= 1) 
   if (diagnostic_level > 1) 
      fprintf('\n'); 
   end 
   if (converged) 
      fprintf('%s\n',message); 
   else 
      fprintf('Exiting: %s\n',message); 
   end 
end 
 
 
output.iterations = iter; 
output.cgiterations = cgiter; 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LIPSOL Subfunctions %%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function message = detectinf(tol,Hist) 
%DETECTINF  Used to detect an infeasible problem. 
%   One or more of these 5 convergence criteria either blew up or stalled: 
%   trerror  total relative error      max([rrb rrf rru rdgap]) 
%   rrb      relative primal residual  norm(A*x-b) / max(1,norm(b)) 
%   rrf      relative dual residual    norm(A'*y+z-w-f) / max(1,norm(f)) 
%   rru      relative upper bounds     norm(spones(s).*x+s-ub) / max(1,norm(ub)) 
%   rdgap    relative objective gap    abs(fval-dualval) / max([1 abs(fval) abs(dualval)]) 
%   where 
%   fval     primal objective          f'*x 
%   dualval  dual objective            b'*y - w'*ub       
% 
%   Note: these convergence criteria are stored in the first 5 rows of Hist, 
%   while fval and dualval are the 6th and 7th rows of Hist respectively. 
 
% minimum of all relative primal residuals (take upper bounds into account) 
minrp = min(Hist(2,:)+Hist(4,:)); 
% minimum of all relative dual residuals 
minrd = min(Hist(3,:));  
 
if minrp < tol 
   message = sprintf(['         the dual appears to be infeasible (and the primal unbounded).' ,... 
              '      \n         (The primal residual < TolFun=%3.2e.)'],tol); 
   return 
end 
if minrd < tol 
   message = sprintf(['         the primal appears to be infeasible (and the dual unbounded).', ... 
              '\n         (The dual residual < TolFun=%3.2e.)'],tol); 
   return 
end 
 
tol10 = 10 * tol; 
sqrttol = sqrt(tol); 
if ((minrp < tol10) & (minrd > sqrttol)) 
   message = sprintf(['         the dual appears to be infeasible (and the primal unbounded) since', ... 
         '\n         the dual residual > sqrt(TolFun)=%3.2e.' ,... 
         '\n         (The primal residual < 10*TolFun=%3.2e.)'], ... 
      sqrttol,tol10); 
   return 
end 
if ((minrd < tol10) & (minrp > sqrttol)) 
   message = sprintf(['         the primal appears to be infeasible (and the dual unbounded) since', ... 
             '\n         the primal residual > sqrt(TolFun)=%3.2e.', ... 
             '\n         (The dual residual < 10*TolFun=%3.2e.)'], ... 
          sqrttol,tol10); 
   return 
end 
 
iter = size(Hist,2); 
fval = Hist(6,iter); 
dualval = Hist(7,iter); 
if ((fval < -1.e+10) & (dualval < 1.e+6)) 
    message = sprintf(['         the dual appears to be infeasible and the primal unbounded since' ,... 
          '\n         the primal objective < -1e+10',... 
          '\n         and the dual objective < 1e+6.']); 
   return 
end 
if ((dualval > 1.e+10) & (fval > -1.e+6)) 
   message = sprintf(['         the primal appears to be infeasible and the dual unbounded since' ,... 
                   '\n         the dual objective > 1e+10',... 
          '\n         and the primal objective > -1e+6.']); 
   return 
end 
 
message = sprintf(['         both the primal and the dual appear to be infeasible.\n']); 
    
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function [Rxz,Rsw,dgap] = complementy(x,z,s,w,Ubounds_exist) 
%COMPLEMENTY  Evaluate the complementarity vectors and gap. 
 
Rxz = x .* z; 
if (Ubounds_exist) 
   Rsw = s .* w; 
else 
   Rsw = []; 
end 
dgap = full(sum([Rxz; Rsw])); 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function [Rb,Rf,rb,rf,ru,Ru] = feasibility(A,b,f,ub,Ubounds_exist,x,y,z,s,w) 
%FEASIBILITY  Evaluate feasibility residual vectors and their norms. 
 
Rb = A*x - b; 
Rf = A'*y + z - f; 
if (Ubounds_exist) 
   Rf = Rf - w; 
   Ru = spones(s).*x + s - ub; 
else 
   Ru = []; 
end 
rb = norm(Rb); 
rf = norm(Rf); 
ru = norm(Ru); 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function [trerror,rrb,rrf,rru,rdgap,fval,dualval] = ... 
   errornobj(b,rb,bnrm,f,rf,fnrm,ub,ru,unrm,dgap,Ubounds_exist,x,y,w) 
%ERRORNOBJ  Calculate the total relative error and objective values. 
 
rrb = rb/max(1,bnrm); 
rrf = rf/max(1,fnrm); 
rru = 0; 
fval = full(f'*x); 
dualval = full(b'*y); 
if (Ubounds_exist)  
   rru = ru/(1+unrm);  
   dualval = dualval - w'*ub; 
end 
if (min(abs(fval),abs(dualval)) < 1.e-4) 
   rdgap = abs(fval-dualval); 
else 
   rdgap = abs(fval-dualval)/max([1 abs(fval) abs(dualval)]); 
end 
trerror = max([rrb rrf rru rdgap]); 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function Y = reciprocal(X) 
%RECIPROCAL  Invert the nonzero entries of a matrix elementwise. 
%   Y = RECIPROCAL(X) has the same sparsity pattern as X 
%	 (except possibly for underflow). 
 
if (issparse(X)) 
   [m, n]  = size(X); 
   [i,j,Y] = find(X); 
   Y = sparse(i,j,1./Y,m,n); 
else 
   Y = 1./X; 
end 
 
Y = min(Y,1e8);  
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function [P,U] = getpu(A,vmid,Dense_cols_exist,idense,ispars) 
%GETPU  Compute Matrix P and U 
 
[m,n] = size(A); 
 
if (Dense_cols_exist) 
   ns = length(ispars); 
   nd = length(idense);  
   Ds = spdiags(vmid(ispars),0,ns,ns); 
   P = A(:,ispars) * Ds * A(:,ispars)'; 
   Dd = spdiags(sqrt(vmid(idense)),0,nd,nd); 
   U = full(A(:,idense) * Dd); 
else 
   P = A * spdiags(vmid,0,n,n) * A'; 
   U = []; 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function [dx,dy,dz,ds,dw,Sherman_OK,Mdense,Rinf,cgiter] = ... 
   direction(A,P,U,Rb,Rf,Ru,Rxz,Rsw,vmid,xn1,sn1,z,w,mu,flag,... 
   Sherman_OK,Dense_cols_exist,Ubounds_exist,Mdense,perm,Rinf,Hist,cgiter) 
%DIRECTION  Compute search directions 
 
if (mu ~= 0) 
   Rxz = Rxz - mu; 
end 
Rf = Rf - Rxz .* xn1; 
if (Ubounds_exist) 
   if (mu ~= 0) 
      Rsw = Rsw - mu; 
   end 
   Rf = Rf + (Rsw - Ru .* w) .* sn1; 
end 
rhs = -(Rb + A * (vmid .* Rf)); 
 
if (Dense_cols_exist) 
   [dy,Sherman_OK,Mdense,Rinf,cgiter] = ... 
      densol(P,U,full(rhs),flag,Sherman_OK,Mdense,perm,Rinf,Hist,cgiter); 
else 
   if (flag == 1) 
      Rinf = cholinc(sparse(P(perm,perm)),'inf'); 
   end 
   warnstr = warning; 
   warning off 
   dy(perm,1) = Rinf \ (Rinf' \ full(rhs(perm))); 
   warning(warnstr) 
end 
 
dx = vmid .* (A' * dy + Rf); 
dz = -(z .* dx + Rxz) .* xn1; 
if (Ubounds_exist) 
   ds = -(dx .* spones(w) + Ru); 
   dw = -(w .* ds + Rsw) .* sn1; 
else 
   ds = []; 
   dw = []; 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function [ap,ad] = ratiotest(dx,dz,ds,dw,Ubounds_exist,x,z,s,w) 
%RATIOTEST  Ratio test 
 
ap = -1/min([dx(find(x))./nonzeros(x); -0.1]); 
ad = -1/min([dz(find(z))./nonzeros(z); -0.1]); 
if (Ubounds_exist) 
   as = -1/min([ds(find(s))./nonzeros(s); -0.1]); 
   aw = -1/min([dw(find(w))./nonzeros(w); -0.1]); 
   ap = min(ap, as); 
   ad = min(ad, aw); 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function [x,Sherman_OK,Mdense,Rinf,cgiter] = ... 
   densol(P,U,b,flag,Sherman_OK,Mdense,perm,Rinf,Hist,cgiter) 
%DENSOL  Solve linear system with dense columns using either the 
%    Sherman-Morrison formula or preconditioned conjugate gradients 
 
bnrm = norm(b); 
x = zeros(size(b)); 
tol1  = min(1.e-2, 1.e-7*(1+bnrm)); 
tol2  = 10*tol1; 
tolcg = tol1; 
 
iter = size(Hist,2); 
if (iter > 0) 
   trerror = Hist(1,iter); 
   tolcg = min(trerror, tolcg); 
end 
 
if (Sherman_OK) 
   [x,Mdense,Rinf] = sherman(P,U,b,flag,Mdense,perm,Rinf); 
end 
resid = norm(b - (P*x + U*(U'*x))); 
if (resid > tol1) 
   Sherman_OK = 0; 
end 
if (resid < tol2) 
   return 
end 
if (~Sherman_OK) 
   if (flag == 1) 
      Rinf = cholinc(sparse(P(perm,perm)),'inf'); 
   end 
   warnstr = warning; 
   warning off 
   [x(perm,1),pcg_flag,pcg_relres,pcg_iter] = pcg(@lipsol_fun,b(perm), ... 
      tolcg/bnrm,250,Rinf',Rinf,x(perm),P(perm,perm),U(perm,:)); 
   cgiter = cgiter + pcg_iter; 
   warning(warnstr) 
end 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
function F = lipsol_fun(x,P,U) 
% Local function to solve using PCG (conjugate gradients) 
%  lipsol_fun = inline('P*x+U*(U''*x)','x','P','U'); 
F = P*x+U*(U'*x); 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
function [x,Mdense,Rinf] = sherman(P,U,b,flag,Mdense,perm,Rinf) 
%SHERMAN  Use the Sherman-Morrison formula to "solve" (P + U*U')x = b 
%   where P is sparse and U is low-rank. 
 
if (flag == 1) 
   nd = size(U,2);  
   PinvU = zeros(size(U)); 
   Rinf = cholinc(sparse(P(perm,perm)),'inf'); 
   warnstr = warning; 
   warning off 
   PinvU(perm,:) = Rinf \ (Rinf' \ U(perm,:)); 
   warning(warnstr) 
   Mdense = eye(nd) + U' * PinvU; 
end 
 
warnstr = warning; 
warning off 
x(perm,1) = Rinf \ (Rinf' \ b(perm)); 
tmp = U * (Mdense \ (U' * x)); 
tmp(perm,1) = Rinf \ (Rinf' \ tmp(perm)); 
warning(warnstr) 
x = x - tmp; 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%% End of LIPSOL and its subfunctions %%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%