linSolPosdef¶

Solves a real symmetric positive definite system of linear equations \(Ax = b\). Using optional arguments, any of several related computations can be performed. These extra tasks include computing the Cholesky factor, L, of A such that \(A = LL^T\), computing the inverse matrix \(A^{-1}\), or computing the solution of \(Ax = b\) given the Cholesky factor, L.

Synopsis¶

linSolPosdef ( 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 symmetric positive definite linear system \(Ax = b\). If no solution was computed, then None is returned.

Optional Arguments¶

factor (Output): An array of size n × n containing the \(LL^T\) factorization of A. The lower‑triangular part of this array contains L and the upper-triangular part contains \(L^T\).
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. Do not use this option with solveOnly.
factorOnly: Compute the Cholesky factorization \(LL^T\) of A. If factorOnly is used, factor is required. The argument b is then ignored, and the returned value of linSolPosdef is None.
solveOnly: Solve \(Ax = b\) given the \(LL^T\) factorization previously computed by linSolPosdef. By default, the solution to \(Ax = b\) is pointed to by linSolPosdef. 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 linSolPosdef is None.

Description¶

The function linSolPosdef solves a system of linear algebraic equations having a symmetric positive definite coefficient matrix A. The function first computes the Cholesky factorization \(LL^T\) of A. The solution of the linear system is then found by solving the two simpler systems, \(y = L^{-1}b\) and \(x = L^{-T}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 linSolPosdef fails if L, the lower-triangular matrix in the factorization, has a zero diagonal element.

Examples¶

Example 1¶

A system of three linear equations with a symmetric positive definite coefficient matrix is solved in this example. The equations are listed below:

\[x_1 − 3x_2 + 2x_3 = 27\]

\[−3x_1 + 10x_2 − 5x_3 = −78\]

\[2x_1 − 5x_2 + 6x_3 = 64\]

from numpy import *
from pyimsl.math.linSolPosdef import linSolPosdef
from pyimsl.math.writeMatrix import writeMatrix

factor = []
a = array([[1.0, -3.0, 2.0],
           [-3.0, 10.0, -5.0],
           [2.0, -5.0, 6.0]])
b = array([27.0, -78.0, 64.0])

# Solve Ax=b for x
x = linSolPosdef(a, b)
writeMatrix("Solution, x, of Ax = b", x)

Output¶

 
       Solution, x, of Ax = b
          1            2            3
          1           -4            7

Example 2¶

This example solves the same system of three linear equations as in the initial example, but this time returns the \(LL^T\) factorization of A.

from numpy import *
from pyimsl.math.linSolPosdef import linSolPosdef
from pyimsl.math.writeMatrix import writeMatrix

factor = []
a = array([[1.0, -3.0, 2.0],
           [-3.0, 10.0, -5.0],
           [2.0, -5.0, 6.0]])
b = array([27.0, -78.0, 64.0])

# Solve Ax=b for x
x = linSolPosdef(a, b, factor=factor)
writeMatrix("Solution, x, of Ax = b", x)
writeMatrix("Cholesky factor L, and trans(L), of A", factor)

Output¶

 
       Solution, x, of Ax = b
          1            2            3
          1           -4            7
 
  Cholesky factor L, and trans(L), of A
             1            2            3
1            1           -3            2
2           -3            1            1
3            2            1            1

Example 3¶

This example solves the same system as in the initial example, but given the Cholesky factors of A.

from numpy import *
from pyimsl.math.linSolPosdef import linSolPosdef
from pyimsl.math.writeMatrix import writeMatrix

factor = [[1.0, -3.0, 2.0],
          [-3.0, 1.0, 1.0],
          [2.0, 1.0, 1.0]]
a = array([[1.0, -3.0, 2.0],
           [-3.0, 10.0, -5.0],
           [2.0, -5.0, 6.0]])
b = array([27.0, -78.0, 64.0])

# Solve Ax=b for x
x = linSolPosdef(a, b, factor=factor, solveOnly=True)
writeMatrix("Solution, x, of Ax = b", x)

Output¶

 
       Solution, x, of Ax = b
          1            2            3
          1           -4            7

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_NONPOSITIVE_MATRIX`	The leading # by # submatrix of the input matrix is not positive definite.
`IMSL_SINGULAR_MATRIX`	The input matrix is singular.
`IMSL_SINGULAR_TRI_MATRIX`	The input triangular matrix is singular. The index of the first zero diagonal element is #.