LUSOL Function
Solves a general system of real or complex linear equations Ax = b.
Usage
result = LUSOL(b[, a])
Input Parameters
b—One-dimensional matrix containing the right-hand side.
a—Two-dimensional matrix containing the coefficient matrix. Element A(i, j) contains the jth coefficient of the ith equation.
Returned Value
result—A one-dimensional array containing the solution of the linear system Ax = b.
Input Keywords
Double—If present and nonzero, double precision is used.
Transpose—If present and nonzero, AH x = b is solved.
Pivot—Specifies a named variable into which the pivot sequence for the factorization, computed by the LUFAC procedure, is stored. Keywords Pivot and Factor must be used together. Keywords Pivot and Condition cannot be used together.
Output Keywords
Factor—Specifies a named variable in which the LU factorization of A, computed by the LUFAC procedure, is stored. The strictly lower-triangular part of this array contains information necessary to construct L, and the upper-triangular part contains U. Keywords Pivot and Factor must be used together. Keywords Factor and Condition cannot be used together.
Condition—Specifies a named variable into which an estimate of the L1 condition number is stored. This keyword cannot be used with keywords Pivot and Factor.
Inverse—Specifies a named variable into which the inverse of the matrix A is stored.
Discussion
Function LUSOL solves a system of linear algebraic equations with a real or complex coefficient matrix A. Any of several related computations can be performed by using keywords. These extra tasks include solving AHx = b or computing the solution of Ax = b given the LU factorization of A. The function first 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 j ≥ i. 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 matrix 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. The solution of the linear system is then found by solving two simpler systems, y = L–1b and x = U–1y. When the solution to the linear system or 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. Function LUSOL fails if U, the upper-triangular part of the factorization, has a zero diagonal element.
Example 1: Solving a System
This example solves a system of three linear equations. This is the simplest use of the function. The equations are as follows:
x0 + 3x1 + 3x2 = 1
x0 + 3x1 + 4x2 = 4
x0 + 4x1 + 3x2 = –1
; Input a matrix containing the coefficients.
RM, a, 3, 3
row 0: 1 3 3
row 1: 1 3 4
row 2: 1 4 3
; Input a vector containing the right-hand side.
RM, b, 3, 1
row 0: 1
row 1: 4
row 2: -1
; Call LUSOL to compute the solution.
x = LUSOL(b, a)
; Print solution and residual.
PM, x, Title = 'Solution'
; PV-WAVE prints the following:
; Solution
; -2.00000
; -2.00000
; 3.00000
PM, a # x - b, Title = 'Residual'
; PV-WAVE prints the following:
; Residual
; 0.00000
; 0.00000
; 0.00000
Example 2: Transpose Problem
In this example, the transpose problem AHx = b is solved.
; Input the matrix containing the coefficients.
RM, a, 3, 3
row 0: 1 3 3
row 1: 1 3 4
row 2: 1 4 3
; Input the vector containing the right-hand side.
RM, b, 3, 1
row 0: 1
row 1: 4
row 2: -1
; Call LUSOL with keyword Transpose set.
x = LUSOL(b, a, /Transpose)
; Print the solution and the residual.
PM, x, Title = 'Solution'
; PV-WAVE prints the following:
; Solution
; 4.00000
; -4.00000
; 1.00000
PM, TRANSPOSE(a) # x - b, Title = 'Residual'
; PV-WAVE prints the following:
; Residual
; 0.00000
; 0.00000
; 0.00000
Example 3: Solving with Multiple Right-hand Sides
This example computes the solution of two systems. Only the right-hand sides differ. The matrix and first right-hand side are given in the initial example. The second right-hand side is the vector c = [0.5, 0.3, 0.4]T. The factorization information is computed by procedure LUFAC and is used to compute the solutions in calls to LUSOL.
; Input the coefficient matrix.
RM, a, 3, 3
row 0: 1 3 3
row 1: 1 3 4
row 2: 1 4 3
; Input the first right-hand side.
RM, b, 3, 1
row 0: 1
row 1: 4
row 2: -1
; Input the second right-hand side.
RM, c, 3, 1
row 0: .5
row 1: .3
row 2: .4
; Call LUFAC to factor the coefficient matrix.
LUFAC, a, pvt, fac
; Call LUSOL with factored form of the coefficient matrix and the
; first right-hand side.
x = LUSOL(b, Factor = fac, Pivot = pvt)
; Print the solution of Ax = b.
PM, x, Title = 'Solution'
; PV-WAVE prints the following:
; Solution
; -2.00000
; -2.00000
; 3.00000
PM, a # x - b, Title = 'Residual'
; PV-WAVE prints the following:
; Residual
; 0.00000
; 0.00000
; 0.00000
; Call LUSOL with factored form of the coefficient matrix and the
; second right-hand side.
y = LUSOL(c, Factor = fac, Pivot = pvt)
; Print the solution of Ax = b.
PM, y, Title = 'Solution'
; PV-WAVE prints the following:
; Solution
; 1.40000
; -0.100000
; -0.200000
PM, a # y - c, $
Title = 'Residual', Format = '(f8.5)'
; PV-WAVE prints the following:
; Residual
; 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.
