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

[m n] = size(A);
R = zeros(n,n);
if max(max(abs(A))) == 0
    Q = eye(m,n);  return
end
Q = A;
for i = 1:n
    r = norm(Q(:,i),2);
    R(i,i) = r;
    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
