# procedure gram_schmidt,  Maple V.4, packable
# Copyright: Erich Kaltofen

# Gram/Schmidt orthogonalization

LsqPkg['Gram_Schmidt'] := proc(A::Matrix)
# returns a matrix Q whose columns are pairwise orthogonal and
#    span the column space of A
# an optional second argument is the co-factor R such that A = QR
  local m,n, i,j,k,
        v,  # list of column in A
        u,  # list of vectors in orthogonal set
        norm_u_sq, # list of || . ||^2 of u's
        w, norm_sq, # temporary values
        Q, R;
	option `Copyright (c) 1997 by Erich Kaltofen`;
        description `optional second arg 'R' such that A = QR`;

  m := LinearAlgebra[RowDimension](A);
  n := LinearAlgebra[ColumnDimension](A);

  if nargs > 1 # name of right factor matrix
     then R := Matrix(n,n);
  fi;

  v := [ LinearAlgebra[Column](A, 1 .. n) ]; # extract columns of A
  u := []; norm_u_sq := []; k := 0; # current length of u

  for i from 1 to n do

      # orthogonalize v[i]
      w := v[i];
      for j from 1 to k do
         # pull w along projection of v[i] with u[j]
         w := w - LinearAlgebra[DotProduct](v[i], u[j], conjugate=false)/norm_u_sq[j] * u[j];
         
         if nargs > 1 then R[j,i] := LinearAlgebra[DotProduct](v[i], u[j], conjugate = false)/norm_u_sq[j]; fi;
      od;

      if nargs > 1 then
         for j from k+1 to n do R[j,i] := 0; od;
      fi;

      if printlevel >= 2 then
         print(`v[` || i || `]=`, v[i]);
         print(`u[` || i || `]=`, w);
         for j from 1 to k do
             print(`<u[` || j || `], u[` || i || `]>=`, LinearAlgebra[DotProduct](u[j], w, conjugate = false));
         od;
      fi;

      # compute norm of orthogonalized v[i]
      norm_sq := LinearAlgebra[DotProduct](w, w, conjugate = false);
      if printlevel >= 2 then print(`||u[` || i || `]||^2 = `, norm_sq); fi;
      if norm_sq <> 0 then
         # linearly independent new orthogonal vector 
         # add to list of u and norm_u_sq
         u := [op(u), w]; norm_u_sq := [op(norm_u_sq), norm_sq]; k := k+1;
         if nargs > 1 then R[k, i] := 1; fi;
      fi;

  od; # for i

  # convert list of vectors to matrix
  if k > 0 then
     Q := Matrix(m,k,u[1]);
     for i from 1 to m do
         for j from 1 to k do
             Q[i,j] := u[j][i];
     od; od;
     if nargs > 1
        then assign(args[2], LinearAlgebra[SubMatrix](R, 1..k, 1..n));
     fi;
     RETURN(Q);
  fi   
end: