function [deltax,deltay]= solvlss(x,y,RHS,muu,...
                            noiter,nnzQ,toler,indsprec)
%%        [deltax,deltay]= solvlss(x,y,RHS,precond,muu,...
%%                          noiter,nnzQ,toler,indsprec)
%%
%  Solve the least squares solution for the Gauss-Newton Direction.
%  Use the lsqr routing of Paige-Saunders.
%
%  Input:
%        x,y,RHS    -   data from the current iteration
%                       x=usvec(X), y is dual variable, RHS=-Fmu(:)
%        muu                        contains current value of muu
%                                   to use for accuracy in lsqr
%        noiter      -  iteration number in calling program mcsdp
%                    -  used to help in determining accuracy for lsqr
%        nnzQ        -  number of nonzeros in Q, also
%                    -  used to help in determining accuracy for lsqr
%        toler       -  requested tolerance for relative 
%                         duality gap in     mcsdp.m
%       indsprec     -  indices needed for diag. preconditioner calculation
%
%  Output:
%           deltax;deltay  - search direction  

global Xu Q q2kron indsu ;


n=size(Q,1);
n2=n^2;
tn=(n*(n+1)/2);
tnm1=n*(n-1)/2;
%%%%%%%%%%%%%%%%%%%%%%5
%%%%%%%%%%%%%%%%%%%%%%5
%%prepare for MS lss iterative solve
damp=0;
 atol   = max(1.0e-1*toler,...
           (1.0e-6)*(min(muu,1))^(1/3)  );  
 btol   = max(1.0e-1*toler,...
           (1.0e-6)*(min(muu,1))^(1/3)  ); 
conlim = 1.0e+10;
itnlim = 10*tn;
show   = 0;
m=n^2;

%xy=[x;y];   % no longer needed to pass these
xy=[];


   timeprec=cputime;
   X=zeros(n);
   X(indsu)=x;
   X=X+X'+eye(n);
   X2=X.*X;
   x2=sum(X2,2);
   y2=y.*y;


   xy=[xy
       sqrt(y2(indsprec(:,1))+ y2(indsprec(:,2))+q2kron)/sqrt(2);
       sqrt(x2)];
   clear y2 x2 X
timeprec=cputime-timeprec;
     
%%%%%%%%%%%%%%%%%%%

%precond=[];         % was indicator for precond. - use now for Qy
Qy=-Q;
Qy(1:n+1:end)=y;         % replaces precond, diag(y)-Q and sparse
Xu=usMat(x,n);
%%[ deltas, istop, itn, rnorm, xnorm, var ] =  ...  %  for var variable
  [ deltas, istop, itn, rnorm, xnorm      ] =  ...
  MSlsqr( n^2, tn, 'MSFprecondfun', Qy, xy, RHS, damp, atol, btol, conlim, itnlim, show );



tnm1=n*(n-1)/2;
  %deltas=deltas./[xy(n+tnm1+1:length(xy))];
  deltas=deltas./xy;
deltax=deltas(1:tnm1);deltay=deltas(tnm1+1:end);