%
%  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 = 1000;
   n = input(' specify 0 < n < 500   ')
   k = input(' specify 0 < k < n  ')
   test = input(' 1 or 2 ')
   v = randn(n,1);
   if (test == 1),
      M = randn(n);
%     v = randn(n,1);
      A = [];
      for j = 1:n,
         v = v/norm(v);
         A = [A v];
         v = M*v;
      end
   else
%
%
% Recommended Values  rho = .1,  tau = .001, .1, 1, 10
%
% 

     rho = input('input value of rho ')
     tau = input('input value of tau ')

      D = diag(randn(n,1));
      [Q,R] = qr(randn(n));

      for j = 1:n, R(j,j) = 0; end
      U = R;
%
%  U is strict upper triangular.
%

      dmax = max(abs(diag(D)));
      D = D/dmax;
%
%  D is diagonal with spec. radius 1 
%


       R =  rho*D + tau*U;
       A = Q*R*Q';
   end %if
   s = svd(A);
   KondA = s(1)/s(n)
   
   [V1,H1,f1] = Arnoldi(A,k,v);
   r1 = diag(H1);
   
   [V2,H2,f2] = ArnoldiMGS(A,k,v);
   r2 = diag(H2);
   
   [V3,H3,f3] = ArnoldiC(A,k,v);
   r3 = diag(H3);
   
   z1 = sqrt(eps)*ones(n,1);
   z2 = eps*ones(n,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(H) for 3 Methods')
   hold
disp('strike key')
pause


%
%  Test orthgonality of V  
%
   
   I = eye(k);
   
   CGSorth = norm(V1'*V1 - I) 
   MGSorth = norm(V2'*V2 - I) 
   CGSForth = norm(V3'*V3 - I) 
%
%  Surface plot of log of orthogonality
%  measure for CGS.
%
   figure
   mesh(log(abs(V1'*V1 - I)));
   colorbar
   title('CGS Orthogonality -- Log Scale')
disp('strike key')
pause
%
%  Surface plot of log of orthogonality
%  measure for MGS.
%
   figure
   mesh(log(abs(V2'*V2 - I)));
   colorbar
   title('MGS Orthogonality -- Log Scale')
disp('strike key')
pause
%
%  Surface plot of log of orthogonality
%  measure for CGSF.    Uncomment below to see it.
%
   figure
   mesh(log(abs(V3'*V3 - I)));
   colorbar
   title('CGSF Orthogonality -- Log Scale')
disp('strike key')
pause

%
%   Test difference  AV - VH - fe_k'   
%
   
   ek = zeros(k,1); ek(k) = 1;
   cgs_Residual = norm(A*V1 - V1*H1 -f1*ek')
   mgs_Residual = norm(A*V2 - V2*H2 -f2*ek')
   cgsf_Residual = norm(A*V3 - V3*H3 -f3*ek')

   s1 = eig(H1);
   s2 = eig(H2);
   s3 = eig(H3);

%  Evals =  [abs(s1) abs(s2) abs(s3)]

   figure
   plot(real(s1),imag(s1),'ro')
   hold on
   plot(real(s2),imag(s2),'b*')
   plot(real(s3),imag(s3),'g+')
   legend('CGS','MGS','CGSf')
   hold off