Hi Nuno,
it is nice to read you in this list :)
On Tue, 9 Feb 2010, Nuno Santos wrote:
> A: matrix( [1, 2, 3], [4, 5, 6], [7, 8, 9] );
> lu_factor(A)$
> get_lu_factors(%);
>
> and that would give me L and U in %[2], and [%3], respectively.
>
> However, the LU decomposition, but I am looking for it's:
> L - lower triangular matrix (not lower triangular matrix with unitary
> diagonal, which is what get_lu_factors returns)
> U - as upper triangular matrix with unitary diagonal. (not upper triangular
> without unitary diagonal, which is what get_lu_factors returns)
What I usually did was to construct a diagonal matrix D with its
diagonal equal to that of U and then rewrite LU = (L.D).(inverse(D).U)
L2=L.D will no longer have ones in the diagonal and U2=inverse(D).U will
have ones in the diagonal. That is done in Maxima in this way:
[P,L,U]: get_lu_factors(lu_factor(A))$
D:genmatrix(lambda ([i,j], if (i=j) then U[i,j] else 0), 3, 3)$
L2: L.D;
U2: invert(D).U;
However, in your particular example that method fails because U has a zero
in the diagonal, so inverse(D) does not exist.
A more clever way of doing it is what Leo has proposed:
On Tue, 2010-02-09 at 15:18 +0000, Leo Butler wrote:
>
> I think if you want to change this behaviour, you should compute the
> factorization of transpose(A), then transpose L & U.
(%i69) A: matrix( [1, 2, 3], [4, 5, 6], [7, 8, 9] )$
(%i70) [P,L,U]: get_lu_factors(lu_factor(transpose(A)))$
(%i71) L2: transpose(U);
(%o71) matrix([1,0,0],[4,-3,0],[7,-6,0])
(%i72) U2: transpose(L);
(%o72) matrix([1,2,3],[0,1,2],[0,0,1])
(%i73) L2.U2;
(%o73) matrix([1,2,3],[4,5,6],[7,8,9])
Regards,
Jaime Villate