% Test of three flavors of QR decomposition
% See Trefethen and Bau, Lecture 9

 m = 30;                          % m = number of rows

 ortherr1 = zeros(m,1);
 ortherr2 = zeros(m,1);
 ortherr3 = zeros(m,1);
 for n=1:m                       % n = number of columns
    [U,X] = qr(randn(m,n),0);    % create a random orthogonal matrix, U
    [V,X] = qr(randn(n),0);      % create a random orthogonal matrix, V
    S = diag(2.^-[1:n]);         % diagonal matrix with exponentially decaying entries
    A = U*S*V';                  % construct the test matrix 
    [Q1,R1] = cgs_qr(A);
    [Q2,R2] = mgs_qr(A);
    [Q3,R3] = slow_householder_qr(A);
    ortherr1(n) = norm(Q1'*Q1-eye(n));
    ortherr2(n) = norm(Q2'*Q2-eye(n));
    ortherr3(n) = norm(Q3(:,1:n)'*Q3(:,1:n)-eye(n));
 end
   
 figure(1), clf
 semilogy(1:m, ortherr1, 'r-','linewidth',2); hold on
 semilogy(1:m, ortherr2, 'b-','linewidth',2); 
 semilogy(1:m, ortherr3, 'k-','linewidth',2); 
 set(gca,'fontsize',14)
 xlabel('number of columns of A, n','fontsize',16)
 ylabel('|| Q''*Q - I ||','fontsize',16)
 title('Departure from Orthogonality in Computed Q Factors')
 legend('Classical Gram-Schmidt', 'Modified Gram-Schmidt', 'Householder',2)

 break

 figure(2)
 pcolor(log10(abs(Q1'*Q1-eye(m))))
 set(gca,'ydir','reverse')
 axis equal, axis([1 30 1 30])
 caxis([-16 0]), colorbar
 title('log_{10} |Q''Q-I|, Classical Gram-Schmidt','fontsize',16)

 figure(3)
 pcolor(log10(abs(Q2'*Q2-eye(m))))
 set(gca,'ydir','reverse')
 axis equal, axis([1 30 1 30])
 caxis([-16 0]), colorbar
 title('log_{10} |Q''Q-I|, Modified Gram-Schmidt','fontsize',16)

 figure(4)
 pcolor(log10(abs(Q3'*Q3-eye(m))))
 set(gca,'ydir','reverse')
 axis equal, axis([1 30 1 30])
 caxis([-16 0]), colorbar
 title('log_{10} |Q''Q-I|, Householder QR','fontsize',16)