Artifact e9bf73b9d8ebae90395aa2de47da4f66092103595384578e8d7187124b85ad4d:
- Executable file
r37/packages/linalg/tadjoint.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 6964) [annotate] [blame] [check-ins using] [more...]
- Executable file
r38/packages/linalg/tadjoint.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 6964) [annotate] [blame] [check-ins using]
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;