function [Q,R] = clgs(A);
% CLGS computes reduced QR factorization by CLASSICAL Gram-Schmidt.
%    This is Trefethen & Bau Algorithm 7.1.  UNSTABLE!  Do not use for
%    serious computation.  See MGS or built-in QR (= Householder).
%
%    [Q,R] = clgs(A)  produces Q with orthonormal columns and R upper
%    triangular so that A = QR.  For m x n and m >= n matrices A with full rank,
%    returned Q is m by n, R is n by n; i.e. produces reduced QR.

myeps = 2*eps; % use machine eps times 2

[m n] = size(A);
scal = max(max(abs(A)));
R = zeros(n,n);

R(1,1) = norm(A(:,1),2);
if abs(R(1,1)) / scal < myeps
    error('nearly-zero first column; A not full rank?')
end
Q(:,1) = A(:,1) / R(1,1);

for j = 2:n
    v = A(:,j);
    for i = 1:(j-1)
        R(i,j) = Q(:,i)' * A(:,j);
        v = v - R(i,j) * Q(:,i);
    end
    R(j,j) = norm(v,2);
    if abs(R(j,j)) / scal < myeps
        error('nearly-zero column; A not full rank?')
    end
    Q(:,j) = v / R(j,j);
end
