eigifp.m

Frequently used algorithms for numerical analysis and signal processing

conv=ones(itermax,1);
rate=ones(10,1);

Lambda_1=zeros(k,1);
Lambda_Max=zeros(k,1);
EVecs=zeros(n,k);
BEVecs=EVecs;
pertbound=0; 
shift=0; 
if ( isempty(L) )
    regenL=1;
    startL=0;
else
    regenL=0;
    startL=1;
end

                    % set initial random vector if not given
if (~ isempty(X))
    x=X(:,1);
    normX=norm(x);
    if (normX == 0 | size(x,1) ~= n)
        error('Input Error: Invalid input of the initial vector.');
    end
else
    x=rand(n,1)-0.5*ones(n,1);
    normX=norm(x);
end 
x=x/normX;

                    % loop for computing k eigenvalues
for l=1:k

        diff=zeros(n,1);        % initialization for each eigenvalue iteration
        rDiff=diff; 
        bDiff=diff;
    
    if ( (~isempty(B)) & isnumeric(B) )
        bx=B*x; 
    elseif ( (~isempty(B)) & isstr(B) )
        bx=feval(B, x);
    else 
        bx=x;  
    end 

    temp=x'*bx; 
    if (temp <= 0)
        error('Error: B is not positive definite.');
    end

    if (isstr(A))
        r=feval(A, x);
    else
        r=A*x;
    end
        
    lambda=(x'*r)/temp;
    r=r-lambda*bx;
    res=norm(r);
    conv(1)=res;
    initialConv=0; 
    tol0=tol(1) + abs(lambda)*tol(2);
    
    matvecCount=1;
    preconCount=0;
    changeCounter=-itermax;
    mValue=m;
    if (mValue < 1)
        mValue=min(n-1,2);
        if (startL == 1), mValue=1; end
        changeCounter=-10;
    end
    priorM= mValue; 
        priorRate=-eps;
        
        if (outputYes ~=0)
        fprintf(1,'\nComputing Eigenvalue %d:\n',l);
        fprintf(1,'  Iteration\tEigenvalue\tResidual\n');
        fprintf(1,'  ======================================================\n');
        fprintf(1,'  %d\t\t%E\t%E\n',1,lambda,res); 
    end
    
                    % iteration for computing l-th eigenvalue
    for iter=2:itermax
        
        if ( res < tol0 )   % if converged, set parameters for next eigen
            lambda_max = normA; 
            initialConv=1;  
            break;
        end
        projSize=mValue+2;  
                    % compute Ritz value and vector by Krylov projection   
                    % call lancb if no preconditioning 
                    % call parnob if using preconditioning
        if (l>1)
                    % with deflation for l > 1      
          if ( startL == 0 )
                    [lambda, x, r, bx, diff, rDiff, bDiff, d2, lambda_max, res] = ...
              lancb(A, B, lambda, x, r, bx, diff, rDiff, bDiff, projSize,  tol,...
                          EVecs(:,1:l-1), BEVecs(:,1:l-1), Lambda_1(1:l-1), Lambda_Max(1:l-1,:));
          else
                    preconCount=preconCount+mValue;
                    [lambda, x, r, bx, diff, rDiff, bDiff, lambda_max, res] = ...
              parnob(A, B, L, lambda, x, r, bx, diff, rDiff, bDiff, projSize, tol,...
                          EVecs(:,1:l-1), BEVecs(:,1:l-1), Lambda_1(1:l-1), Lambda_Max(1:l-1,:));
          end           
        else
                    % no deflation if l=1. 
          if ( startL == 0 )
                    [lambda, x, r, bx, diff, rDiff, bDiff, d2, lambda_max, res] = ...
              lancb(A, B, lambda, x, r, bx, diff, rDiff, bDiff, projSize, tol);             
          else
            preconCount=preconCount+mValue;
                    [lambda, x, r, bx, diff, rDiff, bDiff, lambda_max, res] = ...
              parnob(A, B, L, lambda, x, r, bx, diff, rDiff, bDiff, projSize, tol);
          end           
        end
             
        conv(iter)=res; 
        matvecCount=matvecCount+mValue;         
        if (outputYes ~= 0)
            fprintf(1,'  %d\t\t%E\t%E\n',iter,lambda,res);
        end
        
                    % update tolerance and check convergence
        tol0=tol(1) + abs(lambda)*tol(2);
        if ( res <= tol0 )
            break;
        end

                    % check on convergence rate and update mValue
        changeCounter=changeCounter+1;
        if ( changeCounter >= 19 )
                    rate=[rate(2:10); aveRate(conv, iter-changeCounter+4, iter)];
                    [mValue, priorM, priorRate, fixM]=updateM(rate,mValue,priorM,priorRate,n);
                    changeCounter=(changeCounter+fixM*itermax)*(1-fixM);  
            elseif ( changeCounter >= 10 )
                    rate=[rate(2:10); aveRate(conv, iter-changeCounter+4, iter)];
            end
 
                    % estimate error 
        if ( (startL==0) & (d2>0) )
            err=min((res^2 / d2), res)/(bx'*x); 
        elseif (startL==0)
            err =abs(lambda);
        end

                    % determine whether to switch to preconditioning            
        if ((startL==0) & (usePrecon ~= 0) & (err <= (0.01*abs(lambda))) & (err < (0.0001*max(abs(lambda),normA/normB)) ) )
            startL=1;
            mValue=1;
            changeCounter=0;
            pertbound=0.0002*max(abs(lambda),normA/normB);          
            if (outputYes ~= 0)
                fprintf(1,'   \n');
                fprintf(1,'    Computing ILU factorization: if it takes too long, try the options:\n'); 
                fprintf(1,'      opt.iluThresh, opt.Preconditioner, or opt.usePrecon. \n');         
                fprintf(1,'   \n');
            end
            shift=lambda-err;
            if ( ~isempty(B) )
                L=ildlte(A-shift*B,eta);
            else
                L=ildlte(A-shift*speye(n),eta);
            end
            if (outputYes ~= 0)
                fprintf(1,'    Done!  Starting preconditioned iterations: \n');
                fprintf(1,'      if convergence is not accelerated, try decreasing opt.iluThresh. \n');
                fprintf(1,'   \n');         
            end     
        end
                        % determine whether to switch off preconditioning
        if ((regenL ~= 0) & (startL==1) & (changeCounter>15) & ( abs(shift-lambda) >= min(0.02*abs(shift), pertbound)))
            if ( sum(rate) > -0.2)  
                startL=0;
                changeCounter=-itermax;
                mValue=m;
                if (mValue < 1)
                    mValue=2;
                    changeCounter=0;
                end         
            end 
        end
    end

                        % store eigenpairs and others
    Lambda_1(l)=lambda;
    Lambda_Max(l)=lambda_max;
    temp=sqrt(bx'*x);
    EVecs(:,l)=x/temp; 
    BEVecs(:,l)=bx/temp;

                        % warn if not converged                         
    if ( res >= tol0 )
        if (isnumeric(A) & nnz(A-A') ~= 0)
            error('Input Error: Matrix A is not symmetric');
        end
        if ( (~isempty(B)) & isnumeric(B) & nnz(B-B') ~= 0)
            error('Input Error: Matrix B is not symmetric');
        end
        fprintf(1,'\n');
        fprintf(1,'Warning: Eigenvalue %d not converged to the tolerance within max iteration\n',l);
        fprintf(1,'         Residual = %E , the set tolerance = %E \n',res, tol0);
        fprintf(1,'         Try setting opt.innerit, a larger tolerance, and/or opt.maxit\n');
    end
    if (outputYes ~= 0)
        fprintf(1,'  ------------\n');
        fprintf(1,'  Eigenvalue %d converged. \n',l);       
        fprintf(1,'  # of multiplications by A (and B):      %d. \n',1+matvecCount);        
        if (preconCount > 0)
            fprintf(1,'  # of multiplications by preconditioner: %d.\n',  preconCount);     
        end 
    end

                        %  set next initial vector  
    if ( l<size(X,2) )
        x=X(:,l+1);
        normX=norm(x);
        if (normX == 0)
            fprintf(1,'Input Error: Invalid input of initial vector.\n');
        end
        x=diff;
    elseif (initialConv==1) 
        x=rand(n,1)-0.5*ones(n,1);  
    else
        x=diff;
    end
    x=x-EVecs(:,1:l)*(BEVecs(:,1:l)'*x);
    normX=norm(x);
    if (normX == 0)
        if (l<k) 
            error('Error: Fail to form an initial vector.');
        end 
        return;
    end
    x=x/normX;              

end
if (outputYes ~= 0)
    fprintf(1,'\n-----\n');
    fprintf(1,'  CPU Time:%f.\n',cputime-t);
end

function rate=aveRate(conv, k, iter)
%
% compute average linear convergence rate of conv over steps k+1 to iter 
%
rate=0; 
if ( iter-k < 2 | k<1 ) 
    return;
end

y=log10(conv((k+1):iter));
xAve=(iter-k+1)/2; 
xyAve=((1:iter-k)*y)/(iter-k)-xAve*sum(y)/(iter-k); 
xAve=((iter-k)^2-1)/12;
rate=xyAve/xAve;

        
  
  
function L=ildlte(A, eta)
%
%ILDLT   Threshold incomplete LDL^T  using luinc.m or ilu 
%

n=size(A, 1);
    
                    % compute L U factors 
options.thresh=0;
options.udiag=1;
if (eta ==2)
    options.milu=1; 
    [L, U] = luinc(A, options);
elseif (eta==1) 
    [L, U] = ilu(A);
else 
    options.droptol=eta;
    [L, U] = luinc(A, options);
end 

                    % scale diagonals to get L 
d = diag (U);
d = sqrt(abs(d));  
for k = 1:n, 
    if ( d(k) < 1e-8) 
        d(k) = 1e-8;
    end
end 
L = L * spdiags(d, 0, n, n);


function [L, U] = ilu(A)
%
% ILU   Incomplete LU factorization with no fill-in.
%   No pivoting is used 
%    

n = size(A, 1); 
M = A; 

for k = 1:(n-1), 
    if ( M(k, k) == 0)
        M(k, k) = 1e-8;
    elseif ( abs(M(k, k)) < 1e-8) 
        M(k, k) = 1e-8*sign(M(k, k));
    end 
    ind = find( M(:, k) ); 
    ind = ind(find(ind>k))';
    for i = ind,
            M(i, k) = M(i, k) / M(k, k); 
    end 

    ind_j = find( M(k, :) ); 
    ind_j = ind_j(find(ind_j>k))';

    for j = ind_j, 
            for i = ind,
                if (M(i, j) ~= 0) 
                    M(i, j) = M(i, j) - M(i, k)*M(k,j); 
                end 
            end 
    end 
end 

L=tril(M, -1)+speye(n); 
U=triu(M, 0);


function [mValue, priorM, priorRate, fixM]=updateM(rate, mValue, priorM, priorRate,n);          
%
%  Adaptive update of mValue: inner iteration 
%
fixM=0;
maxm=min(n-1, 128);

                    % update m when rate stagnates  
if ( (max(rate)-min(rate)) < 0.1*(-rate(10)) | min(rate) > 0 )
    k=2;                % increase m by k times, 
                    % use larger k if slower convergence 
        if ((rate(10) > -0.001) & 8*mValue <= maxm), 
            k=8; 
        elseif ((rate(10) > -0.01) & 4*mValue <= maxm), 
            k=4;
        end
                        % increase m by testing acceleration rate 
        incFlag=(rate(10)/priorRate)*(priorM/mValue); 
        if (incFlag > 1.05) 
            if(2*mValue > maxm )
                fixM=-1;
            else
                priorM=mValue; 
                priorRate=rate(10); 
                mValue=k*mValue;
                fixM=1;
            end  
        elseif ((rate(10) > -0.001 ) & 2*mValue <= maxm), 
            mValue=k*mValue;
            fixM=1;  
        elseif ((rate(10) > -0.01 ) & 2*mValue <= maxm), 
            mValue=k*mValue;
            fixM=1;   
        elseif (incFlag < 0.9)
            mValue=priorM;
            fixM=-1;   
        end
end


function lambda=Bisection(T,tol)
%
% estimate two lowest eigenvalues of tridiagonal T by Bisection 
%

a=diag(T);
b=diag(T,1);
k=2; 
lambda=zeros(k,1); 
xb=norm(T, 1)*1.0000001; 
xa=-xb; 
xc=0; 

p=0;
x1=[];
x2=[];
n1=[];
n2=[];
na=0; 
nc=Negcount(a,b,xc);

[p,x1,x2,n1,n2]=put(xa,na,xc,nc,p,x1,x2,n1,n2);
if( nc < k)
    nb=Negcount(a,b,xb);
    [p,x1,x2,n1,n2]=put(xc,nc,xb,nb,p,x1,x2,n1,n2);
end 
q=0;
while ( (p~=0) & (q < k) )
    [low,nlow,up,nup,p,x1,x2,n1,n2]=remove(p,x1,x2,n1,n2);
    if(nlow <= k) 
        if ( abs(up-low)<tol*abs(up) ) 
            ind=q+nup-nlow; 
            ind=min(ind, k); 
            q=q+1;
            lambda(q:ind)=((up+low)/2)*ones((ind-q+1),1);
        else
            mid=(low+up)/2;
            nmid=Negcount(a,b,mid);  
            if ( (nup>nmid) & (nmid <k) )
                [p,x1,x2,n1,n2]=put(mid,nmid,up,nup,p,x1,x2,n1,n2);
            end
            if nmid>nlow  
                [p,x1,x2,n1,n2]=put(low,nlow,mid,nmid,p,x1,x2,n1,n2);
            end
        end
    end 
end 


function [p,x1,x2,n1,n2]=put(xa,na,xb,nb,p,x1,x2,n1,n2)

p=p+1;
x1(p)=xa;
x2(p)=xb;
n1(p)=na;
n2(p)=nb;


function [xa,na,xb,nb,p,x1,x2,n1,n2]=remove(p,x1,x2,n1,n2)
xa=x1(p);
xb=x2(p);
na=n1(p);
nb=n2(p);
p=p-1;

 
function m=Negcount(a,b,z)
%
% NEGCOUNT: evaluates the number of eigenvalues that are smaller than z
%

m=0;
n=size(a,1);
f=size(b,1);

d(1)=a(1)-z;
if d(1)<0
    m=m+1;
elseif d(1)==0,   
        d(1)=eps^2; 
end
for i=2:n
    d(i)=a(i)-z-b(i-1)^2/d(i-1);
    if d(i)<0
        m=m+1;
        elseif d(i)==0
            d(i)=eps^2; 
    end
end


function x=minusA(x)
global Afunc123;

x=feval(Afunc123,x);
x=-x;