function [Q,R,z]=GramSchmidt(Q);
% GRAMSCHMIDT modified Gram-Schmidt with partial reorthogonalization
%   [Q,R,z]=GramSchmidt(Q); computes the QR decomposition of the
%   matrix Q using modified Gram-Schmidt with partial
%   reorthogonalization; z is a vector that counts the
%   reorthogonalization steps per column.
%   Translation of the ALGOL procedure in H. Rutishauser: "Algol 60"

z=[]; [m,n]=size(Q); R=zeros(n);
for k=1:n,
  t=norm(Q(:,k));
  reorth=1;
  u=0;                               % count reorthogonalizations
  while reorth,
    u=u+1;
    for i=1:k-1,
      s= Q(:,i)'*Q(:,k);
      R(i,k)=R(i,k)+s;
      Q(:,k)=Q(:,k)-s*Q(:,i);
    end
    tt=norm(Q(:,k));
    if tt>10*eps*t & tt<t/10,        % if length short reorthogonalize
      reorth=1; t=tt;           
    else
      reorth=0;
      if tt<10*eps*t, tt=0; end      % linearly dependent
    end
  end
  z=[z u];
  R(k,k)=tt;
  if tt*eps~=0, tt=1/tt; else tt=0; end
  Q(:,k)=Q(:,k)*tt;
end
