linSolPosdefBand¶

Solves a real symmetric positive definite system of linear equations Ax = b in band symmetric storage mode. Using optional arguments, any of several related computations can be performed. These extra tasks include computing the \(R^TR\) Cholesky factorization of A, computing the solution of \(Ax = b\) given the Cholesky factorization of A, or estimating the \(L_1\) condition number of A.

Synopsis¶

linSolPosdefBand (a, ncoda, b)

The type double function is linSolPosdefBand.

Required Arguments¶

float a[[]] (Input): Array of size (ncoda + 1) × n containing the n × n positive definite band coefficient matrix in band symmetric storage mode.
int ncoda (Input): Number of upper codiagonals of 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¶

factor, (Output)

An array of size (ncoda + 1) × n containing the \(LL^T\) factorization of A.

These parameters are “Input” if solve is specified. They are “Output” otherwise.

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 \(LL^T\) factorization of A. If factorOnly is used, either factor is required. The argument b is then ignored, and the returned value of linSolPosdefBand is None.

solveOnly

Solve \(Ax = b\) given the \(LL^T\) factorization previously computed by linSolPosdefBand. By default, the solution to \(Ax = b\) is pointed to by linSolPosdefBand. If solveOnly is used, argument factor is required and the argument a is ignored.

Description¶

The function linSolPosdefBand solves a system of linear algebraic equations with a real symmetric positive definite band coefficient matrix A. It computes the \(R^TR\) Cholesky factorization of A. R is an upper triangular band matrix.

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 Higham’s modifications to Hager’s method, as given in Higham (1988). 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 linSolPosdefBand fails if any submatrix of R is not positive definite or if R has a zero diagonal element. These errors occur only if A is very close to a singular matrix or to a matrix which is not positive definite.

The function linSolPosdefBand is partially based on the LINPACK subroutines CPBFA and SPBSL; see Dongarra et al. (1979).

Example 1¶

Solves a system of linear equations \(Ax = b\), where

\[\begin{split}A = \begin{bmatrix} 2 & 0 & -1 & 0 \\ 0 & 4 & 2 & 1 \\ -1 & 2 & 7 & -1 \\ 0 & 1 & -1 & 3 \end{bmatrix} \text{ and } b = \begin{bmatrix} 6 \\ -11 \\ -11 \\ 19 \end{bmatrix}\end{split}\]

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

n = 4
ncoda = 2

# Note that a is in band storage mode
a = array([[0.0, 0.0, -1.0, 1.0],
           [0.0, 0.0, 2.0, -1.0],
           [2.0, 4.0, 7.0, 3.0]])
b = array([6.0, -11.0, -11.0, 19.0])

factor = []
condition = []
x = linSolPosdefBand(a, ncoda, b)
writeMatrix("Solution, x, of Ax = b", x)

Output¶

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

Example 2¶

This example solves the same problem \(Ax = b\) given in the first example. The solution is returned and an estimate of \(\kappa_1(A)\) is computed. Additionally, the \(R^TR\) factorization is returned. Then, knowing that \(\kappa_1(A) = \|A\| \|A^{-1}\|\), the condition number is computed directly and compared to the estimate from Higham’s method.

from __future__ import print_function
from numpy import *
from pyimsl.math.linSolPosdefBand import linSolPosdefBand
from pyimsl.math.vectorNorm import vectorNorm
from pyimsl.math.writeMatrix import writeMatrix

n = 4
ncoda = 2
a = array([[0.0, 0.0, -1.0, 1.0],
           [0.0, 0.0, 2.0, -1.0],
           [2.0, 4.0, 7.0, 3.0]])
b = array([6.0, -11.0, -11.0, 19.0])

factor = []
condition = []
x1 = linSolPosdefBand(a, ncoda, b,
                      factor=factor,
                      condition=condition)
writeMatrix("Solution, x, of Ax = b", x1)

# Find one norm of inverse
inverse_norm = 0.0
e_i = empty((n), dtype=double)
for i in range(0, n):
    for j in range(0, n):
        e_i[j] = 0.0
    e_i[i] = 1.0

    # Determine one norm of each column of inverse
    x = linSolPosdefBand(a, ncoda, e_i,
                         factor=factor,
                         solveOnly=True)
    column_norm = vectorNorm(x, oneNorm=True)

    # The max of the column norms is the norm of inv(A)
    if (inverse_norm < column_norm):
        inverse_norm = column_norm

print("\nHigham's condition estimate = %f" % condition[0])
print("Direct condition estimate   = %f" % (11.0 * inverse_norm))

Output¶

Higham's condition estimate = 8.650485
Direct condition estimate   = 8.650485
 
              Solution, x, of Ax = b
          1            2            3            4
          4           -6            2            9

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.