function [x,k,timev] =sbls (A,b,l,u,rep,opt)
%
%SBLS   Sparse solver for linear least squares with box constrints.
%       %Z%%M% Version %I% %G%
%       Mikael Lundquist, Linkoping University.
%       e-mail: milun@math.liu.se
%
%       x = sbls(A,b,l,u,rep,opt) solves the least squares problem
%           min || Ax-b || , subject to l <= x <= u
%            x            2
%       A        m-by-n sparse matrix, rank(A) <= n.
%       b        m-by-1 vector
%       u,l      upper and lower bound, n-by-1 vector
%       rep = 1  produces a report in each iteration (optional)
%
%       opt = 0  uses normal equations (CSNE) instead of QR factorization
%                when solving the sub problems (optional)

 nu = 3;

if nargin < 6,
  opt = 0;
end
if nargin < 5,
  rep = 0;
end

if rep==1,
  fprintf(' Iter ||x+s-u||      ||AT...||          gap\n')
end



% -- Initiation

 [m,n] = size(A);

 if length(l)==1,
   l=ones(n,1)*l;
 end
 if length(u)==1,
   u=ones(n,1)*u;
 end
 
 

% - Compute column minimum degree ordering for sparse factorization

 q = colmmd(A);
 
 A = A(:,q);


 if (opt == 0),

% - Form matrix of normal equations

	ATA = A'*A;
 end

 
 u = u(q);
 l = l(q);
 if nargin == 7,
   xexact = xexact(q);
 end

% - Shift to make the lower bounds equal to zero

 u = u - l;
 b = b - A*l;

% - Definition of the starting point (x0,w0)

 x = u/2;
 s = u-x;
 res = A'*(A*x-b);
 v = abs(res);
 y = v-res;
% ind=(res>0);
% v(ind)=res(ind);
% y(~ind)=-res(~ind);

clear res

% - Computation of the initial gap

 gap = y'*s+x'*v;
 stop1 = 0;
 stop2 = 0;
 oldgap = gap*2;
% --- Start interior point iterations

 k = 0;
 while ((abs(gap) > n*1.0e-18) | (stop1 > n*1.0e-12) | (stop2 > n*1.0e-12))

   k = k + 1;

   if k==100,break;end,

   if (oldgap/gap<1) & (gap<1e-3),
     fprintf('sbls can not converge.');
     break;
   end
   oldgap=gap;
   
% --- Predictor phase 
% - Computation of the right-hand side and deltax

   vxinv = v./x;
   ysinv = y./s;
   d2 = vxinv + ysinv;
   d= sqrt(d2);

   bd1 = b-A*x;
   bd2 = (-y + ysinv.*(u-x))./d;
   bd = [bd1 ; bd2];

% - Computation of the damped matrix and its Cholesky factor

   AD = [A ; sparse(1:length(d),1:length(d),d)];

   if (opt == 0),
      R = chol(ATA + sparse(1:length(d2),1:length(d2),d2));
   else
     R=qr(AD,0);
   end

   dx = csolve(AD,R,bd);
   
% - Computation of ds, dv and dy

   ds = (u - x - s) - dx;
   dv = - v - vxinv.*dx;
   dy = - y - ysinv.*ds;


% - Computation of the step theta


  theta = 0.9999995 *dth(x,dx,s,ds,v,dv,y,dy);

% --- Corrector phase
 
% - Computation of mu and some useful vectors

  tempgap=(x+theta*dx)'*(v+theta*dv) +(s+theta*ds)'*(y+theta*dy);
  mu = (tempgap/gap)^nu*(gap/(2*n));

  mue = mu*ones(n,1);
  c1 = (mue - dx.*dv)./x;
  c2 = (mue - ds.*dy)./s;

% - Computation of the right-hand side and dx

   bd2 = bd2 + (c1 - c2)./d;
   bd = [bd1 ; bd2 ];
   dx = csolve(AD,R,bd);

% - Computation of ds, dv and dy

   ds = (u - x - s) - dx;
   dv = c1 - v - vxinv.*dx;
   dy = c2 - y - ysinv.*ds;


% - Computation of the step theta

   theta = 0.9999995 *dth(x,dx,s,ds,v,dv,y,dy);
% - Update of x, s, v and y

   x = x + theta * dx;
   s = s + theta * ds;
   v = v + theta * dv;
   y = y + theta * dy;

% - Computation of the stopping criteria and report

   stop1 = norm (x + s - u);
   stop2 = norm (A'*(A*x - b) + y - v);
   gap = x'*v + s'*y;

   if (rep == 1),
     fprintf('  %12e   %12e   %12e\n',stop1,stop2,gap);
   end

 end %while

% - Shift the solution back

	x(q) = x + l;

% end function sbls


% --------------------------------------------------------------------------------

function teta=dth(x,deltax,y,deltay,s,deltas,v,deltav)
%dth    Help function to sbls
%       %Z%%M% Version %I% %G%
%       Mikael Lundquist, Linkoping University.
%       e-mail: milun@math.liu.se
        
i=find(deltax <-eps);if  ~isempty(i), tetax=min(-x(i)./deltax(i));end
i=find(deltay <-eps);if  ~isempty(i), tetay=min(-y(i)./deltay(i));end
i=find(deltas <-eps);if  ~isempty(i), tetas=min(-s(i)./deltas(i));end
i=find(deltav <-eps);if  ~isempty(i), tetav=min(-v(i)./deltav(i));end

teta=min([tetax tetay tetas tetav]);
if isempty(teta),teta=0,
else teta=min([teta 1]);end



% --------------------------------------------------------------------------------

function x = csolve(A,R,b,it)
%CSOLVE Computes the CSNE solution to a least squares problem.
%       %Z%%M% Version %I% %G%
%       Mikael Lundquist, Linkoping University.
%       e-mail: milun@math.liu.se
%
%       x=csolve(A,R,b) computes the CSNE solution to min||Ax-b||
%       by R'*Rx=A'b with one step of iterative refinement.

 [m,n] = size(A);
 x = zeros(n,1);

% - sne solution


 ATb=(b'*A)';
 y=R'\ATb;

 x=R\y;
% - Iterative refinement (not used)
 for i=1:-1,
   r=b-A*x;
   y=R'\(A'*r);
   dx=R\y;
   x=x+dx;
 end

% end of function solve