share/matrix/invert.mac

   1
   2 /* the following routines compute inverses and adjoints of matrices */
   3 /* if ratmx is false [the default] then the elements will not be converted
   4 to cre-form */
   5
   6 adjoint(mat):= block([adj,n], n:length(mat), adj:ident(n),
   7                      if n#1 then
   8                        for i thru n do for j thru n do
   9                            adj[i,j]:(-1)^(i+j)*determinant(minor(mat,j,i)),
  10                      adj)$
  11
  12 invert(mat):= block([adj,ans], adj:adjoint(mat),
  13                         ans:block([scalarmatrixp:true],
  14                                 adj/(row(mat,1).col(adj,1))),
  15                 if scalarmatrixp=true and length(mat)=1
  16 /* row(mat,1).col(adj,1) = determinant(mat) */
  17                          then ans[1,1] else ans)$
  18