123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204 |
- module tadjoint;
- %**********************************************************************%
- % %
- % Calculation of the Triangularizing Adjoint for a given matrix A. The %
- % Triangularizing Adjoint F is a lower triangular matrix. %
- % %
- % Author: Walter Tietze, July 1998 %
- % %
- % The algorithm is due to Arne Storjohann, The Triangularizing Adjoint,%
- % Institut f"ur Wissenschaftliches Rechnen, ETH Z"urich, Switzerland, %
- % http://www.inf.ethz.ch/personal/storjoha %
- % Ref: See Arne's paper in the ISSAC 1997 proceedings %
- % %
- %**********************************************************************%
- symbolic smacro procedure mksq!*mat in_mat;
- %
- % Converts entries in matrix to standard quotients.
- %
- begin scalar tmp_mat,out_mat;
- tmp_mat:= cdr in_mat;
- out_mat:= for each u in tmp_mat collect
- for each v in u collect
- if atom v then v
- else mk!*sq v;
- return 'mat . out_mat;
- end;
- symbolic smacro procedure reval!*mat in_mat;
- %
- % Revals the entries in matrix.
- %
- begin scalar tmp_mat,out_mat;
- tmp_mat:= cdr in_mat;
- out_mat:= for each u in tmp_mat collect
- for each v in u collect
- my_reval v;
- return 'mat . out_mat;
- end;
- symbolic procedure triang_adjoint in_mat;
- %
- % Due to the algorithm of Arne Storjohann this function
- % calculates the lower triangular matrix, the so-called
- % triangularizing adjoint.
- %
- begin scalar u;
- if eqcar(u:=aeval(in_mat),'matrix) then u:=cadr u
- else u:=reval in_mat;
- if not matrixp(u) then
- rederr "Error in triang_adjoint: non matrix input."
- else if not squarep(u) then
- rederr "Error in triang_adjoint: input matrix should be square."
- else return matsm adtriang!* u;
- end;
- put('triang_adjoint,'rtypefn,'getrtypecar);
- symbolic procedure adtriang!* in_mat;
- %
- % Strategy : Calculate the triangularizing adjoint by
- % transforming in_mat to a matrix with lowest 2-potential
- % dimension greater or equal to the dimension of in_mat.
- %
- (begin scalar mat_dim,tmp_mat,l,base;
- integer dim;
- on fast_la;
- dim:=cadr matlength in_mat;
- base:=logb(dim,2);
- mat_dim:= if base-floor(base)=0 then fix base
- else (floor(base) + 1);
- mat_dim:=2**mat_dim;
- if mat_dim>dim then
- tmp_mat:=extend(in_mat,mat_dim - dim,mat_dim - dim,0)
- else tmp_mat:=in_mat;
- tmp_mat:=adjoint!*lu(tmp_mat);
- l:= for i:=1:dim collect i;
- tmp_mat:=sub_matrix(tmp_mat,l,l);
- return tmp_mat;
- end) where !*fast_la = !*fast_la;
- symbolic procedure adjoint!*lu(in_mat);
- %
- % This function calculates iteratively the triangularizing
- % adjoint.
- %
- begin scalar a1,a_tmp, a_tmp1, f1, a4!*,f4!*, subdim0, subdim1, l;
- integer determinant, dim, crrnt_dim;
- dim := cadr matlength in_mat;
- if dim < 2 then return 'mat . list({1});
- crrnt_dim := 1;
- f1:=list('mat,{1});
- while crrnt_dim < dim do
- begin
- subdim0:= for i:=1:crrnt_dim collect i;
- subdim1:= for i:=(crrnt_dim+1):(2*crrnt_dim) collect i;
- a1:=sub_matrix(in_mat,subdim0,subdim0);
- a1:= reval!*mat a1;
- determinant:=0;
- for j:=1:(crrnt_dim) do
- determinant:={'plus,my_reval determinant, my_reval{'times,
- getmat(f1,crrnt_dim,j),getmat(a1,j,crrnt_dim)}};
- if my_reval determinant = 0 then
- << a_tmp := sub_matrix(in_mat,append(subdim0,subdim1),subdim0);
- if rank!-eval(list(a_tmp)) < crrnt_dim then
- << f1:=extend(f1,crrnt_dim,crrnt_dim,0);
- crrnt_dim:= 2*crrnt_dim;
- >>
- else
- << if crrnt_dim=1 then
- <<f1:=extend(f1,1,1,0);
- setmat(f1,2,1,{'times, -1, getmat(in_mat,2,1)});
- crrnt_dim:=2*crrnt_dim;>>
- else <<f1:=compose!*mat(in_mat,f1,subdim0,crrnt_dim);
- crrnt_dim:=2*crrnt_dim;>>;
- >>;
- >>
- else
- <<
- a4!*:= matinverse matsm a1;
- a_tmp:=sub_matrix(in_mat,subdim1,subdim0);
- a4!*:= multm(matsm a_tmp, a4!*);
- a4!* := for each u in a4!* collect
- for each v in u collect
- v:= negsq v;
- a_tmp1:='mat . a4!*;
- a_tmp:=sub_matrix(in_mat,subdim0,subdim1);
- a4!*:= multm(a4!*, matsm a_tmp);
- a4!*:= addm(a4!*,matsm sub_matrix(in_mat,subdim1,subdim1));
- a4!*:= 'mat . a4!*;
- f4!* := adjoint!*lu reval!*mat mksq!*mat a4!*;
- l:= for i:=1:crrnt_dim collect i;
- a_tmp := mult_rows(f4!*,l,determinant);
- a_tmp:=reval!*mat a_tmp;
- f1:=extend(f1,crrnt_dim,crrnt_dim,0);
- f1:=copy_into(a_tmp,f1,crrnt_dim+1,crrnt_dim+1);
- a_tmp:='mat . multm(matsm a_tmp,matsm mksq!*mat a_tmp1);
- a_tmp:=mksq!*mat a_tmp;
- f1:=copy_into(a_tmp,f1,crrnt_dim+1,1);
- crrnt_dim:=crrnt_dim*2;
- >>;
- end;
- return f1;
- end;
- symbolic procedure compose!*mat(in_mat,f1,subdim0,crrnt_dim);
- %
- % Due to the algorithm of Arne Storjohann this function
- % determines the rows of the triangularizing adjoint which
- % can be set equal to zero.
- %
- begin scalar tmp_mat,tmp_row,k;
- for i:=(2*crrnt_dim) step (-1) until crrnt_dim do
- begin
- k:= for j:=1:i collect j;
- if rank!-eval {sub_matrix(in_mat,k,subdim0)} < crrnt_dim then
- << f1:=extend(f1,crrnt_dim,crrnt_dim,0);
- for j:=(i+1):(2*crrnt_dim) do
- begin
- k:= append(k,{j});
- tmp_mat:=sub_matrix(in_mat,k,k);
- tmp_row:=adjoint_lst_row!*(tmp_mat,j);
- for l:=1:j do
- setmat(f1,j,l,nth(cadr tmp_row,l));
- end;
- i:=crrnt_dim-1;
- >>
- end;
- return f1;
- end;
- symbolic procedure adjoint_lst_row!*(in_mat,len);
- %
- % Calculates one row for the triangularizing adjoint in
- % last row of submatrix(len,len) of matrix in_mat.
- %
- begin
- scalar tmp_mat, adj_row,det;
- integer sign;
- if len=1 then return in_mat;
- for j:=1:len do
- begin
- sign := (-1)**(len+j);
- tmp_mat := minor(in_mat, j, len);
- if sign = -1 then
- det:= mk!*sq negsq(simpdet({tmp_mat}))
- else
- det:= mk!*sq simpdet({tmp_mat});
- adj_row:=append(adj_row,{det});
- end;
- return 'mat . {adj_row}
- end;
- endmodule;
- end;
|