function [Q,R]=ModifiedGramSchmidtTwice(A);
% MODIFIEDGRAMSCHMIDTTWICE modified Gram-Schmidt with reorthogonalization
%   [Q,R]=ModifiedGramSchmidtTwice(A); applies the modified
%   Gram-Schmidt procedure with reorthogonalization to the columns in
%   matrix A

[m,n]=size(A); R=zeros(n);
Q=A;
for k=1:n,
  for t=1:2                                       % reorthogonalize 
    for i=1:k-1,
      V=Q(:,i)'*Q(:,k);                           % projections 
      Q(:,k)=Q(:,k)-V*Q(:,i); 
      R(i,k)=R(i,k)+V;
    end
  end
  R(k,k)=norm(Q(:,k)); Q(:,k)=Q(:,k)/R(k,k);
end
