% This script illustrates the superior stability
% of the modified Gram-Schmidt algorithm
% (see Trefethen/Bau p.65).
% It uses two functions clgs and mgs that implement 
% the classical and modified Gram-Schmidt 
% orthogonalization algorithms. 
% clgs/mgs NEED TO BE PROVIDED BY THE USER.
[U, X] = qr(randn(80));    % U random orthogonal matrix
[V, X] = qr(randn(80));    % V random orthogonal matrix
S = diag(2.^(-1:-1:-80));  % S diagonal matrix with
                                        % exponentially graded entries
A = U*S*V;                       % A matrix with these entries
                                         % as singular values

% Compute classical and modified Gram-Schmidt 
[QC, RC] = clgs(A);
[QM, RM] = mgs(A);

% Compute machine epsilon
eps = 1;
while (1+eps) > 1
    eps = eps/2;
end
eps = eps*2;

% Create plot of diagonal entries of R and singular values
[m, n] = size(A);
for j=1:n
    yc(j) = RC(j,j);
    ym(j) = RM(j,j);
    ysing(j) = 2^(-j);
    x(j) = j;
end
yeps = eps*ones(n);
yepshalf = sqrt(eps)*ones(n);

semilogy(x,yc,'or',x,ym,'xg',x,ysing,'--k',x,yeps,'--k',...
    x,yepshalf,'--k','LineWidth',2,'MarkerSize',10)
xlabel('j','FontSize',14);
ylabel('r_{jj}','FontSize',14,'rotation',0);
legend('Classical GS', 'Modified GS');
text(5,2*10^(-9),'(Machine-\epsilon)^{1/2}','FontSize',14);
text(5,5*10^(-17),'Machine-\epsilon','FontSize',14);
text(70,5*10^(-23),'2^{-j}','FontSize',14);