function [Q,R] = mgs(A)
%  Modified Gram-Schmidt for QR decomposition
%
   [m,n] = size(A);  % assuming m>n
   Q = A;
   R = zeros(n,n);
   for k = 1:n 
      R(k,k) = norm(Q(:,k));
      Q(:,k) = Q(:,k)/R(k,k);
      for j = k+1:n
	 R(k,j) = Q(:,k)'*Q(:,j);
	 Q(:,j) = Q(:,j) - R(k,j)*Q(:,k);
      end
   end