tadjoint.red 6.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204
  1. module tadjoint;
  2. %**********************************************************************%
  3. % %
  4. % Calculation of the Triangularizing Adjoint for a given matrix A. The %
  5. % Triangularizing Adjoint F is a lower triangular matrix. %
  6. % %
  7. % Author: Walter Tietze, July 1998 %
  8. % %
  9. % The algorithm is due to Arne Storjohann, The Triangularizing Adjoint,%
  10. % Institut f"ur Wissenschaftliches Rechnen, ETH Z"urich, Switzerland, %
  11. % http://www.inf.ethz.ch/personal/storjoha %
  12. % Ref: See Arne's paper in the ISSAC 1997 proceedings %
  13. % %
  14. %**********************************************************************%
  15. symbolic smacro procedure mksq!*mat in_mat;
  16. %
  17. % Converts entries in matrix to standard quotients.
  18. %
  19. begin scalar tmp_mat,out_mat;
  20. tmp_mat:= cdr in_mat;
  21. out_mat:= for each u in tmp_mat collect
  22. for each v in u collect
  23. if atom v then v
  24. else mk!*sq v;
  25. return 'mat . out_mat;
  26. end;
  27. symbolic smacro procedure reval!*mat in_mat;
  28. %
  29. % Revals the entries in matrix.
  30. %
  31. begin scalar tmp_mat,out_mat;
  32. tmp_mat:= cdr in_mat;
  33. out_mat:= for each u in tmp_mat collect
  34. for each v in u collect
  35. my_reval v;
  36. return 'mat . out_mat;
  37. end;
  38. symbolic procedure triang_adjoint in_mat;
  39. %
  40. % Due to the algorithm of Arne Storjohann this function
  41. % calculates the lower triangular matrix, the so-called
  42. % triangularizing adjoint.
  43. %
  44. begin scalar u;
  45. if eqcar(u:=aeval(in_mat),'matrix) then u:=cadr u
  46. else u:=reval in_mat;
  47. if not matrixp(u) then
  48. rederr "Error in triang_adjoint: non matrix input."
  49. else if not squarep(u) then
  50. rederr "Error in triang_adjoint: input matrix should be square."
  51. else return matsm adtriang!* u;
  52. end;
  53. put('triang_adjoint,'rtypefn,'getrtypecar);
  54. symbolic procedure adtriang!* in_mat;
  55. %
  56. % Strategy : Calculate the triangularizing adjoint by
  57. % transforming in_mat to a matrix with lowest 2-potential
  58. % dimension greater or equal to the dimension of in_mat.
  59. %
  60. (begin scalar mat_dim,tmp_mat,l,base;
  61. integer dim;
  62. on fast_la;
  63. dim:=cadr matlength in_mat;
  64. base:=logb(dim,2);
  65. mat_dim:= if base-floor(base)=0 then fix base
  66. else (floor(base) + 1);
  67. mat_dim:=2**mat_dim;
  68. if mat_dim>dim then
  69. tmp_mat:=extend(in_mat,mat_dim - dim,mat_dim - dim,0)
  70. else tmp_mat:=in_mat;
  71. tmp_mat:=adjoint!*lu(tmp_mat);
  72. l:= for i:=1:dim collect i;
  73. tmp_mat:=sub_matrix(tmp_mat,l,l);
  74. return tmp_mat;
  75. end) where !*fast_la = !*fast_la;
  76. symbolic procedure adjoint!*lu(in_mat);
  77. %
  78. % This function calculates iteratively the triangularizing
  79. % adjoint.
  80. %
  81. begin scalar a1,a_tmp, a_tmp1, f1, a4!*,f4!*, subdim0, subdim1, l;
  82. integer determinant, dim, crrnt_dim;
  83. dim := cadr matlength in_mat;
  84. if dim < 2 then return 'mat . list({1});
  85. crrnt_dim := 1;
  86. f1:=list('mat,{1});
  87. while crrnt_dim < dim do
  88. begin
  89. subdim0:= for i:=1:crrnt_dim collect i;
  90. subdim1:= for i:=(crrnt_dim+1):(2*crrnt_dim) collect i;
  91. a1:=sub_matrix(in_mat,subdim0,subdim0);
  92. a1:= reval!*mat a1;
  93. determinant:=0;
  94. for j:=1:(crrnt_dim) do
  95. determinant:={'plus,my_reval determinant, my_reval{'times,
  96. getmat(f1,crrnt_dim,j),getmat(a1,j,crrnt_dim)}};
  97. if my_reval determinant = 0 then
  98. << a_tmp := sub_matrix(in_mat,append(subdim0,subdim1),subdim0);
  99. if rank!-eval(list(a_tmp)) < crrnt_dim then
  100. << f1:=extend(f1,crrnt_dim,crrnt_dim,0);
  101. crrnt_dim:= 2*crrnt_dim;
  102. >>
  103. else
  104. << if crrnt_dim=1 then
  105. <<f1:=extend(f1,1,1,0);
  106. setmat(f1,2,1,{'times, -1, getmat(in_mat,2,1)});
  107. crrnt_dim:=2*crrnt_dim;>>
  108. else <<f1:=compose!*mat(in_mat,f1,subdim0,crrnt_dim);
  109. crrnt_dim:=2*crrnt_dim;>>;
  110. >>;
  111. >>
  112. else
  113. <<
  114. a4!*:= matinverse matsm a1;
  115. a_tmp:=sub_matrix(in_mat,subdim1,subdim0);
  116. a4!*:= multm(matsm a_tmp, a4!*);
  117. a4!* := for each u in a4!* collect
  118. for each v in u collect
  119. v:= negsq v;
  120. a_tmp1:='mat . a4!*;
  121. a_tmp:=sub_matrix(in_mat,subdim0,subdim1);
  122. a4!*:= multm(a4!*, matsm a_tmp);
  123. a4!*:= addm(a4!*,matsm sub_matrix(in_mat,subdim1,subdim1));
  124. a4!*:= 'mat . a4!*;
  125. f4!* := adjoint!*lu reval!*mat mksq!*mat a4!*;
  126. l:= for i:=1:crrnt_dim collect i;
  127. a_tmp := mult_rows(f4!*,l,determinant);
  128. a_tmp:=reval!*mat a_tmp;
  129. f1:=extend(f1,crrnt_dim,crrnt_dim,0);
  130. f1:=copy_into(a_tmp,f1,crrnt_dim+1,crrnt_dim+1);
  131. a_tmp:='mat . multm(matsm a_tmp,matsm mksq!*mat a_tmp1);
  132. a_tmp:=mksq!*mat a_tmp;
  133. f1:=copy_into(a_tmp,f1,crrnt_dim+1,1);
  134. crrnt_dim:=crrnt_dim*2;
  135. >>;
  136. end;
  137. return f1;
  138. end;
  139. symbolic procedure compose!*mat(in_mat,f1,subdim0,crrnt_dim);
  140. %
  141. % Due to the algorithm of Arne Storjohann this function
  142. % determines the rows of the triangularizing adjoint which
  143. % can be set equal to zero.
  144. %
  145. begin scalar tmp_mat,tmp_row,k;
  146. for i:=(2*crrnt_dim) step (-1) until crrnt_dim do
  147. begin
  148. k:= for j:=1:i collect j;
  149. if rank!-eval {sub_matrix(in_mat,k,subdim0)} < crrnt_dim then
  150. << f1:=extend(f1,crrnt_dim,crrnt_dim,0);
  151. for j:=(i+1):(2*crrnt_dim) do
  152. begin
  153. k:= append(k,{j});
  154. tmp_mat:=sub_matrix(in_mat,k,k);
  155. tmp_row:=adjoint_lst_row!*(tmp_mat,j);
  156. for l:=1:j do
  157. setmat(f1,j,l,nth(cadr tmp_row,l));
  158. end;
  159. i:=crrnt_dim-1;
  160. >>
  161. end;
  162. return f1;
  163. end;
  164. symbolic procedure adjoint_lst_row!*(in_mat,len);
  165. %
  166. % Calculates one row for the triangularizing adjoint in
  167. % last row of submatrix(len,len) of matrix in_mat.
  168. %
  169. begin
  170. scalar tmp_mat, adj_row,det;
  171. integer sign;
  172. if len=1 then return in_mat;
  173. for j:=1:len do
  174. begin
  175. sign := (-1)**(len+j);
  176. tmp_mat := minor(in_mat, j, len);
  177. if sign = -1 then
  178. det:= mk!*sq negsq(simpdet({tmp_mat}))
  179. else
  180. det:= mk!*sq simpdet({tmp_mat});
  181. adj_row:=append(adj_row,{det});
  182. end;
  183. return 'mat . {adj_row}
  184. end;
  185. endmodule;
  186. end;