LUFAC Procedure
Computes the LU factorization of a real or complex matrix.
Usage
LUFAC, a[, pivot[, fac]]
Input Parameters
a—Two-dimensional matrix containing the coefficient matrix. Element A (i, j) contains the jth coefficient of the ith equation.
Output Parameters
pivot—One-dimensional matrix containing the pivot sequence of the factorization.
fac—Two-dimensional matrix containing the LU factorization of A. The strictly lower-triangular part of this matrix contains information necessary to construct L, and the upper-triangular part contains U.
Input Keywords
Double—If present and nonzero, double precision is used.
Transpose—If present and nonzero, ATX=b is solved.
Output Keywords
Condition—Specifies a named variable into which an estimate of the L1 condition number is stored.
L—Specifies a named variable into which the strictly lower-triangular matrix L of the LU factorization is stored.
U—Specifies a named variable into which the upper-triangular matrix U of the LU factorization is stored.
PA—Specifies a named variable into which the matrix resulting from applying the pivot permutation to A is stored.
Inverse—Specifies a named variable into which the inverse of the matrix A is stored.
Discussion
Any of several related computations can be performed by using keywords. These extra tasks include computing the LU factorization of AT, computing an estimate of the L1 condition number, and returning L or U separately.
The LUFAC procedure computes the LU factorization of A with partial pivoting such that L–1PA = U. The matrix U is upper-triangular, while L–1A ≡ Pn – 1 Ln – 2Pn – 2 ... L0 P0 A ≡ U. The factors Pi and Li are defined by the partial pivoting. Each Pi is an interchange of row i with row i ≥ j. Thus, Pi is defined by that value of j. Every Li = mieiT is an elementary elimination matrix. The vector mi is zero in entries 0, ..., i – 1. This vector is stored as column i in the strictly lower-triangular part of the working array containing the decomposition information.
The factorization efficiency is based on a technique of “loop unrolling and jamming” due to Dr. Leonard J. Harding of the University of Michigan, Ann Arbor, Mich. When the inverse of the matrix is sought, an estimate of the L1 condition number of A is computed using the same algorithm as in Dongarra et al. (1979). If the estimated condition number is greater than 1/ε (where ε is the machine precision), a warning message is issued. This indicates that very small changes in A may produce large changes in the solution x. The LUFAC procedure fails if U, the upper-triangular part of the factorization, has a zero diagonal element.
Example 1
This example computes the LU factorization of a matrix and prints it out in the default form with the information needed to construct L and U combined in one array. The matrix is as follows:
; Input the matrix to be factored.
RM, a, 3, 3
row 0: 1 3 3
row 1: 1 3 4
row 2: 1 4 3
; Factor the matrix by calling LUFAC.
LUFAC, a, pvt, fac
; Print the results.
PM, fac, Title = 'LU factors of A'
; PV-WAVE prints the following:
; LU factors of A
; 1.00000 3.00000 3.00000
; -1.00000 1.00000 0.00000
; -1.00000 -0.00000 1.00000
PM, pvt, Title = 'Pivot sequence'
; PV-WAVE prints the following:
; Pivot sequence
; 1
; 3
; 3
Example 2
This example computes the factorization, uses keywords to return the factorization in separate named variables, and returns the original matrix after the pivot permutation is applied.
; Input the matrix to be factored.
RM, a, 3, 3
row 0: 1 3 3
row 1: 1 3 4
row 2: 1 4 3
; Call LUFAC with the keywords L and U.
LUFAC, a, L = l, U = u, PA = pa
; Print the results.
PM, l, Title = 'L'
; PV-WAVE prints the following:
; L
; 1.00000 0.00000 0.00000
; 1.00000 1.00000 0.00000
; 1.00000 0.00000 1.00000
PM, u, Title = 'U'
; PV-WAVE prints the following:
; U
; 1.00000 3.00000 3.00000
; 0.00000 1.00000 0.00000
; 0.00000 0.00000 1.00000
PM, l # u - pa, Title = 'Residual: L # U - PA'
; PV-WAVE prints the following:
; Residual: L # U - PA
; 0.00000 0.00000 0.00000
; 0.00000 0.00000 0.00000
; 0.00000 0.00000 0.00000
Warning Errors
MATH_ILL_CONDITIONED—Input matrix is too ill-conditioned. An estimate of the reciprocal of its L1 condition number is #. The solution might not be accurate.
Fatal Errors
MATH_SINGULAR_MATRIX—Input matrix is singular.
