function[x,residnorm,residual,fnrm,it,npcg,ex,LAMBDA]=sllsbox(A,b,lb,ub,xstart,verb,options)
%SLLSBOX Linear least-squares with bounds
%
% x=slls(A,b,lb,ub) returns the solution to the
% box-constrained linear least-squares problem,
%
%        min { ||Ax - b||^2 : lb <= x <= ub}.
%
% where A is a matrix with more rows than cols. (may be virtual),
%
% [x,val] = sllsbox(b,A,xstart,options,data) returns the value of
% the objective function, ||Ax - b||^2.
%
% [x,val,fnrm] = sllsbox(b,A,xstart,options,data) returns the 
% optimality vector
%
% [x,val,fnrm,it]=sllsbox(b,A,xstart,options,data) returns the
% number of iterations used.
%
% [x,val,fnrm,it,npcg]=sllsbox(b,A,xstart,options,data) returns
% the total number of conjugate gradient iterations used.
%
% [x,val,fnrm,it,npcg,ex]=sllsbox(b,A,xstart,options,data) returns
% the termination code.

%   Copyright (c) 1990-98 by The MathWorks, Inc.
%   $Revision: 1.11 $  $Date: 1998/09/05 15:23:16 $

%   INITIALIZATION
if nargin < 1, error('slls requires at least 1 argument'), end
%if nargin <2, A = []; end
if nargin < 5, xstart = []; end
if nargin < 6, verb = 1 ; end 
if nargin < 7, options = []; end
data = []; 

pcflags = optimget(options,'PrecondBandWidth') ;
tol = optimget(options,'tolfun') ;
itb = optimget(options,'maxiter') ;
pcgtol = optimget(options,'TolPCG');
active_tol = optimget(options,'ActiveConstrTol',sqrt(eps));

% These will be user-settable later.
pcmtx = optimget(options,'Preconditioner','aprecon') ;
mtxmpy = optimget(options,'JacobMult','atamult') ;
showstat = optimget(options,'showstatus','off');
switch showstat
case 'iter'
   showstat = 2;
case {'none','off'}
   showstat = 0;
case 'final'
   showstat = 1;
case 'iterplus'
   showstat = 3;
otherwise
   showstat = 1;
end

c = feval(mtxmpy,-b,A,data,-1); 
n = length(c); it = 1; cvec = c; nbnds = 1;
% In case the defaults were gathered from calling: optimset('lsqlin'):
numberOfVariables = n;
kmax = optimget(options,'MaxPCGIter', max(1,floor(n/2))) ;
if ischar(kmax)
   kmax = eval(kmax);
end

% get options out
typx = optimget(options,'typicalx') ;
% In case the defaults were gathered from calling: optimset('quadprog'):
numberOfVariables = n;
if ischar(typx)
   typx = eval(typx);
end

if n == 0, error('n must be positive'), end
if isempty(lb), lb = -inf*ones(n,1); end
if isempty(ub), ub = inf*ones(n,1); end
arg = (ub >= 1e10); arg2 = (lb <= -1e10);
ub(arg) = inf*ones(length(arg(arg>0)),1);
lb(arg2) = -inf*ones(length(arg2(arg2>0)),1);
if any(ub == lb) 
   errmsg=sprintf('%s\n%s',...
      'Equal upper and lower bounds not permitted in this large-scale method.',...
      'Use equality constraints and the medium-scale method instead.');
   error(errmsg)
elseif min(ub-lb) <= 0
   error('Inconsistent bounds.')
end
lvec = lb; uvec = ub;
if isempty(xstart), xstart = startx(ub,lb); end
if min(min(ub-xstart),min(xstart-lb)) < 0, xstart = startx(ub,lb); end
if isempty(typx), typx = ones(n,1); end
if isempty(verb), verb = 1; end

ds = ones(n,1); DS = sparse(1:n,1:n,ds);

%   SHIFT AND SCALE
[xstart,lb,ub,ds,DS,c] = shiftsc(xstart,lb,ub,typx,mtxmpy,data,cvec,A);

%   MORE INITIALIZATIONS
pcgit = 0; tol1 = tol; tol2 = sqrt(tol1);
dellow = 1.; delup = 10^3; npcg = 0; digits = inf; ex = 0; v = zeros(n,1);
dv = ones(n,1); del = 10*eps; posdef = 0;x = xstart; y = x;sigma = ones(n,1);
oval = inf; [val,g] = fquad(x,c,A,mtxmpy,data,DS);
vval(it,1) = val; vpcg(1,1) = 0; vpos(1,1) = 1;

if ((ub == inf*ones(n,1)) & (lb == -inf*ones(n,1)))
   nbnds = 0;
end
if showstat > 1,
   figtr=display1('init',itb,tol,showstat,nbnds,x,g,lb,ub);
end

%   MAIN LOOP: GENERATE FEAS. SEQ.  x(it) S.T. q(x(it)) IS DECREASING.
while ~ex
   
   %     Interactive exit?
   figtr = findobj('type','figure','Name','Progress Information') ;
   if ~isempty(figtr)
      lsotframe = findobj(figtr,'type','uicontrol',...
         'Userdata','LSOT frame') ;
      if get(lsotframe,'Value'), 
         ex = 10 % New exiting condition  
         if verb > 0
            display('Exiting per request.')
         end
         
      end ;
   end ;
   %
   %     Update and display
   [v,dv] = definev(g,x,lb,ub);
   fnrm = norm(v.*g,inf); 
   vfnrm(it,1) = fnrm;
   r = abs(min(ub-x,x-lb)); 
   degen = min(r+abs(g));
   vdeg(it,1) = min(degen,1);
   if ((ub == inf*ones(n,1)) & (lb == -inf*ones(n,1)))
      degen = -1;
   end
   bndfeas = min(min(x-lb,ub-x));
   if showstat > 1
      display1('progress',it,fnrm,val,pcgit,npcg,...
         degen,bndfeas,showstat,nbnds,x,g,lb,ub,figtr);
   end
   %
   %     TEST FOR CONVERGENCE
   diff = abs(oval-val); vdiff(it,1) = diff/(1+abs(oval));
   if it > 1, digits = (prev_diff)/max(diff,eps); end
   prev_diff = diff; oval = val; vflops(it,1) = flops;
   if diff < tol1*(1+abs(oval)),
      ex = 1;
      if verb > 0
         disp('Optimization terminated successfully:')
         disp(' Relative function value changing by less than OPTIONS.TolFun');
      end
   elseif (  (diff < tol2*(1+abs(oval))) & (digits < 3.5)),
      ex = 2;
      if verb > 0
         disp('Optimization terminated successfully:')
         disp(' Relative function value changing by less than sqrt(OPTIONS.TolFun),');
         disp(' and rate of progress is slow.')
      end
      
      
   elseif ((fnrm < tol1) & (posdef ==1)), 
      ex = 3;
      if verb > 0
         disp('Optimization terminated successfully:')
         disp(' First order optimality with optimality gradient norm less than OPTIONS.TolFun');
      end
      
      
   end
   %
   if ~ex
      %       DETERMINE THE SEARCH DIRECTION
      dd = abs(v); D = sparse(1:n,1:n,full(sqrt(dd).*sigma));
      grad = D*g; normg = norm(grad); 
      delta = max(dellow,norm(v)); delta = min(delta,delup);
      vdelta(it,1) = delta;
      [s,posdef,pcgit] = drqpbox(D,DS,grad,delta,g,dv,mtxmpy,data,...
         pcmtx,pcflags,pcgtol,A,1,kmax);
      npcg = npcg + pcgit; vpos(it+1,1) = posdef; vpcg(it+1,1) = pcgit;
      
      %       DO A REFLECTIVE (BISECTION) LINE SEARCH. UPDATE x,y,sigma.
      strg= s'*(sigma.*g); ox = x;  osig = sigma; ostrg = strg;
      if strg >= 0, 
         ex = 5;
         if verb > 0
            disp('Optimization terminated:')
            disp(' ill-conditioning prevents further optimization.');
         end
      else
         [x,sigma,alpha] = biqpbox(s,c,ostrg,ox,y,osig,lb,ub,oval,posdef,...
            normg,DS,mtxmpy,data,A);
         y = y + alpha*s; 
         
         %          PERTURB x AND y ?
         [pert,x,y] = perturb(x,lb,ub,del,y,sigma);
         
         %          EVALUATE NEW FUNCTION VALUE, GRADIENT. 
         it=it+1; [val,g] = fquad(x,c,A,mtxmpy,data,DS); 
         vval(it,1) = val;
      end
   end
   if it >= itb, 
      it=it -1; 
      ex=4;
      if verb > 0
         disp('Maximum number of function evaluations exceeded;')
         disp('   increase OPTIONS.MaxIter.')
      end     
   end
end

%   RESCALE, UNSHIFT, AND EXIT.
x = unshsca(x,lvec,uvec,DS);

residual = b-feval(mtxmpy,x,A,data,1);; % b-A*x
residnorm = sum(residual.*residual);


g = feval(mtxmpy,-residual,A,data,-1); % A'*(A*x-b) so -residual
g = full(g);
LAMBDA.lower = zeros(length(lvec),1);
LAMBDA.upper = zeros(length(uvec),1);
active_tol = sqrt(eps);
argl = logical(abs(x-lvec) < active_tol);
argu = logical(abs(x-uvec) < active_tol);

LAMBDA.lower(argl) = (g(argl));
LAMBDA.upper(argu) = -(g(argu));
LAMBDA.ineqlin = []; 
LAMBDA.eqlin = [];


if showstat > 1,  
   display1('final',figtr); 
end
if showstat, xplot(it,vval,vfnrm,vflops,vpos,vdeg,vpcg); end




