lipsol.m

数学建模的源代码

         lambda.eqlin(sgrows) = result(sgcols);
      end
   else
      iAineqindex = Aindex( find (Aindex >= Aineqstart & Aindex <= Aineqend) );
      % Subtract off extraAineqvars and numRowsAineq even if those were
      %    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 (diagnostic_level >= 1)
   if (diagnostic_level > 1)
      fprintf('\n');
   end
   if (converged)
      fprintf('%s\n',message);
   else
      fprintf('Exiting: %s\n',message);
   end
end

if (converged)
   exitflag = 1;
elseif (iter == maxiter+1)
   exitflag = 0;
else
   exitflag = -1;
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,1e16);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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
   lipsol_fun = inline('P*x+U*(U''*x)','x','P','U');
   [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 [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 %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%