%SimulatedLanczos
%This example illustrates how the base Q of lanczos is succesive formed,
%by looking on how the eigenvectors occurenses in the base varies. The
%matrix X that is aquired is the Q matrix of lanczos to the diagonal
%egenvalue matrix D of the original matrix. That is, A=PDP'=P(XTX')P'.
%
%Please notice that this is a simulation of the Lanczos-algorithm, which means
%that surten effects as rounding off errors are not entirely represented in
%the code, since the code have access to the true eigenvectors and thus can
%perform perfect matrix multiplications. To simulate rounding off errors
%enter a noise value to 'instability'. There are however some other rounding off
%errors due to the procedure of removing components of alpha and beta, so
%the parameter 'reort' gives a chance to reort aswell.
%
%Exampel configurations:
%-Put gaps in the eigenvalues-spectra, and observe that eigenpairs close to
%the gap converges faster.
%-Use a startvector that is rich in components of desired eigenvectors.
%-Use a startvector that lacks components to vectors that are already
%known. This requires that some noise is inserted, to simulate properly.
%Perhaps in the order of 10^(-7).
%
%PARAMETERS
%eigenvalues: Enter the eigenvalues of the matrix.
%instability: The amount of rounding off error.
%reort: 1 if full reorthogonalisation, 0 i simple Lanczos.
%cStart: The initial distribution of eigenvectors. (will automatically be normalized)
%iterations: stop at iteration 'iterations'

eigenValues=1:100;
instability=10^(-7);
reort=1;
cStart=ones(1,length(eigenValues));
iterations=100;

%-------------------------------------------------
Size=length(eigenValues);
X=ones(Size,iterations);
X(:,1)=cStart/norm(cStart);
H=zeros(iterations);

for i=2:iterations
    X(:,i)=eigenValues'.*X(:,i-1);
        
    if reort
        reortUntil=1; %do full reorthogonalization
    else
        reortUntil=max([1 i-2]); %no reorthogonalization
    end
    
    for j=i-1:-1:reortUntil
        a=X(:,j)'*X(:,i);  %a=alpha if j=i-1, and a=beta if j=i-2
        X(:,i)=X(:,i)-a*X(:,j);
    end
    %simulate numeric instability
    X(:,i)=X(:,i)+instability*rand(Size,1);
    X(:,i)=X(:,i)/norm(X(:,i));
end

imagesc(abs(X()'))
%imagesc(log(abs(X()'))) %might give a better picture of the coefficient sizes
