linSolGenBand¶

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

Synopsis¶

linSolGenBand (a,  nlca, nuca, b)

Required Arguments¶

float a[[]] (Input): Array of size (nlca + nuca + 1) 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.
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)

int pPvt (Output)

An array of length n containing the pivot sequence for the factorization.

float 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 linSolGenBand is None.

solveOnly

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

blockingFactor (Input)

The blocking factor. blockingFactor must be set no larger than 32.

Default: blockingFactor = 1

Description¶

The function linSolGenBand solves a system of linear algebraic equations with a real band matrix A. It first computes the LU factorization of A based on the blocked LU factorization algorithm given in Du Croz et al. (1990). Level-3 BLAS invocations are replaced with inline loops. The blocking factor blockingFactor has the default value of 1, but can be reset to any positive value not exceeding 32.

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 linSolGenBand fails if U, the upper triangular part of the factorization, has a zero diagonal element.

Examples¶

Example 1¶

This example demonstrates the simplest use of this function by solving a system of four linear equations. The equations are as follows:

\[2x_1 − x_2 = 3\]

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

\[−x_3 + 2x_4 = 11\]

\[2x_3 + x_4 = −2\]

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

nlca = 1
nuca = 1
a = array([[0.0, -1.0, -2.0, 2.0],
           [2.0, 1.0, -1.0, 1.0],
           [-3.0, 0.0, 2.0, 0.0]])
b = [3.0, 1.0, 11.0, -2.0]
x = linSolGenBand(a, nlca, nuca, b)
writeMatrix("Solution of Ax = b", x)

Output¶

 
                Solution of Ax = b
          1            2            3            4
          2            1           -3            4

Example 2¶

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

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

factor = {}
nlca = 1
nuca = 1
a = array([[0.0, -1.0, -2.0, 2.0],
           [2.0, 1.0, -1.0, 1.0],
           [-3.0, 0.0, 2.0, 0.0]])
b = [3.0, 1.0, 11.0, -2.0]
x1 = linSolGenBand(a, nlca, nuca, b,
                   factor=factor)
writeMatrix("Solution of Ax = b", x1)
# "a" argument is ignored with solveOnly
x2 = linSolGenBand(None, nlca, nuca, b,
                   factor=factor,
                   solveOnly=True,
                   transpose=True)
writeMatrix("Solution of trans(A)x = b", x2)

Output¶

 
                Solution of Ax = b
          1            2            3            4
          2            1           -3            4
 
             Solution of trans(A)x = b
          1            2            3            4
         -6           -5           -1            0

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.