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 argument b is then ignored, and the returned value of linSolGen is None.

solveOnly

Solve $Ax = b$ given the LU factorization previously computed by linSolGen. By default, the solution to $Ax = b$ is pointed to by linSolGen. If solveOnly is used, argument factor is required, and the argument a is ignored.

inverseOnly

Compute the inverse of the matrix A. If inverseOnly is used, inverse is required. The argument b is then ignored, and the returned value of linSolGen is None.

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

$$L^{-1} = L_{n-1}P_{n-1} \ldots L_1P_1$$

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:

$$x_1 + 3x_2 + 3x_3 = 1$$

$$x_1 + 3x_2 + 4x_3 = 4$$

$$x_1 + 4x_2 + 3x_3 = −1$$

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:

$$L^{-1} = L_2 P_2 L_1 P_1$$

$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.

$$\begin{split}L^{-1} = L_2P_2L_1P_1 = \begin{bmatrix} 1 & 0 & 0 \\ -1 & 0 & 1 \\ -1 & 1 & 0 \end{bmatrix}\end{split}$$

U is the upper triangle of factor:

$$\begin{split}U = \begin{bmatrix} 1 & 3 & 3 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}\end{split}$$

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.