 function  [V, R]=mgs(X)
%
%-----------------------------------------------------------
%  Modified Gram-Schmit on X with selective re-orthogonalization
%-----------------------------------------------------------
%
%  Input:
%
%      X      (n-by-k)
%
%  Output:
%
%      V      n-by-k  Orthogonalized vectors from columns of X
%      N      Number of Orthogonalized vectors returned
%
%                                 RCL 06/28/2009
 [n,k]=size(X); 
%
 V = zeros(n,k); R = zeros(k,k); 
 nrm2 = norm(X(:,1)); R(1,1)=nrm2; V(:,1)=X(:,1)/nrm2;
 for j=2:k,
   vh = X(:,j); nrm2 = norm(vh);
   tolorth=1e-4*nrm2;
%  by MGS
   for i=1:j-1,
         R(i,j)=( V(:,i) )' * vh;
         vh = vh - R(i,j)*V(:,i);
   end
   nrm2=norm(vh); 
%  >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
%  Perform re-orthogonalization once by MGS when
%  deemed necessary.
   if nrm2 <= tolorth
      for i=1:j-1,
         tmp=( V(:,i) )' * vh;
         vh = vh - tmp*V(:,i);
         R(i,j)=R(i,j)+tmp;
      end
      nrm2=norm(vh); 
   end
   V(:,j)=vh/nrm2; R(j,j)=nrm2;
 end
