function [Q,R] = mgs(A)
  [m,n]=size(A); Q=A; R=eye(n);
  for i=1:n
    R(i,i)=norm(Q(:,i));
    if (R(i,i)<eps) break; endif;
    Q(:,i)=Q(:,i)/R(i,i);
    for j=i+1:n
      R(i,j)=Q(:,i)'*A(:,j);
      Q(:,j)=Q(:,j)-R(i,j)*Q(:,i);
    end;
  end;
end