function [Q,R]=clgs(A);
% CLGS computes QR factors by CLASSICAL Gram-Schmidt.
%    Essentially Trefethen & Bau Algorithm 7.1.
%    See MGS for more accurate version and QR for Matlab built in version.
%
%    [Q,R] = clgs(A)  produces Q with orthonormal columns and R upper
%    triangular so that A = QR.  For m by n A with full rank (so m>=n).
%    Returned Q is m by n, R is n by n, that is, "economy size".
%ELB 10/03

myeps=2*eps; % use machine eps times 2
[m n]=size(A);
scal=max(max(abs(A)));
R=zeros(n,n);
if scal==0, Q=eye(m,n); return, end % case of zero matrix
r=norm(A(:,1),2);
R(1,1)=r;
if abs(r)/scal<myeps, error('nearly zero col; A not full rank?'), end
Q(:,1)=A(:,1)/r;
for j=2:n
    vo=A(:,j); v=vo;
    for i=1:(j-1)
        r=Q(:,i)'*vo;
        R(i,j)=r;
        v=v-r*Q(:,i);
    end
    r=norm(v,2);
    R(j,j)=r;
    if abs(r)/scal<myeps, error('nearly zero col; A not full rank?'), end
    Q(:,j)=v/r;
end