#
#   function [Y,nonnull]=grschm(X,G)
#
#   Gram-Schmidt orthonormalization of (the columns of) X
#   with respect to the scalar product with Gram matrix G.
#   This second argument is optional, if not entered it is
#   assumed to be the identity matrix, otherwise it must
#   be an (n,n) symmetric positive definite matrix (where
#   n is the length of the columns of X).
#
#   Zero-length columns (compared with epsilon) are set to
#   zero exactly and moved to the last positions. The number
#   of resulting nonnull columns is also returned.
#
#
function [Y,nonnull]=grschm(X,G)
[n,m]=size(X);
if nargin<2
     G=eye(n);
else
   [n1,n2]=size(G);
   if n1!=n|n2!=n
      disp('*** grmschm: wrong G');
      return;
   endif
endif
epsilon=1.e-4;
j=1;
nonnull=m;
while j<=nonnull
   j1=j-1;
   if j1>0
      for k=1:j1
         proj=X(:,j)'*G*X(:,k);
         X(:,j)=X(:,j)-proj*X(:,k);
      endfor
   endif
   ssq=X(:,j)'*G*X(:,j);
   if abs(ssq)<epsilon
      disp('*** grschm: zero-length column');
      X(:,j)=zeros(n,1);
      nonnull--;
      if j<nonnull
         J1=j:nonnull;
         J2=[(j+1):nonnull j];
         X(:,J1)=X(:,J2);
      endif
   else
      X(:,j)=X(:,j)/sqrt(ssq);
      j++;
   endif
endwhile
Y=X;
endfunction