function [Q,R] = gsqr(A)
%
%  Gram-Schmidt orthogonalization for a QR decomposition:
%  [Q,R] = gsqr(A) returns a QR decomposition of A.
%
%  ON ENTRY :
%    A         an n-by-k real matrix, n >= k.
%
%  ON RETURN :
%    Q         an orthogonal matrix: Q'Q = I, whose column span 
%              is the same as the column span of A;
%    R         an upper triangular matrix: Q*R = A.
%
%  EXAMPLE :
%    a = rand(4,3);    % a random 4-by-3 matrix
%    [q,r] = gsqr(a);
%    a - q*r           % check residual
%    q'*q              % check orthogonality
%
k = size(A,2);
Q = A;
R = zeros(k,k);
for j = 1:k
  p = A(:,j);
  for i = 1:j-1
    R(i,j) = Q(:,i)'*A(:,j);
    p = p - Q(:,i)*R(i,j);
  end;
  Q(:,j) = p/norm(p,2);
  R(j,j) = Q(:,j)'*A(:,j);
end;