linSolGen¶
Solves a real general system of linear equations \(Ax = b\). Using optional arguments, any of several related computations can be performed. These extra tasks include computing the LU factorization of A using partial pivoting, computing the inverse matrix \(A^{-1}\), solving \(A^Tx = b\), or computing the solution of \(Ax = b\) given the LU factorization of A.
Synopsis¶
linSolGen (a, b)
Required Arguments¶
- float
a[[]]
(Input) - Array of size n × n containing the matrix.
- float
b[]
(Input) - Array of size n containing the right-hand side.
Return Value¶
The solution x of the linear system \(Ax = b\). If no solution was
computed, then None
is returned.
Optional Arguments¶
transpose
Solve \(A^Tx = b\).
Default: Solve \(Ax = b\)
factor
, pPvt
, pFactor
(Output)
pPvt
(Output)- An array of length n containing the pivot sequence for the factorization.
pFactor
(Output)- An array of size n × n containing the LU factorization of A with column pivoting. The lower-triangular part of this array contains information necessary to construct L, and the upper-triangular part contains U.
These parameters are input if
solve
is specified. They are output otherwise.
inverse
(Output)- An array of size n × n containing the inverse of the matrix A.
condition
(Output)- A scalar containing an estimate of the \(L_1\) norm condition number
of the matrix A. This option cannot be used with the option
solveOnly
. factorOnly
- Compute the LU factorization of A with partial pivoting. If
factorOnly
is used,factor
is required. The argumentb
is then ignored, and the returned value oflinSolGen
isNone
. solveOnly
- Solve \(Ax = b\) given the LU factorization previously computed by
linSolGen
. By default, the solution to \(Ax = b\) is pointed to bylinSolGen
. IfsolveOnly
is used, argumentfactor
is required, and the argumenta
is ignored. inverseOnly
- Compute the inverse of the matrix A. If
inverseOnly
is used,inverse
is required. The argumentb
is then ignored, and the returned value oflinSolGen
isNone
.
Description¶
The function linSolGen
solves a system of linear algebraic equations
with a real coefficient matrix A. It first computes the LU factorization
of A with partial pivoting such that \(L^{-1}A = U\). Let F be the
matrix pFactor
returned by optional argument factor
. The triangular
matrix U is stored in the upper triangle of F. The strict lower triangle
of F contains the information needed to reconstruct \(L^{-1}\) using
The factors \(P_i\) and \(L_i\) are defined by partial pivoting.
\(P_i\) is the identity matrix with rows i and pPvt[i‑1]
interchanged. \(L_i\) is the identity matrix with \(F_{ji}\) , for
\(j = i + 1, …, n\), inserted below the diagonal in column i.
The factorization efficiency is based on a technique of “loop unrolling and
jamming” by Dr. Leonard J. Harding of the University of Michigan, Ann Arbor,
Michigan. The solution of the linear system is then found by solving two
simpler systems, \(y = L^{-1}b\) and \(x = U^{-1}y\). When the
solution to the linear system or the inverse of the matrix is sought, an
estimate of the \(L_1\) 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 function linSolGen
fails if U,
the upper triangular part of the factorization, has a zero diagonal element.
Examples¶
Example 1¶
This example solves a system of three linear equations. This is the simplest use of the function. The equations follow below:
from numpy import *
from pyimsl.math.linSolGen import linSolGen
from pyimsl.math.writeMatrix import writeMatrix
a = array([[1.0, 3.0, 3.0],
[1.0, 3.0, 4.0],
[1.0, 4.0, 3.0]])
b = [1.0, 4.0, -1.0]
# Solve A*x = b for x
x = linSolGen(a, b)
writeMatrix("Solution, x, of Ax = b", x)
Output¶
Solution, x, of Ax = b
1 2 3
-2 -2 3
Example 2¶
This example solves the transpose problem \(A^Tx = b\) and returns the LU factorization of A with partial pivoting. The same data as the initial example is used, except the solution \(x = A^{-T}b\) is returned in an array allocated in the main program. The L matrix is returned in implicit form.
from numpy import *
from pyimsl.math.linSolGen import linSolGen
from pyimsl.math.writeMatrix import writeMatrix
a = array([[1.0, 3.0, 3.0],
[1.0, 3.0, 4.0],
[1.0, 4.0, 3.0]])
b = [1.0, 4.0, -1.0]
factor = {}
# Solve trans(A)*x = b for x
x = linSolGen(a, b, transpose=True, factor=factor)
writeMatrix("Solution, x, of Ax = b", x)
writeMatrix("LU factors of A", factor['pFactor'])
writeMatrix("Pivot sequence", factor['pPvt'])
Output¶
Solution, x, of Ax = b
1 2 3
4 -4 1
LU factors of A
1 2 3
1 1 3 3
2 -1 1 0
3 -1 0 1
Pivot sequence
1 2 3
1 3 3
Reconstruction of \(L^{-1}\) and U from factor
:¶
\(P_i\) is the identity matrix with row i and row pivot
[i-1]
interchanged.
pivot = 1, 3, 3 |
|
\[\begin{split}P_1 =
\begin{bmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{bmatrix}\end{split}\]
|
row 1 and row pivot [0],
or row 1, are interchanged,
which is still the identity
matrix. |
\[\begin{split}P_2 =
\begin{bmatrix}
1 & 0 & 0 \\
0 & 0 & 1 \\
0 & 1 & 0
\end{bmatrix}\end{split}\]
|
row 2 and row pivot[1] , or
row 3, are interchanged. |
\(L_i\) is the identity matrix with \(F_{ji,}\) for \(j = i + 1,
n\), inserted below the diagonal in column i, where F is factor
:
\[\begin{split}\text{factor} =
\begin{bmatrix}
1 & 3 & 3 \\
-1 & 1 & 0 \\
-1 & 0 & 1
\end{bmatrix}\end{split}\]
|
|
\[\begin{split}L_1 =
\begin{bmatrix}
1 & 0 & 0 \\
-1 & 1 & 0 \\
-1 & 0 & 1
\end{bmatrix}\end{split}\]
|
second and third elements
of column 1 of factor
are inserted below the
diagonal in column 1. |
\[\begin{split}L_2 =
\begin{bmatrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{bmatrix}\end{split}\]
|
third element of column 2
of factor is inserted
below the diagonal in
column 2. |
U is the upper triangle of factor
:
Example 3¶
This example computes the inverse of the 3 × 3 matrix A of the initial
example and solves the same linear system. The matrix product \(C =
A^{-1}A\) is computed and printed. The function matMulRect
is used to
compute C. The approximate result \(C = I\) is obtained.
from numpy import *
from pyimsl.math.linSolGen import linSolGen
from pyimsl.math.matMulRect import matMulRect
from pyimsl.math.writeMatrix import writeMatrix
a = array([[1.0, 3.0, 3.0],
[1.0, 3.0, 4.0],
[1.0, 4.0, 3.0]])
b = [1.0, 4.0, -1.0]
c = [0.5, 0.3, 0.4]
inv = []
# Solve A*x = b for x
x = linSolGen(a, b, inverse=inv)
# Print solution
writeMatrix("Solution, x, of Ax = b", x)
writeMatrix("Input A", a)
writeMatrix("Inverse of A", inv)
c = matMulRect("A*B",
aMatrix=inv,
bMatrix=a)
writeMatrix("Product, inv(A)*A", c)
Output¶
Solution, x, of Ax = b
1 2 3
-2 -2 3
Input A
1 2 3
1 1 3 3
2 1 3 4
3 1 4 3
Inverse of A
1 2 3
1 7 -3 -3
2 -1 0 1
3 -1 1 0
Product, inv(A)*A
1 2 3
1 1 0 0
2 0 1 0
3 0 0 1
Example 4¶
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 with the first solution and is used to compute the second solution. The factorization work done in the first step is avoided in computing the second solution.
from numpy import *
from pyimsl.math.linSolGen import linSolGen
from pyimsl.math.writeMatrix import writeMatrix
a = array([[1.0, 3.0, 3.0],
[1.0, 3.0, 4.0],
[1.0, 4.0, 3.0]])
b = [1.0, 4.0, -1.0]
c = [0.5, 0.3, 0.4]
factor = {}
# Solve A*x = b for x
x = linSolGen(a, b, factor=factor)
writeMatrix("Solution, x, of Ax = b", x)
# Solve for A*y = c for y
y = linSolGen(a, c,
solveOnly=True,
factor=factor)
writeMatrix("Solution, y, of Ay = c", y)
Output¶
Solution, x, of Ax = b
1 2 3
-2 -2 3
Solution, y, of Ay = c
1 2 3
1.4 -0.1 -0.2
Warning Errors¶
IMSL_ILL_CONDITIONED |
The input matrix is too ill-conditioned.
An estimate of the reciprocal of its
\(L_1\) condition number is
“rcond ” = #. The solution might
not be accurate. |
Fatal Errors¶
IMSL_SINGULAR_MATRIX |
The input matrix is singular. |