lanczos1.m

工程计算MATLAB code to calculate the reorthogonalized sine tapers input: N = the length of the time se

function [X1] = lanczos1(X1, M1, P, s)
 
%%MATLAB code to implement the block Lanczos procedure
%%input:  X1 = [n-by-p] initial guess for p eigenvectors of A=P*M*P
%%        M1 = [1-by-n] the first row of the Toeplitz matrix M
%%        P  = [n-by-n] projection matrix P
%%        s  = block dimension of the Lanczos procedure
%%output: X1 = [n-by-(s+1)*p] matrix obtained by the block Lanczos
%%             procedure 
 
p=length(X1(1,:));
Y=P*X1;
Y=cycmult(M1, Y);
Y=P*Y;
Y=X1'*Y;
for j=1:s
        Y2=P*X1(:,(((j-1)*p+1):(j*p)));
        Y2=cycmult(M1, Y2);
        Y2=P*Y2;
        if (j==1)
            R=Y2-X1(:,(((j-1)*p+1):(j*p)))*Y;
        else
            R=Y2-X1(:,(((j-1)*p+1):(j*p)))*Y-X1(:,(((j-2)*p+1):((j-1)*p)))*B';
        end
        R=R-X1*(X1'*R);        
        [X2,B]=qr(R);
        X2=X2(:,(1:p));
        X2=X2-X1*(X1'*X2);
        B=B((1:p),(1:p));
        Y=P*X2;
        Y=cycmult(M1, Y);
        Y=P*Y;
        Y=X2'*Y;
        X1=[X1, X2];
end