function [Q,R,ierr] = QRcgs(A);
%
%  Input:   A  -  an m by n matrix with m .ge. n
%
%  Output:  Q,R -  The Modified Gram Schmidt QR factorization
%
%           ierr = 0   successful
%                  j   A not full rank detected at step j 
%
%  Factors A = QR ,  where A is m by n
%  Q'Q = I_n, R is n by n upper triangular
%
%  Q is m by n
%
%  D.C. Sorensen 
%  26 Oct 00
%-----------------------------------------------------------------

   [m,n] = size(A);
   ierr = 0;
   R = zeros(n); Q = zeros(m,n);
   
   rho = norm(A(:,1)); if ( rho == 0), ierr = 1; return; end
   Q(:,1) = [A(:,1)/rho]; R(1,1) = rho; 
   for j = 2:n,
       q = A(:,j);
       for i = 1:j-1
          R(i,j) = Q(:,i)'*q;
          q = q - Q(:,i)*R(i,j);
       end
       rho = norm(q); if ( rho == 0), ierr = j; return; end
       Q(:,j) = q/rho; R(j,j) = rho;
   end