%%***************************************************************************
%% steplength: compute xstep such that  X + xstep*dX >= 0.
%%
%% [xstep] = steplength(blk,X,dX,Xchol,invXchol,randnstate);
%%
%% SDPT3: version 3.0 
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%***************************************************************************

    function [xstep,Xchol,invXchol] = steplength(blk,X,dX,Xchol,invXchol,randnstate);

%%
    if (nargin == 3); randnstate = 0; end; 
    randn('state',randnstate); 
%%
    for p = 1:size(blk,1)
       pblk = blk(p,:); 
       numblk = length(pblk{2}); 
       pblksize = sum(pblk{2}); 
       if strcmp(pblk{1},'s') 
          if (max(pblk{2}) >= 200);  
             use_lanczos = 1; 
          else
             use_lanczos = 0; 
          end; 
          if (use_lanczos)
             tol = 1e-3;     
             maxit = max(min(pblksize,30),round(sqrt(pblksize))); 
             initvec = randn(pblksize,1);
             if issparse(Xchol{p}); 
                invXchol{p} = inv(Xchol{p});
     	        Xchol{p} = []; 
             end
             [lam,delta,res] = lanczos(Xchol{p},invXchol{p},dX{p},maxit,tol,initvec);
             %%
             %% Note: lam <= actual largest eigenvalue <= lam + delta.
             %%   
	     d = lam+delta; 
          else
             if isempty(invXchol{p}); 
                invXchol{p} = inv(Xchol{p}); 
             end
	     tmp = Prod2(pblk,dX{p},invXchol{p},0); 
             M = Prod2(pblk,invXchol{p}',tmp,1); 
             d = blkeig(pblk,-M); 
          end
          tmp = max(d) + 1e-15*max(abs(d)); 
          if (tmp > 0);  
             xstep(p) = 1/max(tmp);
          else;          
             xstep(p) = 1e6; 
          end;
       elseif strcmp(pblk{1},'q')
          aa = mexqops(pblk,dX{p},dX{p},2); 
          bb = mexqops(pblk,dX{p},X{p},2); 
          cc = mexqops(pblk,X{p},X{p},2);
          dd = bb.*bb - aa.*cc; 
          tmp = min(aa,bb); 
          idx = find(dd > 0 & tmp < 0); 
          steptmp = 1e6*ones(numblk,1); 
          if ~isempty(idx)
             steptmp(idx) = -(sqrt(dd(idx))+bb(idx))./aa(idx);       
          end
          idx = find(abs(aa) < eps & bb < 0); 
          if ~isempty(idx)
             steptmp(idx) = -cc(idx)./(2*bb(idx)); 
          end
          xstep(p) = min(steptmp); 
       elseif strcmp(pblk{1},'l')
          idx = find(dX{p} < 0); 
          if ~isempty(idx); 
             xstep(p) = min(-X{p}(idx)./dX{p}(idx));  
          else 
             xstep(p) = 1e6;
          end
       end
    end
    xstep = min(xstep); 
%%===========================================================================
%% lanczos: find the largest eigenvalue of 
%%          invXchol'*dX*invXchol via the lanczos iteration.
%%
%% [lam,delta] = lanczos(invXchol,dX,maxit,tol,v)
%%
%% lam:  an estimate of the largest eigenvalue.
%% lam2: an estimate of the second largest eigenvalue.
%% res:  residual norm of the largest eigen-pair.
%% res2: residual norm of the second largest eigen-pair.
%%
%%***************************************************************************

   function [lam,delta,res] = lanczos(Xchol,invXchol,dX,maxit,tol,v) 
                                            
   n = length(dX); 
   V = zeros(n,maxit+1); H = zeros(maxit+1,maxit); 
   v = v/norm(v); 
   V(:,1) = v; 
%%
%% lanczos iteration. 
%%
   for k = 1:maxit
      if ~isempty(invXchol) 
         w = dX*(invXchol*v); 
         w = -(w'*invXchol)'; 
      else         
         w = mexbsolve(Xchol,v,1); 
         w = -mexbsolve(Xchol,dX*w,2); 
      end
      wold = w;
      if (k > 1); 
         w = w - H(k,k-1)*V(:,k-1); 
      end;
      alp = w'*V(:,k); 
      w   = w - alp*V(:,k); 
      H(k,k) = alp; 
      %%
      %% one step of iterative refinement if necessary. 
      %%
      if (norm(w) <= 0.8*norm(wold));
         s = (w'*V(:,1:k))'; 
         w = w - V(:,1:k)*s;
         H(1:k,k) = H(1:k,k) + s;
      end; 
      nrm = norm(w); 
      v = w/nrm; 
      V(:,k+1) = v; 
      H(k+1,k) = nrm;  H(k,k+1) = nrm; 
      %%
      %% compute ritz pairs and test for convergence
      %%
      if (rem(k,5) == 0) | (k == maxit); 
         Hk = H(1:k,1:k); Hk = 0.5*(Hk+Hk'); 
         [Y,D] = eig(Hk); 
         eigH  = real(diag(D)); 
         [dummy,idx] = sort(eigH);
         res_est = abs(H(k+1,k)*Y(k,idx(k)));
         if (res_est <= 0.1*tol) | (k == maxit);
            lam = eigH(idx(k));  
            lam2 = eigH(idx(k-1)); 
            z   = V(:,1:k)*Y(:,idx(k));
            z2   = V(:,1:k)*Y(:,idx(k-1));
            if ~isempty(invXchol); 
               res = norm(-invXchol'*(dX*(invXchol*z)) -lam*z); 
               res2 = norm(-invXchol'*(dX*(invXchol*z2)) -lam*z2);   
            else
               res = norm(-mexbsolve(Xchol,dX*mexbsolve(Xchol,z,1),2) -lam*z); 
               res2 = norm(-mexbsolve(Xchol,dX*mexbsolve(Xchol,z2,1),2) -lam*z2);   
            end
            tmp = lam-lam2 - res2; 
            if (tmp > 0); beta = tmp; else; beta = eps; end;  
            delta = min(res,res^2/beta); 
            if (delta <= tol); break; end;
         end; 
      end; 
   end;
%%
%%***************************************************************************
