%
%  This script tests the three variants of Gram-Schmidt
%  QR factorization.
%
%    CGS - Classical Gram-Schmidt
%
%    MGS - Modified Gram-Schmidt
%
%    CGSF - Classical Gram-Schmidt with DGKS Ortogonality Fix

%
%   Compute the three different versions of QR
%
   m = 100;
   n = 80;
   M = randn(m);
   v = randn(m,1);
   A = [];
   for j = 1:n,
      v = v/norm(v);
      A = [A v];
      v = M*v;
   end
   s = svd(A);
   KondA = s(1)/s(n)
   
   [Q1,R1,ierr] = QRcgs(A);
   r1 = diag(R1);
   
   [Q2,R2,ierr] = QRmgs(A);
   r2 = diag(R2);
   
   [Q3,R3,ierr] = QRcgsF(A);
   r3 = diag(R3);
   
   z1 = sqrt(eps)*ones(80,1);
   z2 = eps*ones(80,1);
   
%
%  Graph results on semilog scale
%
   semilogy(r1,'ro');
   hold
   semilogy(r2,'+');
pause
   semilogy(r3,'g*');
   semilogy(z1,'-.');
   semilogy(z2,'-');
   semilogy(s,'k')
   legend('CGS','MGS','CGSF','Sqrt(Eps)', 'Eps','Svals',3)
   title('Comparison Diag(R) for 3 Methods')
   hold
pause


%
%  Test orthgonality of Q  
%
   
   I = eye(80);
   
   CGSorth = norm(Q1'*Q1 - I) 
   MGSorth = norm(Q2'*Q2 - I) 
   CGSForth = norm(Q3'*Q3 - I) 
%
%  Surface plot of log of orthogonality
%  measure for CGS.
%
   figure
   mesh(log(abs(Q1'*Q1 - I)));
   colorbar
   title('CGS Orthogonality -- Log Scale')
pause
%
%  Surface plot of log of orthogonality
%  measure for MGS.
%
   figure
   mesh(log(abs(Q2'*Q2 - I)));
   colorbar
   title('MGS Orthogonality -- Log Scale')
pause
%
%  Surface plot of log of orthogonality
%  measure for CGSF.    Uncomment below to see it.
%
   figure
   mesh(log(abs(Q3'*Q3 - I)));
   colorbar
   title('CGSF Orthogonality -- Log Scale')
pause

%
%   Test difference  A - QR   
%
   
   cgs_Residual = norm(A - Q1*R1)
   mgs_Residual = norm(A - Q2*R2)
   cgsf_Residual = norm(A - Q3*R3)