%  LP Solver: sphslfbf.m
%
%  See user guide at the end of this file.
%
 if exist('toler') ~= 1
   toler = 1.e-6;
 end;
 if exist('beta') ~= 1
   beta = .995;
 end;
 if exist('bindx') ~= 1
   bindx = [];
 end;
 if exist('findx') ~= 1
   findx = [];
 end;
 [nb,x]  = size(bindx);
 [nf,x]  = size(findx);
 norbc = 1+max([c;b])-min([c;b]);
 c = c/norbc;
 b = b/norbc;
 if nb > .5
     u=u/norbc;
 end
 if nf > .5,
   A = [A  -A(:,findx)];
   c = [c ; - c(findx)];
 end;
 [m,n] = size(A);
 [i,j,v] = find(A);
 ee = ones(n,1);
%
% Initialization
%
% b    = b/(norm(b)+1);
% c    = c/(norm(c)+1);
 x    = ((norm(b)+1)/sqrt(n))*ones(n,1);
 s    = ((norm(c)+1)/sqrt(n))*ones(n,1);
 if nb > 0.5,
   z    = ((norm(b)+1)/sqrt(n))*ones(nb,1);
   w    = ((norm(c)+1)/sqrt(n))*ones(nb,1);
   eez  = ones(nb,1);
 else
   z   = 1;
   w   = 0;
   eez = 0;
   rb  = 0;
   u   = 0;
 end;
 y    = zeros(m,1);
 tau0  = 1;
 kappa0= (norm(c)+1)*(norm(b)+1)/n;
% tau0  = (norm(b)+1)/sqrt(n);
% kappa0= (norm(c)+1)/sqrt(n);
 tau   = tau0;
 kappa = kappa0;
%
% Compute residuals
%
 mu0  = (x'*s+z'*w+tau*kappa)/(n+nb+1);;
 mu   = mu0;
 rp  = tau*b - A*x;
%
 if nb > 0.5,  
   rb  = x(bindx) + z - tau*u; 
 end;
%
 rd  = tau*c - A'*y -s;
 rd(bindx) = rd(bindx) + w;
 obp = c'*x;
 obd = b'*y-u'*w;
 rg  = obp - obd + kappa;
 zh   = (obp-obd)/tau;
%
% Iteration
%
 gamma=1/sqrt(n+nb+1);
 go   = 1;
 iter=1;
 while go >= toler,
%  [x,y,s,tau,kappa]=sphslfsubb(A,b,c,x,y,s,tau,kappa,mu,rp,rd,rg,beta,gamma);
   sphslfsubf
%
%  Compute new residuals
%
%   tau/kappa
   mu  = (x'*s+z'*w+tau*kappa)/(n+nb+1),
   rp  = tau*b - A*x;
   if nb > 0.5,  
     rb  = x(bindx) + z - tau*u; 
   end;
   rd  = tau*c - A'*y -s;
   rd(bindx) = rd(bindx) + w;
   obp = c'*x;
   obd = b'*y-u'*w;
   rg  = obp - obd + kappa;
%   pause
   go = max([norm([rp;rb])/(tau+norm([x;z])); norm(rd)/(tau+norm([s;w]))]);
   go = max([go; abs(obp-obd)/(tau+abs(obd))]);
%
%  Check infeasibility
% 
   if (tau*kappa0/(tau0*kappa) < toler) & (mu/mu0 < toler/n),
     if (obp < 0) & (obd > 0),
      disp('Both primal and dual are infeasible.'); go=-1;
     end;
     if obp < 0, 
      disp('Dual is infeasible.')  ; go=-1; 
     end;
     if obd > 0,
      disp('Primal is infeasible.'); go=-1;
     end;
   end;
   zh   = [zh (obp-obd)/tau];
%   
%  Adjust centering weight
%
   gamma = 1/sqrt(n+nb+1);
   if nb > 0.5, 
     if min([x.*s;z.*w;tau*kappa])/mu >= (1-beta),
       gamma = 0;
     end;
   else
     if min([x.*s;tau*kappa])/mu >= (1-beta),
       gamma = 0;
     end;
   end;
%
   iter = iter + 1;
 end;
%
% Make an output
%
% b    = b/(1-norm(b));
% c    = c/(1-norm(c));
 s(bindx) = s(bindx) - w;
 if go >= 0,
   x = norbc*x/tau;
   z = norbc*z/tau;
   s = norbc*s/tau;
   y = norbc*y/tau;
   w = norbc*w/tau;
 end
 if nf > .5,
   n = n - nf;
   x(findx) = x(findx) - x(n+1:n+nf);
   x = x(1:n);
   s = s(1:n);
   A = A(:,1:n);
   c = c(1:n);
 end;
 c = c*norbc;
 b = b*norbc;
 if nb > .5
     u=u*norbc;
 end;
 clear norbc ee eez rp rb rd rg i j v
%return
%
%  This program solves linear program:
%
%      minimize    c'*x
%      subject to  A*x=b, 
%                  x_f free, x_p >= 0, u >= x_b >= 0.
%
%  Input 
%      A: sparse constraint matrix
%      b: constraint right-hand column vector
%      c: objective column vector
%
%  Optional Input
%      u: upp bound vector for a set of nonnegative variable x_b, whose
%         indices are represented by bindx.
%      bindx: the index set for those upp-bounded nonnegative variables
%             Default value: []
%      findx: the index set for free variables x_f
%             Default value: []
%      toler: stopping tolerance, the objective value close to the      
%             optimal one in the range of the tolerance. 
%             Default value: 1.0e-6
%      beta : step size, 0 < beta < 1. 
%             Default value: .995
%
%
%  Output
%     x  : optimal solution
%     y  : optimal dual solution (shadow price) for equality constraints
%     w  : optimal dual solution for uppbound constraints
%     s  : optimal dual slack solution
%     zh : (infeasible) primal-dual gap history vs iteration. 
%     
%  Subroutine called : sphslfsubf
%
%
%  How to use it? Just type "sphslfbf" and hit the "return".
%
%
%  This program is the MATLAB implementation of the homogeneous and
%  selfdual interior-point algorithm for LP. The result is a
%  strict complementarity solution, i.e., a solution in the relative
%  interior of the LP optimal face.
%
%  Technical References
%
%  Y. Ye, M. J. Todd and S. Mizuno, "An $O(\sqrt{n}L)$-iteration
%  homogeneous and self-dual linear programming algorithm," (1992),
%  to appear in Math. of OR.
%
%  
%  X. Xu, P. Hung and Y. Ye, "A simplified homogeneous and self-dual
%  linear programming algorithm and its implementation," manuscript,
%  Department of Management Sciences, The University of Iowa
%  (Iowa City, IA 52242, 1993).
%
%  Changes from basic version: 
%     
%   1) Adaptively adjust gamma by checking min(x.*s)/mu. If it is
%      bounded from below by constant (1-beta)/10, then gamma=0.
%
%   2) Set eta = 1 to speed up feasibility improvement.
%
%   3) Sparse version.
%
%   4) Handle up bounds x_i <= u_i, for i in bindx.
%
%   5) Handle free variables x_i for i in findx.
%
%   6) cleaned on 11/16/93