function Q = gram_schmidt ( X )
% Q = gram_schmidt ( X ) computes the Gram-Schmidt
% orthogonalization of the columns of X regarded
% as vectors.  The resulting orthonormal set
% is returned as the columns of Q

nq=0;  % nq will be the number of columns of Q
for k=(1:size(X,2))
  for l=1:nq
    rlk=Q(:,l)'*X(:,k);
    X(:,k)=X(:,k)-rlk*Q(:,l);
  end
  rkk=norm(X(:,k));
  if rkk>0
    nq=nq+1;
    Q(:,nq)=X(:,k)/rkk;
  end
end