eigifp.m

Frequently used algorithms for numerical analysis and signal processing

            temp= Wb(:,k)'*r;
        else
            temp= V(:,k)'*r;
        end
        r = r - temp*V(:,k);
    end
                        % reorthogonalization if m > 6 
    if( m > 6)  
        for k = 1:(i-1),        
            if ( ~ isempty(B) )
                temp= Wb(:,k)'*r;
            else
                temp= V(:,k)'*r;
            end
            r = r - temp*V(:,k);
            end
        end
    if (norm(r) == 0)
        m=i;
        break 
    end
                    % normalize and save new basis vector 
                    % as well as Bx and (A-lambda B)x 
    if ( (~isempty(B)) & isnumeric(B) )
        bx=B*r;
    elseif ( (~isempty(B)) & isstr(B) )
        bx=feval(B, r);
    else 
        bx=r; 
    end
    
    temp=bx'*r;
    if (temp <= 0)
        error('Error: B is not positive definite.');
    end
    temp = sqrt(temp);
    
    V(:,i) = r;
    if (isstr(A))
        r=feval(A, V(:,i))-lambda*bx;
    else
        r=A*V(:,i)-lambda*bx;
    end
    V(:,i) = V(:,i) / temp;
    r=r/temp;
    if ( ~ isempty(B) )
        Wb(:,i)=bx/temp;
    end
    Wr(:, i) = r;   
                    % deflation for converged e-vectors
    if ( nargin == 16 )
        count=size(BEVecs,2);
        for j=1:count
            sigma=-lambda_1(j) + lambda_max(count);
            r=r  + sigma*(BEVecs(:,j)'*V(:,i))*BEVecs(:,j);
        end 
    end
    Am(1:i,i)=V(:,1:i)'*r;
end 

                    % add the diff vector to the basis
                    % and complete the projection Am 
diffNorm=sqrt(bDiff'*diff);
for k = 1:m-1,                       
    if ( ~ isempty(B) )
        temp = Wb(:,k)'*diff;
            bDiff = bDiff - temp*Wb(:,k);
    else
        temp = V(:,k)'*diff;
    end
        diff = diff - temp*V(:,k);
    rDiff = rDiff - temp*Wr(:,k);
end
                    % reorthogonalization if m > 6 
if( m > 6) 
    for k = 1:m-1,                  
        if ( ~ isempty(B) )
            temp = Wb(:,k)'*diff;
                    bDiff = bDiff - temp*Wb(:,k);
        else
            temp = V(:,k)'*diff;
        end
            diff = diff - temp*V(:,k);
        rDiff = rDiff - temp*Wr(:,k);
    end
end 
if ( ~ isempty(B) )
    temp = bDiff'*diff;
        if (temp < 0)
        if ( isnumeric(B) )
            bx=B*x;
        else
            bx=feval(B, x);
        end
        temp=bx'*x;   
        if (temp < 0)       
            error('Error: B is not positive definite.');
        end
    end
        temp=sqrt(temp); 
else
    temp = norm(diff);
end

                    % check and add diff only if it's significant
if ( temp <= 1e-8*diffNorm | temp==0 )
        m=m-1;
elseif ( temp <= 1e-2*diffNorm )    % recompute (A-lambda B)diff if necessary
    V(:,m)=diff/temp; 
    if ( (~isempty(B)) & isnumeric(B) )
        bDiff=B*V(:,m);
        Wb(:,m)=bDiff;
    elseif ( (~isempty(B)) & isstr(B) )
        bDiff=feval(B, V(:,m));
        Wb(:,m)=bDiff;
    else 
        bDiff=V(:,m); 
    end
    if (isstr(A))
        rDiff=feval(A, V(:,m))-lambda*bDiff;
    else
        rDiff=A*V(:,m)-lambda*bDiff;
    end    
    Wr(:,m)=rDiff;
    if (nargin == 16)
        count=size(BEVecs,2);
        for j=1:count
            sigma=-lambda_1(j) + lambda_max(count);
            rDiff=rDiff+sigma*(BEVecs(:,j)'*V(:,m))*BEVecs(:,j);
        end
    end     
    Am(1:m,m)=V(:,1:m)'*rDiff;    
else     
    V(:,m)=diff/temp; 
    Wr(:,m)=rDiff/temp;
    if ( ~ isempty(B) ), Wb(:,m)=bDiff/temp;   end
    r=Wr(:,m);
    if (nargin == 16)
        count=size(BEVecs,2);
        for j=1:count
            sigma=-lambda_1(j) + lambda_max(count);
            r=r+sigma*(BEVecs(:,j)'*V(:,m))*BEVecs(:,j);
        end
    end     
    Am(1:m,m)=V(:,1:m)'*r;
end
Am=triu(Am) + triu(Am,1)';

                    % compute Ritz value and vector of projection
[U, D]=eig(Am(1:m,1:m));
[delta, Ieig]=sort(diag(D));    

U(:,Ieig(1))=U(:,Ieig(1))/norm(U(:,Ieig(1))); 
x=V(:,1:m)*U(:,Ieig(1));
if ( ~ isempty(B) ),  bx=Wb(:,1:m)*U(:,Ieig(1)); end    
r=Wr(:,1:m)*U(:,Ieig(1));
                    % post processing Ritz vector
if (nargin == 16)               
    for j=1:count
            temp=BEVecs(:,j)'*x;
            x=x-EVecs(:,j)*temp;
            if ( ~ isempty(B) ),
                    bx=bx-BEVecs(:,j)*temp;
            end
            r=r-(lambda_1(j)-lambda)*BEVecs(:,j)*temp;
    end
end
if ( isempty(B) ), bx=x; end
sigma=(x'*r)/(x'*bx);
lambda=lambda+sigma; 
r=r-sigma*bx; 
                    % update new diff and related vectors
                    % for the next iteration   
U(1,Ieig(1)) = -(U(2:m,Ieig(1))'*U(2:m,Ieig(1)))/U(1,Ieig(1)); 
diff=V(:,1:m)*U(1:m,Ieig(1));
if ( ~ isempty(B) ),  
        bDiff=Wb(:,1:m)*U(1:m,Ieig(1)); 
else
        bDiff=diff;
end 
rDiff=Wr(:,1:m)*U(1:m,Ieig(1)); 
rDiff=rDiff-sigma*bDiff;

temp=norm(x,1); 
res=norm(r,1)/temp; 
lam_max=lambda+delta(size(delta,1));     
tol0=tol(1) + abs(lambda)*tol(2);
if ( res < tol0 | m>10 )        % recompute bx and r if necessary 
    if ( (~isempty(B)) & isnumeric(B) )
        bx=B*x;  
    elseif ( (~isempty(B)) & isstr(B) )
        bx=feval(B, x);
    else
        bx=x; 
    end 
    if (isstr(A))
        r=feval(A, x);
    else
        r=A*x;
    end
    lambda = (x'*r)/(x'*bx); 
    r=r-lambda*bx;
    res=norm(r,1)/temp;
    tol0=tol(1) + abs(lambda)*tol(2);   
end

                    % use a new diff if converged. 
if ( res < tol0 )
        diff = V(:,1:m)*U(:,Ieig(2));
end


function [lambda, x, r, bx, diff, rDiff, bDiff, d2, lam_max, res] = ...
     lancb(A, B, lambda, x, r, bx, diff, rDiff, bDiff, m, tol, EVecs, BEVecs, lambda_1, lambda_max)
%  
%   generate othonormal basis V of Krylov subspace by m steps of lanczos
%   and then compute the Ritz value and Ritz vectors 
%
%   BEVecs = B*EVecs  with EVecs =  converged eigenvectors
%

                    % initialization and normalization 
n=size(x,1); 
V = zeros(n,m);
Wr = V; 
Am=zeros(m,m);
beta=norm(x);
V(:,1) = x/beta; 
r=r/beta; 
bx=bx/beta;
Wr(:,1) = r;
if ( ~ isempty(B) ), 
    Bm=zeros(m,m); 
    Wb=V; 
    Wb(:,1) = bx;
end
beta=0;
                    % loop for Lanczos iteration 
for i = 1:(m-2)                          
    r=r-beta*x;     
    x=V(:,i);   
    alpha=x'*r; 
    r=r-alpha*x;
                    % reorthogonalization if m > 6
    if( m > 6 ) 
        for k = 1:i,                
            temp= V(:,k)'*r;
            r = r - temp*V(:,k);
        end
    end
                    % construct projection      
    beta=norm(r);
    Am(i,i)=alpha; 
    Am(i+1, i)=beta;
    Am(i, i+1)=beta;    
    if ( ~ isempty(B) ), Bm(1:i,i)=V(:,1:i)'*bx; end
    
    if (beta == 0)           
             break;
    end 
        V(:,i+1)=r/beta; 
        
                    % generate new basis vector     
    if ( (~isempty(B)) & isnumeric(B) )
        bx=B*V(:,i+1);  
        Wb(:,i+1)=bx;
    elseif ( (~isempty(B)) & isstr(B) )
        bx=feval(B, V(:,i+1));
        Wb(:,i+1)=bx;
    else 
        bx=V(:,i+1);  
    end 
    if (isstr(A))
        r=feval(A, V(:,i+1))-lambda*bx;
    else
        r=A*V(:,i+1)-lambda*bx;
    end
    Wr(:, i+1) = r; 
                    % deflation for converged e-vectors
    if (nargin == 15)
        count=size(BEVecs,2);
        for j=1:count
            sigma=-lambda_1(j) + lambda_max(count);
            r=r+sigma*(BEVecs(:,j)'*V(:,i+1))*BEVecs(:,j);
        end
    end
end 
                    % complete projection for the last vec
if (beta == 0) 
    m=i+1; 
else        
    Am(m-1, m-1) = V(:,m-1)'*r;
    if ( ~ isempty(B) ), Bm(1:m-1,m-1)=V(:,1:m-1)'*bx; end
end 

                    % add the diff vec to the basis
                    % and complete projection 
diffNorm=norm(diff);
for k = 1:(m-1),                         
    temp= V(:,k)'*diff;
    diff = diff - temp*V(:,k);
    rDiff = rDiff - temp*Wr(:,k);
    if ( ~ isempty(B) ), bDiff=bDiff-temp*Wb(:,k);  end
end
if ( m > 6)             % reorthogonalization if m>6 
    for k = 1:(m-1),                         
        temp= V(:,k)'*diff;
        diff = diff - temp*V(:,k);
        rDiff = rDiff - temp*Wr(:,k);
        if (~isempty(B)), bDiff=bDiff-temp*Wb(:,k);  end
    end
end 
temp = norm(diff); 

                    % check and add diff only if it's significant
if ( temp <= 1e-8*diffNorm | temp==0 )
        m=m-1; 
    if ( ~ isempty(B) ), 
            Bm=triu(Bm) + triu(Bm,1)';
        end
elseif ( temp <= 1e-2*diffNorm )
    V(:,m)=diff/temp; 
        if ( (~isempty(B)) & isnumeric(B) )
        bDiff=B*V(:,m);
            Wb(:,m)=bDiff;
    elseif ( (~isempty(B)) & isstr(B) )
        bDiff=feval(B, V(:,m));
            Wb(:,m)=bDiff;
    else 
        bDiff=V(:,m); 
    end
    if (isstr(A))
        rDiff=feval(A, V(:,m))-lambda*bDiff;
    else
        rDiff=A*V(:,m)-lambda*bDiff;
    end    
    Wr(:,m)=rDiff;
    if (nargin == 15)
        count=size(BEVecs,2);
        for j=1:count
            sigma=-lambda_1(j) + lambda_max(count);
            rDiff=rDiff+sigma*(BEVecs(:,j)'*V(:,m))*BEVecs(:,j);
        end
    end     
    Am(1:m,m)=V(:,1:m)'*rDiff;    
    Am(m,1:m)=Am(1:m,m)'; 
    if ( ~ isempty(B) ),  
        Bm(1:m,m)=V(:,1:m)'*Wb(:,m);
        Bm=triu(Bm) + triu(Bm,1)'; 
    end
else    
    V(:,m)=diff/temp; 
    Wr(:,m)=rDiff/temp;
    if ( ~ isempty(B) ), Wb(:,m)=bDiff/temp;   end
    r=Wr(:,m);
    if (nargin == 15)
        count=size(BEVecs,2);
        for j=1:count
            sigma=-lambda_1(j) + lambda_max(count);
            r=r  + sigma*(BEVecs(:,j)'*V(:,m))*BEVecs(:,j);
        end
    end
        
    Am(1:m,m)=V(:,1:m)'*r;
    Am(m,1:m)=Am(1:m,m)'; 
    if ( ~ isempty(B) ),  
        Bm(1:m,m)=V(:,1:m)'*Wb(:,m);
        Bm=triu(Bm) + triu(Bm,1)'; 
    end
end

                    % compute Ritz value and vector and
                    % approx gap d2 used for error estimate 
if (~ isempty(B) )
    [U, D]=eig(Am(1:m,1:m),Bm(1:m,1:m),'chol');
    [delta, Ieig]=sort(diag(D));                    
    EVals2=Bisection(Am(1:m,1:m), 0.1); 
    EVals2=sort(EVals2);
    if ( EVals2(2) <= 0 ) 
        d2=-1; 
    elseif ( EVals2(1) > 0 )
        d2=eps;
    else 
        d2=max(abs(EVals2));
    end     
else
    [U, D]=eig(Am(1:m,1:m));
    [delta, Ieig]=sort(diag(D));                
    EVals2=delta(1:2);
    if ( EVals2(2) <= 0 ) 
        d2=-1; 
    elseif ( EVals2(1) > 0 )
        d2=eps;
    else 
        d2=max(abs(EVals2));
    end                     
end

U(:,Ieig(1))=U(:,Ieig(1))/norm(U(:,Ieig(1))); 

x=V(:,1:m)*U(:,Ieig(1));
if ( ~ isempty(B) ),  bx=Wb(:,1:m)*U(:,Ieig(1)); end    
r=Wr(:,1:m)*U(:,Ieig(1));

if (nargin == 15)                   % post processing Ritz vector
    for j=1:count
            temp=BEVecs(:,j)'*x;
            x=x-EVecs(:,j)*temp;
            if ( ~ isempty(B) ),
                    bx=bx-BEVecs(:,j)*temp;
            end
            r=r-(lambda_1(j)-lambda)*BEVecs(:,j)*temp;
    end
end
if ( isempty(B) ), bx=x; end
sigma=(x'*r)/(x'*bx);
lambda=lambda+sigma; 
r=r-sigma*bx;

                    % update new diff and related vectors
                    % for the next iteration  
U(1,Ieig(1)) = -(U(2:m,Ieig(1))'*U(2:m,Ieig(1)))/U(1,Ieig(1)); 
diff=V(:,1:m)*U(1:m,Ieig(1));
if ( ~ isempty(B) ),  
        bDiff=Wb(:,1:m)*U(1:m,Ieig(1)); 
else
        bDiff=diff;
end 
rDiff=Wr(:,1:m)*U(1:m,Ieig(1)); 
rDiff=rDiff-sigma*bDiff;

temp=norm(x,1); 
res=norm(r,1)/temp; 
lam_max=lambda+delta(size(delta,1));    
tol0=tol(1) + abs(lambda)*tol(2);
if ( res < tol0 | m> 10)        % recompute bx and r if necessary 
    if ( (~isempty(B)) & isnumeric(B) )
        bx=B*x;  
    elseif ( (~isempty(B)) & isstr(B) )
        bx=feval(B, x);
        else
            bx=x; 
    end 
    if (isstr(A))
        r=feval(A, x);
    else
        r=A*x;
    end
    lambda = (x'*r)/(x'*bx); 
    r=r-lambda*bx;
    res=norm(r,1)/temp;
    tol0=tol(1) + abs(lambda)*tol(2);   
end


if (res < tol0 )            % use new diff if converged 
        diff = V(:,1:m)*U(:,Ieig(2));
end


 
function [Lambda_1, EVecs] = ...
    ifree(A, B, n, m, tol, itermax, k, X, eta, L, usePrecon, normA, normB,outputYes);
%
%  The main function that carries out the (outer) iteration. 
% 
t=cputime;
                    % initialization