disp(' ');
disp(' Lanczos_tridiagonal.m   ver 1.0   July 10, 2010');
disp(' By Tom Irvine  Email: tomirvine@aol.com');
disp('  ');
%
%  Reference: Humar page 631
%
clear A;
clear X;
clear V;
clear Xhat;
clear Yhat;
clear T;
clear alpha;
clear beta;
clear size;
%
disp(' The input matrix must be real, symmetric. ');
disp(' ');
disp(' Enter the name of the matrix ');
A=input(' ');
%
n=max(size(A));
%
  V=rand(1,n);
% V=ones(n,1)*0.5;
%    
Xhat=zeros(n,n);
Xhat(:,1)=V;
beta(1)=sqrt(Xhat(:,1)'*Xhat(:,1));
X(:,1)=Xhat(:,1)/beta(1);
%
AX=A*X(:,1);
%
alpha(1)=X(:,1)'*AX;
Xhat(:,2)=AX-alpha(1)*X(:,1);
%
for(i=2:n)
   beta(i)=sqrt(Xhat(:,i)'*Xhat(:,i)); 
   X(:,i)=Xhat(:,i)/beta(i);
   AX=A*X(:,i);
   alpha(i)=X(:,i)'*AX;
   if(i==n)
       break;
   end
%   
   Xhat(:,i+1)=AX-beta(i)*X(:,i-1)-alpha(i)*X(:,i);
%
   Yhat=zeros(n,n);
   Yhat(:,1)=X(:,1)/norm(X(:,1));
%
   for ig=2:i
       for j=ig:i
           X(:,j)=X(:,j)-(Yhat(:,ig-1)'*X(:,j))*Yhat(:,ig-1);
       end
       Yhat(:,ig)=X(:,ig)/norm(X(:,ig));
   end   
   clear X;
   X=Yhat;
%
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Modified Gram-Schmidt Orthogonalization
%
Yhat=zeros(n,n);
Yhat(:,1)=X(:,1)/norm(X(:,1));
%
for i=2:n
    for j=i:n
        X(:,j)=X(:,j)-(Yhat(:,i-1)'*X(:,j))*Yhat(:,i-1);
    end
    Yhat(:,i)=X(:,i)/norm(X(:,i));
end   
clear X;
X=Yhat;
%
disp(' ');
disp(' Lanczos Vectors ');
X
disp(' ');
disp(' Tridiagonal Matrix ');
T=X'*A*X;
T