% script file: compare
%
% compare different orthogonalization scheme.
%
% Dependency
%    ./BlkLanCom.m  componentwise block tridiagonalizaion algorithm 
%    ./BlkLanNorm.m block tridiagonalizaion algorithm using norm
%                   calculation
%    ./LanTri.m     tridiagonalization
%    ./errchk.m     check errors
%    ./testinfo.m   output the test informatio
%
matIn = 'mat128.mat';           % the input matrix file
testOut = 'testrep.txt';        % the output test report file
matOut = '';                    % the output test data file
%
load(matIn);                    % load the matrix        
if (strcmp(testOut, '') ~= 1)
    fid = fopen(testOut, 'w');
else
    fid = 1;
end
n = size(A,1);
blkSize = [2;4;8;16;32];
%
% constants
global ALGS
ALGS = 3;
SCALAR = 1;     % testinfo.m assumes SCALAR must be the first.
NORM   = 2; 
COM    = 3; 
%
% data structure for test 
testData = struct('tag',       '', ...
   'matFile',   matIn, ...               % the matrix file                     
   'matSize',   n, ...                   % matrix size
   'blkSize',   0, ...                   % block size                     
   'blkTriTime',0.0, 'triTime', 0.0, ... % time for block tridiag and tridiag
   'blkTriVec', 0,   'triVec',  0, ...   % vectors selected for orthogonalization                   
   'orthErr',   0.0, 'err',     0.0);    % error of orthogonality and factorization 
%
info = repmat(testData,ALGS,length(blkSize)); 
% initialize the test data
for j = 1:length(blkSize)          
    info(SCALAR,j).tag = 'Single-vector:';
    info(NORM,j).tag   = 'Block Norm:   ';
    info(COM,j).tag    = 'Componentwise:';
end
%
% scalar Lanczos tridiagonalization
t = cputime;
[a,b,Q,nIter,scalarTriVec] = LanMPO(A,ones(n,1),n);
scalarTriTime = cputime - t;
% check errors
[scalarOrthErr, scalarErr] = errchk (A,nIter,a,b,Q);
%
% record the scalar Lanczos test for later comparision
for j = 1:length(blkSize);
    info(SCALAR,j).triVec = scalarTriVec;
    info(SCALAR,j).triTime = scalarTriTime;
    info(SCALAR,j).orthErr = scalarOrthErr;
    info(SCALAR,j).err = scalarErr;
end
clear a b Q scalarTriVec scalarTriTime scalarOrthErr scalarErr;
%  
for j = 1:length(blkSize);
    bs = blkSize(j);    
    % create a random starting matrix for block lanczos
    [S, V] = qr(rand(n,bs),0);
    clear V;
    nSteps = n/bs;
%       
    for k = 2:ALGS
        % block tridiagonalization
        info(k,j).blkSize = bs;    
        t = cputime;
        switch k
            case NORM 
                fprintf('\n\nNow is Block Norm scheme... ...');
                [M,B,Q,nIter,info(k,j).blkTriVec] = BlkLanNorm(A,S,nSteps);
            case COM    
                fprintf('\nNow is Componentwise scheme... ...');
                [M,B,Q,nIter,info(k,j).blkTriVec] = BlkLanCom(A,S,nSteps);
        end
        info(k,j).blkTriTime = cputime - t;
%           
        % tridiagonalization
        t = cputime;;
        [a,b,P,info(k,j).triVec] = LanTri(M,B,ones(nIter*bs,1));
        info(k,j).triTime = cputime - t;        
        % check errors
        [info(k,j).orthErr, info(k,j).err] = errchk (A,nIter*bs,a,b,P,Q);
%
        clear M B Q a b P;
    end %for k = 2:ALGS
%
    % record the test report    
    testinfo(fid,info(:,j));
    clear S;
    fprintf('\nBlock size = %d. DONE!!!', bs);
end % for j = 1:length(blkSize);
%
if (fid ~= 1)
    fclose(fid);
end
if (strcmp(matOut, '') ~= 1)
    save(matOut, 'info');
end
fprintf('\nDONE!!!');