function [Q,R] = mgs(A);
% MGS computes reduced QR decomposition

% ELB 3/10/09: better rank-deficiency checking; compare mgs2.m

[m n] = size(A);
tol = (max([m n]) + 10) * eps;
R = zeros(n,n);
scal = max(max(abs(A)));
if scal == 0
    Q = eye(m,n);
    return, end
Q = A;
for i = 1:n
    r = norm(Q(:,i),2);
    R(i,i) = r;
    if abs(r)/scal < tol
        error(['rank deficient: nearly zero internal column; i=' num2str(i)])
    end
    w = Q(:,i)/r;
    Q(:,i) = w;
    for j = i+1:n
        r = w'*Q(:,j);
        R(i,j) = r;
        Q(:,j) = Q(:,j)-r*w;
    end
end