linSolPosdefBandComplex¶

Solves a complex Hermitian 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^HR\) 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¶

linSolPosdefBandComplex (a, ncoda, b)

Required Arguments¶

complex 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.
complex 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 \(R^HR\) factorization of 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 \(R^HR\) factorization of A. If factorOnly is used, either factor is required. The argument b is then ignored, and the returned value of linSolPosdefBandComplex is None.
solveOnly: Solve \(Ax = b\) given the \(R^HR\) factorization previously computed by linSolPosdefBandComplex. By default, the solution to \(Ax = b\) is pointed to by linSolPosdefBandComplex. If solveOnly is used, argument factor is required and the argument a is ignored.

Description¶

The function linSolPosdefBandComplex solves a system of linear algebraic equations with a real symmetric positive definite band coefficient matrix A. It computes the \(R^HR\) Cholesky factorization of A. Argument 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 linSolPosdefBandComplex 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 linSolPosdefBandComplex is based partially on the LINPACK subroutines SPBFA and CPBSL; see Dongarra et al. (1979).

Examples¶

Example 1¶

Solve a linear system \(Ax = b\) where

\[\begin{split}A = \begin{bmatrix} 2 & -1+i & 0 & 0 & 0 \\ -1-i & 4 & 1+2i & 0 & 0 \\ 0 & 1-2i & 10 & 4i & 0 \\ 0 & 0 & -4i & 6 & 1+i \\ 0 & 0 & 0 & 1-i & 9 \end{bmatrix} \text{ and } b = \begin{bmatrix} 1+5i \\ 12-6i \\ 1-16i \\ -3-3i \\ 25+16i \end{bmatrix}\end{split}\]

from numpy import *
from pyimsl.math.linSolPosdefBandComplex import linSolPosdefBandComplex
from pyimsl.math.writeMatrixComplex import writeMatrixComplex

n = 5
ncoda = 1

# Note that a is in band storage mode
a = array([[0.0 + 0.0j, -1.0 + 1.0j, 1.0 + 2.0j, 0.0 + 4.0j, 1.0 + 1.0j],
           [2.0 + 0.0j, 4.0 + 0.0j, 10.0 + 0.0j, 6.0 + 0.0j, 9.0 + 0.0j]])
b = array([1.0 + 5.0j, 12.0 - 6.0j, 1.0 - 16.0j, -3.0 + -3.0j, 25.0 + 16.0j])
factor = []
condition = []
x = linSolPosdefBandComplex(a, ncoda, b)
writeMatrixComplex("Solution, x, of Ax = b", x, column=True)

Output¶

 
   Solution, x, of Ax = b
(          2,          1)
(          3,          0)
(         -1,         -1)
(          0,         -2)
(          3,          2)

Example 2¶

This example solves the same problem \(Ax = b\) given in the first example. The solution is returned in user-allocated space and an estimate of \(\kappa_1(A)\) is computed. Additionally, the \(R^HR\) 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.linSolPosdefBandComplex import linSolPosdefBandComplex
from pyimsl.math.vectorNormComplex import vectorNormComplex
from pyimsl.math.writeMatrixComplex import writeMatrixComplex

n = 5
ncoda = 1

# Note that a is in band storage mode
a = array([[0.0 + 0.0j, -1.0 + 1.0j, 1.0 + 2.0j, 0.0 + 4.0j, 1.0 + 1.0j],
           [2.0 + 0.0j, 4.0 + 0.0j, 10.0 + 0.0j, 6.0 + 0.0j, 9.0 + 0.0j]])
b = array([1.0 + 5.0j, 12.0 - 6.0j, 1.0 - 16.0j, -3.0 - 3.0j, 25.0 + 16.0j])

factor = []
condition = []
x1 = linSolPosdefBandComplex(a, ncoda, b,
                             factor=factor,
                             condition=condition)
writeMatrixComplex("Solution, x, of Ax = b", x1, column=True)

# Find one norm of inverse
inverse_norm = 0.0
e_i = empty((n), dtype=complex)
for i in range(0, n, 1):
    for j in range(0, n, 1):
        e_i[j] = (0.0 + 0.0j)
    e_i[i] = (1.0 + 0.0j)
    x = linSolPosdefBandComplex(a, ncoda, e_i,
                                factor=factor,
                                solveOnly=True)
    column_norm = vectorNormComplex(x, oneNorm=True)
    if (inverse_norm < column_norm):
        inverse_norm = column_norm
print("\nHigham's condition estimate = %f" % condition[0])
print("Direct condition estimate   = %f" % ((14.0 + sqrt(5.0)) * inverse_norm))

Output¶

Higham's condition estimate = 19.377722
Direct condition estimate   = 19.377722
 
   Solution, x, of Ax = b
1  (          2,          1)
2  (          3,          0)
3  (         -1,         -1)
4  (          0,         -2)
5  (          3,          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_NONPOSITIVE_MATRIX`	The leading # by # submatrix of the input matrix is not positive definite.
`IMSL_SINGULAR_MATRIX`	The input matrix is singular.