reduce-historical/r37/packages/linalg/gramschm.red

module gramchmd;

%**********************************************************************%
%                                                                      %
% Computation of the Gram Schmidt Orthonormalisation process. The      %
% input vectors are represented by lists.                              %
%                                                                      %
% Authors: Karin Gatermann (used symbolically in her symmetry package).%
%          Matt Rebbeck (first few lines of code that make it          %
%                        available from the user level).   May 1994.   %
%                                                                      %
%**********************************************************************%



put('gram_schmidt,'psopfn,'gram_schmidt1); % To allow variable input.
symbolic procedure gram_schmidt1(vec_list);
  %
  % Can take a list of lists(which are representing vectors) or any 
  % number of arguments each being a list(again which represent the 
  % vectors).
  %
  % Karin used lists of standard quotient elements as vectors so a bit 
  % of fiddling is required to get the input/output right.
  %
  begin
    scalar gs_list;
    % Deal with the possibility of the user entering a list of lists.
    if pairp vec_list and pairp car vec_list and caar vec_list = 'list 
     and pairp cdar vec_list and pairp cadar vec_list 
      and caadar vec_list = 'list 
       then vec_list := cdar vec_list;
    vec_list := convert_to_sq(vec_list);
    % This bit does all the real work.
    gs_list := gram!+schmid(vec_list);
    return convert_from_sq(gs_list);
  end;



symbolic procedure convert_to_sq(vec_list);
  %
  % Takes algebraic list and converts to sq form for input into 
  % GramSchmidt.
  %
  begin
    scalar sq_list;
    sq_list := for each list in vec_list collect
      for each elt in cdr list collect simp!* elt;
    return sq_list;
  end;



symbolic procedure convert_from_sq(sq_list);
  %
  % Converts sq_list to a readable (from algebraic mode) form.
  %
  begin
    scalar gs_list;
    gs_list := 'list.for each elt1 in sq_list collect 
                'list.for each elt in elt1 collect prepsq elt;
    return gs_list;
  end;


%
% All the rest is Karin's.
%

symbolic procedure vector!+p(vector1);
  %
  % returns the length of a vector
  % vector -- list of sqs
  %
  begin
    if length(vector1)>0 then return t;
  end;



symbolic procedure get!+vec!+dim(vector1);
  %
  % returns the length of a vector
  % vector -- list of sqs
  %
  begin
    return length(vector1);
  end;



symbolic procedure get!+vec!+entry(vector1,elem);
  %
  % returns the length of a vector
  % vector -- list of sqs
  %
  begin
    return nth(vector1,elem);
  end;



symbolic procedure mk!+vec!+add!+vec(vector1,vector2);
  %
  % returns a vector= vector1+vector2 (internal structure)
  %
  begin
    scalar ent,res,h;
    res:=for ent:=1:get!+vec!+dim(vector1) collect
    <<
      h:= addsq(get!+vec!+entry(vector1,ent),
             get!+vec!+entry(vector2,ent));
      h:=subs2 h where !*sub2=t;
      h
    >>;
    return res;
  end;



symbolic procedure mk!+squared!+norm(vector1);
  %
  % returns a scalar= sum vector_i^2 (internal structure)
  %
  begin
    return mk!+inner!+product(vector1,vector1);
  end;



symbolic procedure my!+nullsq!+p(scal);
  %
  % returns true, if ths sq is zero
  %
  begin
    if null(numr( scal)) then return t;
  end;



symbolic procedure mk!+null!+vec(dimen);
  %
  % returns a vector of zeros
  %
  begin
    scalar nullsq,i,res;
    nullsq:=(nil ./ 1);
    res:=for i:=1:dimen collect nullsq;
    return res;
  end;



symbolic procedure null!+vec!+p(vector1);
  %
  % returns a true, if vector is the zero vector
  begin
    if my!+nullsq!+p(mk!+squared!+norm(vector1)) then return t;
  end;



symbolic procedure mk!+normalize!+vector(vector1);
  %
  % returns a normalized vector (internal structure)
  %
  begin
    scalar scalo,vecres;
    scalo:=simp!* {'sqrt, mk!*sq(mk!+squared!+norm(vector1))};
    if my!+nullsq!+p(scalo) then 
     vecres:= mk!+null!+vec(get!+vec!+dim(vector1)) 
     else
     <<
       scalo:=simp prepsq scalo;
       scalo:=quotsq((1 ./ 1),scalo);
       vecres:= mk!+scal!+mult!+vec(scalo,vector1);
     >>;
    return vecres;
  end;



symbolic procedure mk!+Gram!+Schmid(vectorlist,vector1);
  %
  % returns a vectorlist of orthonormal vectors
  % assumptions: vectorlist is orthonormal basis, internal structure
  %
  begin
    scalar i,orthovec,scalo,vectors;
    orthovec:=vector1;
    for i:=1:(length(vectorlist)) do
    <<
      scalo:= negsq(mk!+inner!+product(orthovec,nth(vectorlist,i)));
      orthovec:=mk!+vec!+add!+vec(orthovec,
         mk!+scal!+mult!+vec(scalo,nth(vectorlist,i)));
    >>;
    orthovec:=mk!+normalize!+vector(orthovec);
    if null!+vec!+p(orthovec) then vectors:=vectorlist 
     else vectors:=add!+vector!+to!+list(orthovec,vectorlist);
    return vectors;
  end;



symbolic procedure Gram!+Schmid(vectorlist);
  %
  % returns a vectorlist of orthonormal vectors
  %
  begin
    scalar ortholist,i;
    if length(vectorlist)<1
      then rederr "Error in Gram Schmidt: no input.";
    if vector!+p(car vectorlist) then ortholist:=nil
     else rederr "Error in Gram_schmidt: empty input.";
    for i:=1:length(vectorlist) do
     ortholist:=mk!+Gram!+Schmid(ortholist,nth(vectorlist,i));
    return ortholist;
  end;



symbolic procedure add!+vector!+to!+list(vector1,vectorlist);
  %
  % returns a list of vectors consisting of vectorlist 
  % and the vector1 at the end
  % internal structure
  begin
    return append(vectorlist,list(vector1));
  end;



symbolic procedure mk!+inner!+product(vector1,vector2);
  %
  % returns the inner product of vector1 and vector2 
  % (internal structure)
  %
  begin
    scalar z,sum,vec2;
    if not(get!+vec!+dim(vector1) = get!+vec!+dim(vector2)) then rederr
   "Error in Gram_schmidt: each list in input must be the same length.";
    sum:=(nil ./ 1);
    if !*complex then vec2:=mk!+conjugate!+vec(vector2) 
     else vec2:=vector2;
    for z:=1:get!+vec!+dim(vector1) do 
     sum:=addsq(sum,multsq(
            get!+vec!+entry(vector1,z),
            get!+vec!+entry(vec2,z)
                           )
                );
    sum:=subs2 sum where !*sub2=t;
    return sum;
  end;



symbolic procedure mk!+scal!+mult!+vec(scal1,vector1);
  %
  % returns a vector= scalar*vector (internal structure)
  %
  begin
    scalar entry,res,h;
    res:=for each entry in vector1 collect
    <<
      h:=multsq(scal1,entry);
      h:=subs2 h where !*sub2=t;
      h
    >>;
    return res;
  end;


endmodule;  % gram_schmidt.

end;