function [Q,R,ierr] = QRcgsF(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
%
%  Maintains orthogonality using DGKS fix
%
%  Q is m by n
%
%  D.C. Sorensen 
%  26 Oct 00
%  Modified 28 Sep 01  by DCS
%    included additional test 
%    for very pathological cases
%    of rank deficient matrices.
%
%-----------------------------------------------------------------

   [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;
%
%         Apply DGKS correction
%
          c = Q'*q;
          q = q - Q*c;
          r = r + c;

       rho = norm(q);  if ( rho == 0), ierr = j; return; end
       q = q/rho; 
%
%      The following test asks if A(:,j) is numerically in Range(Q).
%      If the test is passed then rho is indistinguishable from 0
%      when compared to norm(r) and the component q is meaningless
%
       if ( rho < 10*eps*norm(r)),
%
%          Construct a random q and orthogonalize against Q
%          Since rho is indistinguishable from 0, q*rho will be of
%          negligible size for any unit vector q
%      
           q = randn(m,1);
           c = Q'*q;
           q = q - Q*c;
              c = Q'*q;
              q = q - Q*c;
           q = q/norm(q);
       end
          
       Q = [Q,q]; R = [R r ; zeros(1,j-1) rho];

   end