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