function [eigvec,steps]=EIGLRLag(c,s,f,d,cutoff,k)

% [eigvec,steps]=EIGLRLag(c,s,f,d,cutoff,k) computes eigenvalues of a 
% symmetric positive definite diagonal-plus-semiseparable matrix A, 
% constructed by the Givens representation with Givens rotations
% c_1,...,c_(n-1) and s_1,...,s_(n-1), vector f, and diagonal d
%
% cutoff is the parameter for the check for the convergence. The default 
% value is 1e-14.
% 
% k is the number of eigenvalues we want, by default all eigenvalues are 
% computed
%
% eigvec is the vector with the computed eigenvalues, steps is the number
% of the required iterations
%
% The method uses LR method with Laguerre shifts. The algorithm is
% describbed in 
% 
% Bor Plestenjak, Ellen Van Camp, and Marc Van Barel:
% A Cholesky LR algorithm for the positive definite symmetric 
% diagonal-plus-semiseparable eigenproblem 
%

% Bor Plestenjak
% Department of Mathematics, University of Ljubljana
% Jadranska 19, 1000 Ljubljana, SLovenia
% e-mail: bor.plestenjak@fmf.uni-lj.si
% 10.03.2005

eigvec=[];    % computed eigenvalues
found=0;      % number of eigenvalues computed so far 
steps=0;      % number of iterations
n=length(f);  
if nargin<5,  cutoff=1e-14; end    % cutoff criterium for the convergence

% Initial test if the matrix is positive definite at all
[f2,d2,flag] = CholSSD(c,s,f,d);
if flag==0
    error('Matrix is not positive definite')
end

a=1; b=n;      % interval of eigenvalue indices we are computing
shift=0;       % shift for the next LR step
lastshift=0;   % last good shift (where A-rho*I was s.p.d.)
while found < n
    m=length(f);
    if b>2
        [c2,s2,f2,d2,newshift,lastshift,tmpsteps]=EIGLRPhase(c(1:b-1),s(1:b-1),f(1:b),d(1:b),shift,lastshift,cutoff);
        steps=steps+tmpsteps;
        [c2,s2,f2,d2,eigvec]=CutEig(c2,s2,f2,d2,eigvec,cutoff);
        neweig=length(eigvec)-found;
        found=found+neweig;
        if b<m
            c=[c2; c(b:end)]; s=[s2; s(b:end)]; f=[f2; f(b+1:end)]; d=[d2; d(b+1:end)];
        else
            c=c2; s=s2; f=f2; d=d2;
        end
        b=b-neweig;
        for k=1:b-1
            if abs(s(k))<cutoff
                b=k;
                break;
            end
        end
    else
        if b==1
            eigvec=[eigvec, f(1)+d(1)];
            neweig=length(eigvec)-found;
            found=found+neweig;
            if b<m
                c=c(2:end); s=s(2:end); f=f(2:end); d=d(2:end);
                b=b-1;
            end
        end
        if b==2
            [c2,s2,f2,d2,eigvec]=CutEig(c(1:1),s(1:1),f(1:2),d(1:2),eigvec,cutoff);
            neweig=length(eigvec)-found;
            found=found+neweig;
            if b<m
                c=c(3:end); s=s(3:end); f=f(3:end); d=d(3:end);
                b=b-2;
            end
        end
    end
    if b<=0
        b=length(f);
    end
    shift=newshift;
end

eigvec=eigvec(:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary method - EIGLRPhase
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [c,s,f,d,shift,lastgood,steps]=EIGLRPhase(c,s,f,d,shift,lastgood,cutoff)

% In this function a loop of the Cholesky LR steps is performed until one 
% s(i) can be reduced to zero

steps=0;      % number of LR steps
overshoot=0;  % the last shift was such that A-rho I was not s.p.d.
level=0;      % first we try to reduce the shift by factor 1-1e-6, if this is not good we go back to zero

while min(abs(s)) > cutoff
    if overshoot==0   
        [c2,s2,f2,d,flag,newshift]=LRSSDLag(c,s,f,d,shift);
        % if the shift is not good for the first step we take shift=0
        if flag==1
            lastgood=shift; % new good shift
            shift=shift+newshift*(1-1e-4); % for precaution we do not use full shift
            c=c2; s=s2; f=f2;
            level=0;
            if (newshift/shift)<1e-6
                overshoot=1;
            end
        else
            shift=lastgood;
            overshoot=1;
        end
    else
        [c2,s2,f2,d,flag]=LRSSDLag(c,s,f,d,shift);
        if (flag==0) % if even the last shift is not good enough we go back to 0
            level=level+1;
            if level==1
                shift=(1-1e-6)*shift;
            else
                shift=0;
            end
            overshoot=0;
        else
            c=c2; s=s2; f=f2;
        end
    end
    steps = steps+1;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary method - LRSSDLag
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [c2,s2,f2,d1,flag,deltax]=LRSSDLag(c,s,f,d1,shift)

% [c2,s2,f2,d1,flag,deltax]=LRSSD(c,s,f,d) computes one step of LR method 
% with a given shift for a symmetric positive definite diagonal plus 
% semiseparable matrix A constructed by the so called Givens representation
% with Givens rotations c_1,...,c_(n-1) and s_1,...,s_(n-1), vector f, 
% and diagonal d. The matrix A - shift*I must be s.p.d.
%
% The result is a DPSS matrix B = VV' where V=chol(A). It is
% B=giv2full(c2,s2,f2,d). 
% 
% The output flag is 0 when A - shift*I is not s.p.d. and 1 otherwise.
% deltax is the modification of the current approximation for the
% eigenvalues obtained by the Laguerre method

n  = length(f);
c2 = zeros(n-1,1);
s2 = zeros(n-1,1);
f2 = zeros(n,1);
de = zeros(n,1);
d  = d1-shift;
deltax = 0;

if n<3, flag=0; c2=c; s2=s; f2=f; return, end % matrix is 2x2 or 1x1

% Cholesky decomposition
[fn,dn,flag,g,y,cz,ss] = CholSSD(c,s,f,d);

if flag==0, return, end % matrix is not s.p.d. 

% Laguerre shift for the next step
if nargout>5
    [S1,S2] = LaguerreStep(c,s,fn,dn,y);
    deltax = real(n/(S1+sqrt((n-1)*(n*S2-S1^2))));  % Laguerre shift
end

% Computation of V^T*V
de(n) = fn(n)+dn(n);
for k = n-1:-1:1
    de(k) = fn(k)+c(k)*dn(k);
end

f2(n) = de(n)^2 - d(n);
q = 0;
for k = n-1:-1:1
    q = ss(k)*(q+de(k+1)^2);
    a = cz(k) + ss(k)*g(k) ;
    b = fn(k)*sqrt(q);
    if (b == 0),
       c2(k) = 1; s2(k) = 0; f2(k) = a;
    else
       if (abs(a) >= abs(b))
          t = b/a; r = sqrt(1 + t*t);
          c2(k) = 1/r;
          s2(k) = t*c2(k);
          f2(k) = a*r;
       else
          t = a/b; r = sqrt(1 + t*t);
          s2(k) = 1/r;
          c2(k) = t*s2(k);
          f2(k) = b*r;
       end
    end
end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary method - LaguerreStep
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [y1,y2]=LaguerreStep(c,s,f,d,h)

% [y1,y2]=LaguerreStep(c,s,f,d) computes y1=norm(W,'fro')^2 and 
% y2=norm(WW^T,'fro')^2 where
% W=inv(V) and V=tril(Giv(c,s,f,d))
% h is a vector of the diagonal elements of V

n = length(f);
cn = zeros(n-1,1);
sn = zeros(n-1,1);
ws = zeros(n-1,1);
csk = zeros(n-1,1);

% when c(k)=0 we perturb the rotation slightly to avoid a compuation of 0/0
for k = 1:n-1
    if c(k)==0
        c(k)=1e-20;
    end
end

% Givens rotations in W
% first for c(n-1), s(n-1)
x = c(n-1)*h(n);
y = s(n-1)*d(n-1);
if (y == 0),
   cn(n-1) = 1; sn(n-1) = 0;
else
   if (abs(x) >= abs(y))
      t = y/x; r = sqrt(1 + t*t);
      cn(n-1) = 1/r;
      sn(n-1) = t*cn(n-1);
   else
      t = x/y; r = sqrt(1 + t*t);
      sn(n-1) = 1/r;
      cn(n-1) = t*sn(n-1);
   end
end
% then for k=1,...,n-2
for k=n-2:-1:1
    csk(k) = c(k+1)*s(k);
    x  = c(k)*cn(k+1)*h(k+1);
    y = csk(k)*d(k);
    if (y == 0),
       cn(k) = 1; sn(k) = 0;
    else
       if (abs(x) >= abs(y))
          t = y/x; r = sqrt(1 + t*t);
          cn(k) = 1/r;
          sn(k) = t*cn(k);
       else
          t = x/y; r = sqrt(1 + t*t);
          sn(k) = 1/r;
          cn(k) = t*sn(k);
       end
    end
end

% diagonal and subdiagonal elements of W
wd = 1./h;
csk(n-1) = s(n-1);
for k=1:n-1
    ws(k) = -csk(k)*f(k)*wd(k)*wd(k+1);
end

% WW^T
q = 0;
y1 = 0;
y2 = 0;
for k=1:n-2
    beta = ws(k)/cn(k+1);
    md = cn(k)^2*q + wd(k)^2;
    ms = cn(k)*sn(k)*q + beta*wd(k);
    q = sn(k)^2*q + beta^2;
    y2 = y2 + md^2 + 2*ms^2;
    y1 = y1 + md;
end
beta = ws(n-1); 
md = cn(n-1)^2*q + wd(n-1)^2;
ms = cn(n-1)*sn(n-1)*q + beta*wd(n-1);
q = sn(n-1)^2*q + beta^2;
y2 = y2 + md^2 + 2*ms^2;
y1 = y1 + md;
md = q + wd(n)^2;
y2 = y2 + md^2;
y1 = y1 + md;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary method - CholSSD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [fn,dn,flag,g,y,cz,ss]=CholSSD(c,s,f,d)

% [fn,dn,g,flag]=CHOLSSD(c,s,f,d) computes Cholesky decomposition 
% of a symmetric diagonal plus semiseparable matrix A
% constructed by the so called Givens representation with Givens rotations
% c_1,...,c_(n-1) and s_1,...,s_(n-1), vector f, and diagonal d.
%
% It returns [c,s,fn,dn,g] such that V=tril(giv2full(c,s,fn,dn)) is 
% the Cholesky factor for A. 
%
% flag = 1 when A is spd and flag = 0 when Cholesky fails
% g is a vector g=fn.^2 that can be used in the computation of one LR step for the matrix A
% y is a vector of diagonal elements of Cholesky factor (c(k)fn(k)+dn(k))
% cz is a vector of intermediate results that are used in the LR algorithm
% ss is a vector of squares of s (also used in the LR algorithm)

n  = length(f);
fn = zeros(n,1);
dn = zeros(n,1);
cz = zeros(n,1);
ss = zeros(n-1,1);
g  = zeros(n,1);
y  = zeros(n,1);
flag  = 0;

% k=1 
z = f(1);
cz(1) = c(1)*f(1); 
tmp = cz(1)+d(1);
if tmp<=0, return, end
y(1) = sqrt(tmp);
fn(1) = z/y(1);
dn(1) = d(1)/y(1);
g(1)  = fn(1)^2; 
ss(1) = s(1)^2;
q  = ss(1)*g(1);

for k = 2:n-1
    z = f(k) -  c(k)*q;
    cz(k) = c(k)*z;
    tmp = d(k)+cz(k);
    if tmp<=0, return, end
    y(k) = sqrt(tmp);
    fn(k) = z/y(k);
    dn(k) = d(k)/y(k);
    g(k)  = fn(k)^2;
    ss(k) = s(k)^2;
    q = ss(k)*(q+g(k));
end

% k=n
z = f(n) - q;
cz(n) = z;
tmp = d(n)+cz(n);
if tmp<=0, return, end
y(n) = sqrt(tmp);
fn(n) = y(n);
dn(n) = 0;
%fn(n)= z/y(n);
%dn(n)=d(n)/y(n);
gn(n) = tmp;

flag = 1; % A is spd

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Auxiliary method - CutEig
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [c,s,d,dg,eigvec]=CutEig(c,s,d,dg,eigvec,cutoff)

% this function will cut off eigenvalues
% which are already converged in the first steps
% of the algorithm

%    Software of the MaSe - Group
%    mase@cs.kuleuven.ac.be
%    Revision Date: 16/06/2004
% Adapted Bor Plestenjak 23.12.2004

n=length(d);
l=n-1;
if n==2
  B=[c(1)*d(1)+dg(1),s(1)*d(1);s(1)*d(1),d(2)+dg(2)];
  eigvec=[eigvec, sort(eig(B))'];
else
  while (l>1) & ( abs(s(n-1)*d(n-1)) < cutoff*(abs(d(n))+abs(dg(n)) ))
    eigvec=[eigvec,d(n)+dg(n)];
    n=n-1;
    l=l-1;
    d=d(1:n);
    dg=dg(1:n);
    d(n)=c(n)*d(n);
    s=s(1:n-1);
    c=c(1:n-1);
  end
end
