%%***************************************************************************
%% steplength: compute xstep such that X + xstep*dX >= 0.
%%
%% [xstep] = steplength(blk,X,dX,Xchol,invXchol);
%%
%% SDPT3: version 3.1
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 16 Sep 2004
%%***************************************************************************

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

%%
for p = 1:size(blk,1)
pblk = blk(p,:);
numblk = length(pblk{2});
pblksize = sum(pblk{2});
if (any(isnan(dX{p})) | any(isinf(dX{p}))); xstep = 0; break; end;
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)));
[lam,delta,res] = lanczosfun(Xchol{p},-dX{p},maxit,tol);
%%
%% 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) = 1e12;
end
elseif strcmp(pblk{1},'q')
aa = qops(pblk,dX{p},dX{p},2);
bb = qops(pblk,dX{p},X{p},2);
cc = qops(pblk,X{p},X{p},2);
dd = bb.*bb - aa.*cc;
tmp = min(aa,bb);
idx = find(dd > 0 & tmp < 0);
steptmp = 1e12*ones(numblk,1);
if ~isempty(idx)
steptmp(idx) = -(bb(idx)+sqrt(dd(idx)))./aa(idx);
end
idx = find(abs(aa) < eps & bb < 0);
if ~isempty(idx)
steptmp(idx) = -cc(idx)./(2*bb(idx));
end
%%
%% also need first component to be non-negative
%%
ss = 1 + [0, cumsum(pblk{2})];
ss = ss(1:length(pblk{2}));
dX0 = dX{p}(ss);
X0 = X{p}(ss);
idx = find(dX0 < 0 & X0 > 0);
if ~isempty(idx)
steptmp(idx) = min(steptmp(idx),-X0(idx)./dX0(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) = 1e12;
end
elseif strcmp(pblk{1},'u')
xstep(p) = 1e12;
end
end
xstep = min(xstep);
%%***************************************************************************
%%***************************************************************************
%% lanczos: find the largest eigenvalue of
%% invXchol'*dX*invXchol via the lanczos iteration.
%%
%% [lam,delta] = lanczosfun(Xchol,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] = lanczosfun(Xchol,dX,maxit,tol,v)

if (norm(dX,'fro') < 1e-13)
lam = 0; delta = 0; res = 0;
return;
end
n = length(dX);
if (nargin < 5);
state = randn('state');
randn('state',0);
v = randn(n,1);
randn('state',state);
end
if (nargin < 4); tol = 1e-3; end
if (nargin < 3); maxit = 30; end
V = zeros(n,maxit+1); H = zeros(maxit+1,maxit);
v = v/norm(v);
V(:,1) = v;
if issparse(Xchol); Xcholtransp = Xchol'; end
%%
%% lanczos iteration.
%%
for k = 1:maxit
if issparse(Xchol)
w = dX*mextriangsp(Xcholtransp,v,1);
w = mextriangsp(Xchol,w,2);
else
w = dX*mextriang(Xchol,v,1);
w = mextriang(Xchol,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 issparse(Xchol)
tmp = dX*mextriangsp(Xcholtransp,z,1);
res = norm(mextriangsp(Xchol,tmp,2) -lam*z);
tmp = dX*mextriangsp(Xcholtransp,z2,1);
res2 = norm(mextriangsp(Xchol,tmp,2) -lam*z2);
else
tmp = dX*mextriang(Xchol,z,1);
res = norm(mextriang(Xchol,tmp,2) -lam*z);
tmp = dX*mextriang(Xchol,z2,1);
res2 = norm(mextriang(Xchol,tmp,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
%%***************************************************************************