function [Q,R,ierr] = QRcgs(A);
%
%  Input:   A  -  an m by n matrix with m .ge. n
%
%  Output:  Q,R -  The Classical 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;
   
   rho = norm(A(:,1)); if ( rho == 0), ierr = 1; return; end
   Q = [A(:,1)/rho]; R = [rho]; 
   for j = 2:n,
       r = Q'*A(:,j);
       q = A(:,j) - Q*r;
       rho = norm(q);  if ( rho == 0), ierr = j; return; end
       q = q/rho; 
       Q = [Q,q]; R = [R r ; zeros(1,j-1) rho];
   end