function [EZout,EVout]=eigensolver(EXin,FLAG)

%-------------------------------------------------------------------------%
% This code computes eigenvalues and eigenvectors of a real symmetric     %
% matrix. Program returns number of iterations required for convergence   %
%and approximate number of flops.                                         %
%-------------------------------------------------------------------------%
% Usage: [lambda,U]=eigensolver(A,Option)                                 % 
%       A is real symmetric matrix                                        %
%       Option = 'ALL' - computes all eigenvalues and eigenvectors        %
%                'EVALUES' - computes all eigenvalues                     %
%                'L_EVECTOR' - computes largest eigenvalue and eigenvetor %
%                0 - computes largest eigenvalue                          %
%                                                                         %
%       lambda -- vector of eigenvalues                                   %
%       U      -- matrix of eigenvectors                                  %
%       ( A*U = U*diag(lambda) )                                          %
%-------------------------------------------------------------------------%
% For Options:                                                            %
%         'ALL' and 'EVALUES' -- algorithm uses shifted QR algorithm      %
%         with Wilkinson shift.                                           %
%                                                                         %
%         'L_EVECTOR' and 0   -- algorithm uses Rayleigh quotient         %
%         iterations.                                                     %
%          To ensure convergence to largest eigenvalue,                   %
%          Gershgorin's theorem is utilized to estimate largest           %
%          eigenvalue as starting shift.                                  %
%          To obtain an estimate for starting eigenvector a few inverse   %
%          iterations on random vector are done.                          %
%-------------------------------------------------------------------------%
%   Convergence of shifted QR slows if eigenvalues of matrix are close.   %
%                                                                         %
%   Convergence of RQI appears to be very dependent on initial vector as  %
%   well as on initial shift. Especially true if largest and second       %
%   largeest eigenvalues are close (A few steps of inverse iterations to    %
%   estimate initial vector for RQI seem to improve performance)          %                                                      
%-------------------------------------------------------------------------%

    thresh=1e-12;
    [n,m]=size(EXin);
    NN=4;  %# of iterations to prepare estimate eigenvector for RQI      %   
%-------------------------------------------------------------------------%
function [out,QQ]=tridiag(in,FLAG)
        [n,m]=size(in);
        aux=in;
        v=zeros(n,n);
        if strcmp(FLAG,'L_EVECTOR')||strcmp(FLAG,'ALL')
            QQ=eye(n);
        end
        
        for k=1:n-2
            Q=eye(n);
            e1=zeros(n-k,1);
            e1(1)=1;
            x=aux(k+1:n,k);
            nx=norm(x,2);
            v(k+1:n,k)=(sign(x(1))*nx*e1(1:n-k)+x(1:n-k));
            v(k+1:n,k)=v(k+1:n,k)/norm(v(k+1:n,k),2);
            aux(k+1:n,k:n)=aux(k+1:n,k:n)-2*v(k+1:n,k)*(v(k+1:n,k)'*aux(k+1:n,k:n));
            aux(1:n,k+1:n)=aux(1:n,k+1:n)-2*(aux(1:n,k+1:n)*v(k+1:n,k))*v(k+1:n,k)';
            if strcmp(FLAG,'L_EVECTOR')||strcmp(FLAG,'ALL')
                Q(k+1:n,k+1:n)=eye(n-k)-2*v(k+1:n,k)*v(k+1:n,k)';
                QQ=QQ*Q';
            else
                QQ=0;
            end
        end       
        out=aux;
    end 
%-------------------------------------------------------------------------%
    function [out,Q,k]=shQR(in,FLAG)
        [n,m]=size(in);
        aux=in;
        [aux,Q0]=tridiag(aux,FLAG);
        k=0;
        Q1=eye(n);
        while norm(diag(aux,1),inf)/norm(diag(aux),inf)>thresh
            check=norm(diag(aux,1),inf)/norm(diag(aux),inf);
            k=k+1;
            mu=wilkshift(aux);
            [U,R]=QRmod(aux-mu*eye(n));
            aux=R*U+mu*eye(n);
            if strcmp(FLAG,'ALL')
                Q1=Q1*U;
            end
        end
        out=diag(aux);
        if strcmp(FLAG,'ALL')
            Q=Q0*Q1;
        else
            Q='      eigenvector option disabled      ';
        end
    end
%-------------------------------------------------------------------------% Backward substitution
    function out=back(inM,rhs)
        [n,m]=size(inM);
        rhs(n)=rhs(n)/inM(n,n);
        for j=n-1:-1:1
            rhs(j)=(rhs(j)-inM(j,j+1:n)*rhs(j+1:n))/inM(j,j);
        end
        out=rhs;
    end
%-------------------------------------------------------------------------%Inverse Iterations
    function out=inviter(in,est)
        [n,m]=size(in); 
        u=randn(n,1);
        for j=1:NN
            [q,r]=QRmod(in-est*eye(n));
            v=q'*u;
            u=back(r,v);
            u=u/norm(u,2);
        end
        out=u;
    end
%-------------------------------------------------------------------------% Raileigh iterations
    function [out,vout,k]=rail(in,FLAG)
        k=0;
        [n,m]=size(in);
        [aux,Q0]=tridiag(in,FLAG);
        l1=linspace(0,0,n);
        for j=1:n
            l1(j)=norm(aux(1:n,j),1);
        end
        l=norm(l1,inf);
        lb=0; v=inviter(aux,l);
        while abs(l-lb)>thresh
                k=k+1;
                [Q,R]=QRmod(aux-l*eye(n));
                y=Q'*v;
                w=back(R,y);
                v=w/(norm(w,2));            
                lb=l;
                l=v'*aux*v;
        end
        out=l;
        if strcmp(FLAG,'L_EVECTOR')
            vout=Q0*v;
        else
            vout='      eigenvector option disabled      ';
        end
    end
%-------------------------------------------------------------------------%
    function out=wilkshift(in)
        [n,m]=size(in);
        d=0.5*(in(n-1,n-1)-in(n,n));
        out=in(n,n)+d-sign(d)*sqrt(d^2+in(n-1,n)*in(n,n-1));
    end
%-------------------------------------------------------------------------%
    function [c s]=givens(a,b)
        if b==0
            c=1;s=0;
        else
            if abs(b)>abs(a)
                r=-a/b;s=1/sqrt(1+r^2);c=s*r;
            else
                r=-b/a;c=1/sqrt(1+r^2);s=c*r;
            end
        end
    end
%-------------------------------------------------------------------------%
    function [QQ,R]=QRmod(A)
        [n,m]=size(A);
        QQ=eye(n);
        for j=1:n-1
            Q=eye(n);
            [c s]=givens(A(j,j),A(j+1,j));
            W=[c s;-s c];
            Q(j:j+1,j:j+1)=W;
            A(j:j+1,j:n)=W'*A(j:j+1,j:n);
            QQ=QQ*Q;
        end
        R=A;
    end
%-------------------------------------------------------------------------%
if strcmp(FLAG,'ALL')
    [EZout EVout k]=shQR(EXin,'ALL');
    Iterations_count=k
    flops=30*k*n^2
    %norm(EXin*EVout-EVout*diag(EZout))
elseif strcmp(FLAG,'EVALUES')
    [EZout EVout k]=shQR(EXin,0);
    Iterations_count=k
    flops=14*k*n^2
elseif  strcmp(FLAG,'L_EVECTOR')
    [EZout EVout k]=rail(EXin,'L_EVECTOR');
    Iterations_count=k
    flops=6*k*n^2
else     
    [EZout EVout k]=rail(EXin,0);
    Iterations_count=k
    flops=4*k*n^2
end
end