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