linSolGenBandComplex¶

Solves a complex general band 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, solving \(A^Hx = b\), or computing the solution of \(Ax = b\) given the LU factorization of A.

Synopsis¶

linSolGenBandComplex (a, nlca, nuca, b)

Required Arguments¶

complex a[[]] (Input): Array of size (nlca + nuca + 1) × n containing the n × n banded coefficient matrix in band storage mode.
int nlca (Input): Number of lower codiagonals in a.
int nuca (Input): Number of upper codiagonals in a.
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, None is returned.

Optional Arguments¶

transpose

Solve \(A^Hx = 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 (2nlca + nuca + 1) × n containing the LU factorization of A with column pivoting.

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 linSolGenBandComplex is None.
solveOnly: Solve \(Ax = b\) given the LU factorization previously computed by linSolGenBandComplex. By default, the solution to \(Ax = b\) is pointed to by linSolGenBandComplex. If solveOnly is used, argument factor is required and argument a is ignored.

Description¶

The function linSolGenBandComplex solves a system of linear algebraic equations with a complex band matrix A. It first computes the LU factorization of A using scaled partial pivoting. Scaled partial pivoting differs from partial pivoting in that the pivoting strategy is the same as if each row were scaled to have the same \(L_∞\) norm. The factorization fails if U has a zero diagonal element. This can occur only if A is singular or very close to a singular matrix.

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 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 linSolGenBandComplex fails if U, the upper triangular part of the factorization, has a zero diagonal element. The function linSolGenBandComplex is based on the LINPACK subroutine CGBFA; see Dongarra et al. (1979). CGBFA uses unscaled partial pivoting.

Examples¶

Example 1¶

The following linear system is solved:

\[\begin{split}\begin{bmatrix} -2-3i & 4 & 0 & 0 \\ 6+i & -0.5+3i & -2+2i & 0 \\ 0 & 1+i & 3-3i & -4-1 \\ 0 & 0 & 2i & 1-i \end{bmatrix} \begin{bmatrix} x_0 \\ x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} -10-5i \\ 9.5+5.5i \\ 12-12i \\ 8i \end{bmatrix}\end{split}\]

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

nlca = 1
nuca = 1
factor = {}

# Note that a is in band storage mode
a = array([[0.0 + 0.0j, 4.0 + 0.0j, -2.0 + 2.0j, -4.0 - 1.0j],
           [-2.0 - 3.0j, -0.5 + 3.0j, 3.0 - 3.0j, 1.0 - 1.0j],
           [6.0 + 1.0j, 1.0 + 1.0j, 0.0 + 2.0j, 0.0 + 0.0j]])
b = array([-10.0 - 5.0j, 9.5 + 5.5j, 12.0 - 12.0j, 0.0 + 8.0j])
x = linSolGenBandComplex(a, nlca, nuca, b, factor=factor)
writeMatrixComplex("Solution, x, of Ax = b", x, column=True)

Output¶

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

Example 2¶

This example solves the problem \(Ax = b\) using the data from the first example. This time, the factorizations are returned and then the problem \(A^Hx = b\) is solved without recomputing LU.

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

nlca = 1
nuca = 1
factor = {}

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

# Solve Ax=b and return LU
x1 = linSolGenBandComplex(a, nlca, nuca, b, factor=factor)

# Use precomputed LU to solve ctrans(A)x = b
writeMatrixComplex("solution of Ax = b", x1, column=True)
x2 = linSolGenBandComplex(a, nlca, nuca, b,
                          factor=factor, transpose=True)
writeMatrixComplex("solution of ctrans(A)x = b", x2, column=True)

Output¶

 
     solution of Ax = b
(          3,         -0)
(         -1,          1)
(          3,         -0)
(         -1,          1)
 
 solution of ctrans(A)x = b
(       5.58,      -2.91)
(      -0.48,      -4.67)
(      -6.19,       7.15)
(      12.60,      30.20)

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.