lipsol.m

线性规划算法

      return 
   end 
   xzrcol = zeros(size(izrcol))... 
      + (f(izrcol) < 0).*ub(izrcol)... 
      + (f(izrcol) > 0).*lb(izrcol); 
   inzcol = find(~zrcol); 
   if computeLambda 
      % assign a multiplier according to -f(i), leave the other zero  
      lowermultnonzero = ( f(1:n_orig,1) >= 0 ) & zrcol(1:n_orig,1);  
      uppermultnonzero = ( f(1:n_orig,1) < 0 ) &  zrcol(1:n_orig,1);   
      lambda.lower(lbindex(lowermultnonzero)) = f(lowermultnonzero);  
      lambda.upper(ubindex(uppermultnonzero)) = -f(uppermultnonzero); 
      lbindex = lbindex(inzcol);                         
      ubindex = ubindex(inzcol); 
      zindex = zindex(inzcol); 
      windex = windex(inzcol); 
   end  
   A = A(:,inzcol); 
   f = f(inzcol); 
   lb = lb(inzcol);                         
   ub = ub(inzcol); 
   if (diagnostic_level >= 4) 
      fprintf(['  Deleting %i all-zero columns from' ... 
            ' constraint matrix.\n'],nnz(zrcol)); 
   end 
   data_changed = 1; 
   n = size(f,1); 
end 
 
% solve singleton rows 
Sgtons_exist = 0; 
singleton = (rnnzct == 1); 
nsgrows = nnz(singleton); 
if (nsgrows >= max(1, .01*size(A,1))) 
   Sgtons_exist = 1; 
   isgrows = find(singleton); 
   iothers = find(~singleton); 
   if (diagnostic_level >= 4) 
      fprintf('  Solving %i singleton variables immediately.\n', ... 
         nsgrows); 
   end 
   Atmp = A(isgrows,:); 
   Atmp1 = spones(Atmp); 
   btmp = b(isgrows); 
   if (nsgrows == 1) 
      isolved  = find(Atmp1);  
      insolved = find(Atmp1 == 0);  
      xsolved  = b(isgrows)/Atmp(isolved); 
   else 
      colnnzct = sum(Atmp1); 
      isolved  = find(colnnzct); 
      insolved = find(colnnzct == 0); 
      [ii, jj] = find(Atmp);  
      Atmp = Atmp(ii,jj); 
      btmp = btmp(ii); 
      xsolved  = btmp./diag(Atmp); 
      if (any(colnnzct > 1)) 
         repeat = diff([0; jj]) == 0; 
         for i = 1:length(xsolved) - 1 
            if repeat(i+1) & xsolved(i+1) ~= xsolved(i) 
               if (diagnostic_level >= 1) 
                  fprintf(['Exiting due to infeasibility:' ... 
                        '   Singleton variables in equality constraints are not feasible.\n']); 
               end 
               xsol = []; 
               fval = []; 
               lambda.ineqlin = []; 
               lambda.eqlin = []; 
               lambda.upper = []; 
               lambda.lower = []; 
               exitflag = -1; 
               output = []; 
               return 
            end 
         end 
         ii = find(~repeat); 
         jj = ii; 
         Atmp = Atmp(ii,jj); 
         btmp = btmp(ii); 
         xsolved  = btmp./diag(Atmp); 
      end 
   end 
    
   % check that singleton variables are within bounds 
   if ((any(xsolved < lb(isolved))) | ... 
         (any(xsolved > ub(isolved)))) 
      if (diagnostic_level >= 1) 
         fprintf(['Exiting due to infeasibility:' ... 
               '   %i singleton variables in the equality constraints are not' ... 
               ' within bounds.\n'], ... 
            sum((xsolvedub(isolved)))); 
      end 
      xsol = []; 
      fval = []; 
      lambda.ineqlin = []; 
      lambda.eqlin = []; 
      lambda.upper = []; 
      lambda.lower = []; 
      exitflag = -1; 
      output = []; 
      return 
   end 
    
   if computeLambda 
      % Compute which multipliers will need to be computed 
      %   sgrows: what eqlin lambdas we are solving for (row in A) 
      %   sgcols: what column in A that singleton is in 
      sgAindex = Aindex(singleton); 
      sgrows = sgAindex( (sgAindex >= Aeqstart & sgAindex <= Aeqend) ) ... 
         - extraAineqvars - numRowsAineq; 
 
      % Only want to extract out indices for the original variables, not slacks. 
      %  Must also subtract out the removed variables from zero columns and fixed variables. 
      [iii,jjj]=find(Atmp1'); 
      sgcols = lbindex( iii( iii <= (n_orig - nnz(zrcol) - nnz(fixed)) ) ); 
      if length(sgrows) ~= length(sgcols) 
         errmsg = ... 
            sprintf(['Trying to compute lagrange multipliers. \n',... 
            'When removing singletons in equality constraints, sgrows does not match sgcols.\n ',... 
            'Please report this error to The MathWorks.']); 
         error(errmsg) 
      end 
      Aindex = Aindex(iothers); 
      yindex = yindex(iothers); 
      lbindex = lbindex(insolved); 
      ubindex = ubindex(insolved); 
      zindex = zindex(insolved); 
      windex = windex(insolved); 
   end 
   % reform the unsolved part of the problem 
   b = b(iothers,1) - A(iothers,isolved)*xsolved; 
   A = A(iothers, insolved); 
   f = f(insolved); 
   lb = lb(insolved); 
   ub = ub(insolved); 
   data_changed = 1; 
end 
 
n = length(ub); 
nbdslack = n; 
boundsInA = 0; 
if (isempty(A)) % constraint matrix has been cleared out 
   boundsInA = 1; 
   ubfinite = find(isfinite(ub)); 
   lbfinite = find(isfinite(lb)); 
   lubfinite = length(ubfinite); 
   llbfinite = length(lbfinite); 
   Aineq = [sparse(1:lubfinite,ubfinite,1,lubfinite,n); ... 
         sparse(1:llbfinite,lbfinite,-1,llbfinite,n)]; 
   % Aineqend and Aeqend are less than Aeq(Aineq)start since these are lb,ub constraints 
   %   No singletons, due to slacks; no zero rows; full rank due to slacks; 
   %   could have fixed variables 
   bineq = [ub(ubfinite); -lb(lbfinite)]; 
   if (diagnostic_level >= 4) 
      fprintf(['  No constraints: converting %d finite upper bounds and' ... 
            ' %d finite lower bounds into %d inequality constraints.\n'], ... 
         lubfinite,llbfinite,lubfinite+llbfinite); 
   end 
    
   mineq = size(Aineq,1); 
   m = mineq; 
   f = [f; sparse(mineq,1)]; 
   A = [Aineq speye(mineq)]; 
   b = bineq; 
   lb = [lb; sparse(mineq,1)]; 
   ub = [ub; repmat(Inf,mineq,1)]; 
   if computeLambda 
      Aindex = 1:(lubfinite+llbfinite); 
      Aineqstart = 1; Aineqend = 0;  % empty 
      Aeqstart = 1; Aeqend = 0;  % empty 
   end 
   if (diagnostic_level >= 4) 
      fprintf(['  Adding %d slack variables, one for each inequality' ... 
            ' constraint (bounds), to the existing\n  set of %d variables,'],mineq,n); 
   end 
   n = n + mineq; 
   if (diagnostic_level >= 4) 
      fprintf([' resulting in %d equality constraints on %d variables.\n'], ... 
         mineq,n); 
   end    
end 
 
% shift nonzero lower bounds 
Lbounds_non0 = any(lb ~= 0); 
if (Lbounds_non0) 
   b = b - A*lb; 
   data_changed = 1; 
   if (diagnostic_level >= 4) 
      fprintf(['  Shifting %i non-zero lower bounds to zero.\n'], ... 
         full(sum(lb ~= 0))); 
   end 
end 
% find upper bounds 
nub = 0; 
iubounds = (ub ~= Inf);    
if computeLambda 
   ubindex = ubindex(iubounds);   
   windex(~iubounds) = 0;        % not removing from ub/w, just zeroing out 
end 
Ubounds_exist = full(any(iubounds)); 
if (Ubounds_exist) 
   ub(~iubounds) = 0; 
   ub = sparse(iubounds .* (ub-lb)); 
   nub = nnz(ub); 
end 
m = size(A,1); 
 
if computeLambda 
   % now ignore the lb inf bounds 
   [stmp,itmp]=setdiff(lbindex,find(lbnotinf)); 
   zindex(itmp) = 0;  % zero out inf bounds 
   lbindex = intersect(find(lbnotinf),lbindex); 
end 
 
% Scale the problem 
badscl = 1.e-4; 
col_scaled = 0; 
absnzs = abs(nonzeros(A)); 
thescl = min(absnzs)/max(absnzs); 
 
if (thescl < badscl) 
   if (diagnostic_level >= 4) 
      fprintf(['  Scaling problem by square roots of infinity norms of rows and' ... 
            ' columns of constraint matrix.\n']); 
   end 
    
   % ----- scaling vectors ------ 
   absA = abs(A); 
   colscl = full(sqrt(max(absA, [], 1)')); 
   rowscl = full(sqrt(max(absA',[], 1)')); 
    
   % ----- column scaling ----- 
   if (Ubounds_exist) 
      ub = ub .* colscl; 
   end 
    
   colscl = reciprocal(colscl);  
   A = A * spdiags(colscl,0,n,n); 
    
   f = f .* colscl; 
   col_scaled = 1; 
    
   % ----- row scaling ----- 
   rowscl = reciprocal(rowscl); 
   A = spdiags(rowscl,0,m,m) * A; 
   b = b .* rowscl; 
   bnrm = norm(b); 
   if (bnrm > 0)  
      q = median([1 norm(f)/bnrm 1.e+8]); 
      if (q > 10) 
         A = q * A; 
         b = q * b; 
      end 
   end 
   data_changed = 1; 
end 
 
% Not all variables are fixed, singleton, etc, i.e. more variables left to solve for 
if ~isempty(A) 
    
   idense = []; 
   ispars = 1:n; 
   Dense_cols_exist = 0; 
   nzratio = 1; 
    
   if (m > 2000) 
      nzratio = 0.05; 
   elseif (m > 1000) 
      nzratio = 0.10; 
   elseif (m >  500) 
      nzratio = 0.20; 
   end 
    
   if (nzratio < 1) 
      checking = (sum(spones(A))/m <= nzratio); 
      if (any(checking == 0)) 
         Dense_cols_exist = 1; 
         % checking will be mostly ones, (1-checking) will be very sparse 
         idense = find(sparse(1-checking));   % Dense  column indices 
         ispars = find(checking);             % Sparse column indices     
      end 
   end 
    
   if ((diagnostic_level >= 4) & ~isempty(idense)) 
      fprintf(['  Separating out %i dense' ... 
            ' columns from constraint matrix.\n'],length(idense)); 
   end 
    
   tau0 = .9995;		% step percentage to boundary 
   phi0 = 1.e-5;		% initial centrality factor 
    
   m = size(A,1); 
   nt = n + nub; 
   bnrm = norm(b); 
   fnrm = norm(f); 
   Hist = []; 
    
   if (Ubounds_exist) 
      unrm = norm(ub); 
   else 
      unrm = []; 
   end 
    
   if (Dense_cols_exist) 
      perm = colmmd(A(:,ispars)'); 
   else 
      perm = colmmd(A'); 
   end 
    
   y = zeros(m,1);  
   pmin = max(bnrm/100, 100); 
   dmin = fnrm*.425; 
   dmin = max(dmin, pmin/40); 
   pmin = min(pmin, 1.e+4); 
   dmin = min(dmin, 1.e+3); 
    
   Mdense = []; 
   if (Dense_cols_exist) 
      Sherman_OK = 1; 
   else 
      Sherman_OK = 0; 
   end 
   cgiter = 0;    
   [P,U] = getpu(A,ones(n,1),Dense_cols_exist,idense,ispars); 
   if (Dense_cols_exist) 
      Rinf = []; 
      [x,Sherman_OK,Mdense,Rinf,cgiter] = ... 
         densol(P,U,full(b),1,Sherman_OK,Mdense,perm,Rinf,Hist,cgiter); 
      x = A' * x; 
   else 
      Rinf = cholinc(sparse(P(perm,perm)),'inf'); 
      warnstr = warning; 
      warning off 
      temp(perm,1) = Rinf \ (Rinf' \ full(b(perm))); 
      warning(warnstr) 
      x = A' * temp; 
   end 
   pmin = max(pmin, -min(x)); 
   x = max(pmin,x); 
    
   z = full((f+dmin).*(f > 0) + dmin*(-dmin < f & f <= 0) - f.*(f <= -dmin)); 
    
   if (Ubounds_exist) 
      s = spones(ub).*max(pmin, ub-x); 
      w = spones(ub).*( dmin*(f > 0) + ... 
         (dmin - f).*(-dmin < f & f <= 0) - 2*f.*(f <= -dmin) ); 
   else 
      s = []; 
      w = []; 
   end 
    
   [Rxz,Rsw,dgap] = complementy(x,z,s,w,Ubounds_exist); 
    
   [Rb,Rf,rb,rf,ru,Ru] = feasibility(A,b,f,ub,Ubounds_exist,x,y,z,s,w); 
   if any(isnan([rb rf ru])) 
       if (diagnostic_level >= 1) 
           fprintf(['Exiting: cannot converge because the primal residual, dual residual,  \n' ... 
                    '         or upper-bound feasibility is NaN.\n']); 
       end 
       xsol = []; 
       fval = []; 
       lambda.ineqlin = []; 
       lambda.eqlin = []; 
       lambda.upper = []; 
       lambda.lower = []; 
       exitflag = -1; 
       output = []; 
       return 
   end 
    
   [trerror,rrb,rrf,rru,rdgap,fval,dualval] = ... 
      errornobj(b,rb,bnrm,f,rf,fnrm,ub,ru,unrm,dgap,Ubounds_exist,x,y,w);