Function: lu_factor (<M>, <field>) Return a list of the form
[<LU>, <perm>, <fld>], or
[<LU>, <perm>, <fld>, <lower-cnd> <upper-cnd>], where
(1) The matrix <LU> contains the factorization of <M> in a packed form. Packed form means three things: First, the rows of <LU> are permuted according to the list <perm>. If, for example, <perm> is the list
[3,2,1], the actual first row of the <LU> factorization is the third row of the matrix <LU>. Second, the lower triangular factor of m is the lower triangular part of <LU> with the diagonal entries replaced by all ones. Third, the upper triangular factor of <M> is the upper triangular part of <LU>.
(2) When the field is either
complexfield, the numbers <lower-cnd> and <upper-cnd> are lower and upper bounds for the infinity norm condition number of <M>. For all fields, the condition number might not be estimated; for such fields,
lu_factor returns a two item list. Both the lower and upper bounds can differ from their true values by arbitrarily large factors. (See also
The argument <M> must be a square matrix.
The optional argument <fld> must be a symbol that determines a ring or field. The pre-defined fields and rings are:
generalring - the ring of Maxima expressions, (b)
floatfield - the field of floating point numbers of the type double, (c)
complexfield - the field of complex floating point numbers of the type double, (d)
crering - the ring of Maxima CRE expressions, (e)
rationalfield - the field of rational numbers, (f)
runningerror - track the all floating point rounding errors, (g)
noncommutingring - the ring of Maxima expressions where multiplication is the non-commutative dot operator.
When the field is
runningerror, the algorithm uses partial pivoting; for all other fields, rows are switched only when needed to avoid a zero pivot.
Floating point addition arithmetic isnt associative, so the meaning of field differs from the mathematical definition.
A member of the field
runningerror is a two member Maxima list of the form
[x,n],where <x> is a floating point number and
n is an integer. The relative difference between the true value of
x is approximately bounded by the machine epsilon times
n. The running error bound drops some terms that of the order the square of the machine epsilon.
There is no user-interface for defining a new field. A user that is familiar with Common Lisp should be able to define a new field. To do this, a user must define functions for the arithmetic operations and functions for converting from the field representation to Maxima and back. Additionally, for ordered fields (where partial pivoting will be used), a user must define functions for the magnitude and for comparing field members. After that all that remains is to define a Common Lisp structure
mring. The file
mring has many examples.
To compute the factorization, the first task is to convert each matrix entry to a member of the indicated field. When conversion isnt possible, the factorization halts with an error message. Members of the field neednt be Maxima expressions. Members of the
complexfield, for example, are Common Lisp complex numbers. Thus after computing the factorization, the matrix entries must be converted to Maxima expressions.
(%i1) w[i,j] := random (1.0) + %i * random (1.0); (%o1) w := random(1.) + %i random(1.) i, j (%i2) showtime : true$ Evaluation took 0.00 seconds (0.00 elapsed) (%i3) M : genmatrix (w, 100, 100)$ Evaluation took 7.40 seconds (8.23 elapsed) (%i4) lu_factor (M, complexfield)$ Evaluation took 28.71 seconds (35.00 elapsed) (%i5) lu_factor (M, generalring)$ Evaluation took 109.24 seconds (152.10 elapsed) (%i6) showtime : false$
(%i7) M : matrix ([1 - z, 3], [3, 8 - z]); [ 1 - z 3 ] (%o7) [ ] [ 3 8 - z ] (%i8) lu_factor (M, generalring); [ 1 - z 3 ] [ ] (%o8) [[ 3 9 ], [1, 2], generalring] [ ----- - z - ----- + 8 ] [ 1 - z 1 - z ] (%i9) get_lu_factors (%); [ 1 0 ] [ 1 - z 3 ] [ 1 0 ] [ ] [ ] (%o9) [[ ], [ 3 ], [ 9 ]] [ 0 1 ] [ ----- 1 ] [ 0 - z - ----- + 8 ] [ 1 - z ] [ 1 - z ] (%i10) % . % . %; [ 1 - z 3 ] (%o10) [ ] [ 3 8 - z ]
There are also some inexact matches for
?? lu_factor to see them.
(%o1) true (%i2)
m : matrix([1,2],[3,4]);
mlu : lu_factor(m);
[p,l,u] : get_lu_fact...