lipsol.m

数学建模的源代码

      lambda.lower = [];
      exitflag = -1;
      output = [];
      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 >= 0 ) & zrcol);
      uppermultnonzero = (( f < 0 ) &  zrcol);
      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;
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((xsolved<lb(isolved))|(xsolved>ub(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,n] = size(A);

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)'));
   rowscl = full(sqrt(max(absA')'));
   
   % ----- 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,n] = size(A);
   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]))
      xsol = [];
      fval = [];
      lambda.ineqlin = [];
      lambda.eqlin = [];
      lambda.upper = [];
      lambda.lower = [];
      exitflag = -1;
      output = [];
      return
   end