Stefano Ferri wrote:
> Well, I would certainly do that, but since I don't know the details of
> the implementation I couldn't give more than an example...
>
>
There is a small text announcement.txt in share/linearalgebra which explains
what all this is about. It says:
5) LU factorization
There is code for LU factorization. It uses partial pivoting when
the matrix is numeric.
Recall M=LU where L is lower triangular with 1 on the diagonal,
and U is upper triangular. Then det(M) is easy to compute as the
product of the diagonal U elements, and inversion is also easier:
To solve MX=B one solves first LY=B, then UX=Y, both being triangular
hence easy. Of course vanishing diagonal U elements cause problems, hence
the "pivot" property.
In the above announcement, such inversion is described with the example:
niobe% maxima
(%i1) m : hilbert_matrix(5)$
(%i2) m : lu_factor(m)$
(%i3) lu_backsub(m, matrix([1],[1],[1],[1],[1]));
[ 5 ]
[ ]
[ - 120 ]
[ ]
(%o3) [ 630 ]
[ ]
[ - 1120 ]
[ ]
[ 630 ]
(%i4) hilbert_matrix(5) . %;
[ 1 ]
[ ]
[ 1 ]
[ ]
(%o4) [ 1 ]
[ ]
[ 1 ]
[ ]
[ 1 ]
As you see applying m on %o3 reproduces b that is matrix([1],[1],[1],[1],
[1]), hence %o3 is the solution of mx=b.
Finally linsolve_by_lu defined in lu.lisp, just checks its arguments
and calls lu_factor and lu_backsub as above. You can also specify a field
for computations.
(%i11) m : hilbert_matrix(5)$
(%i12) linsolve_by_lu(m,matrix([1],[1],[1],[1],[1]));
[ 5 ]
[ ]
[ - 120 ]
[ ]
(%o12) [[ 630 ], false]
[ ]
[ - 1120 ]
[ ]
[ 630 ]
Hope that helps documenting the function in question. In general for big
inversion or determinant problems, doing firts the LU factorization seems to
help keeping the problem manageable.
--
Michel Talon