linSolGenComplex¶

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

Synopsis¶

linSolGenComplex (a, b)

Required Arguments¶

complex a[[]] (Input): Array of size n × n containing the matrix.
complex b[] (Input): Array of length 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^Hx = b\)

Default: Solve \(Ax = b\)

factor, int pPvt, complex pFactor (Output)

int pPvt (Output)

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

complex pFactor (Output)

An array of size n × n containing the LU factorization of A with column pivoting. The lower-triangular part of this array contains information necessary to construct L, and the upper-triangular part contains U.

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 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 linSolGenComplex is None.
solveOnly: Solve \(Ax = b\) given the LU factorization previously computed by linSolGenComplex. By default, the solution to \(Ax = b\) is pointed to by linSolGenComplex. If solveOnly is used, argument factor is required and argument a is ignored.
inverseOnly: Compute the inverse of the matrix A. If inverseOnly is used, inverse is required. Argument b is then ignored, and the returned value of linSolGenComplex is None.

Description¶

The function linSolGenComplex solves a system of linear algebraic equations with a complex coefficient matrix A. It first computes the LU factorization of A with partial pivoting such that \(L^{-1}A = U\). Let F be the matrix pFactor returned by optional argument factor. The triangular matrix U is stored in the upper triangle of F. The strict lower triangle of F contains the information needed to reconstruct \(\text{L}^{‑1}\) using

\[L^{-1} = L_{n-1}P_{n-1} \ldots L_1P_1\]

The factors \(P_i\) and \(L_i\) are defined by partial pivoting. \(P_i\) is the identity matrix with rows i and pPvt[i-1] interchanged. \(L_i\) is the identity matrix with \(F_{ji}\) , for \(j = i + 1,…, n\), inserted below the diagonal in column i.

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

Examples¶

Example 1¶

This example solves a system of three linear equations. The equations are:

\[(1 + i) x_1 + (2 + 3i) x^2 + (3 − 3i) x_3 = 3 + 5i\]

\[(2 + i) x_1 + (5 + 3i) x_2 + (7 − 5i) x_3 = 22 + 10i\]

\[(−2 + i) x_1 + (−4 + 4i) x_2 + (5 + 3i) x_3 = −10 + 4i\]

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

factor = {}
a = array([[1.0 + 1.0j, 2.0 + 3.0j, 3.0 - 3.0j],
           [2.0 + 1.0j, 5.0 + 3.0j, 7.0 - 5.0j],
           [-2.0 + 1.0j, -4.0 + 4.0j, 5.0 + 3.0j]])
b = array([3.0 + 5.0j, 22.0 + 10.0j, -10.0 + 4.0j])
x = linSolGenComplex(a, b)
writeMatrixComplex("Solution, x, of Ax = b", x, column=True)

Output¶

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

Example 2¶

This example solves the conjugate transpose problem \(A^Hx = b\) and returns the LU factorization of A using partial pivoting. This example differs from the first example in that the solution array is allocated in the main program.

from numpy import *
from pyimsl.math.linSolGenComplex import linSolGenComplex
from pyimsl.math.writeMatrix import writeMatrix
from pyimsl.math.writeMatrixComplex import writeMatrixComplex

factor = {}
a = array([[1.0 + 1.0j, 2.0 + 3.0j, 3.0 - 3.0j],
           [2.0 + 1.0j, 5.0 + 3.0j, 7.0 - 5.0j],
           [-2.0 + 1.0j, -4.0 + 4.0j, 5.0 + 3.0j]])
b = array([3.0 + 5.0j, 22.0 + 10.0j, -10.0 + 4.0j])

# Solve ctrans(A)*x = b for x
x = linSolGenComplex(a, b, transpose=True, factor=factor)
writeMatrixComplex("Solution, x, of ctrans(A)x = b", x, column=True)
writeMatrixComplex("LU factors of A", factor['pFactor'])
writeMatrix("Pivot sequence", factor['pPvt'])

Output¶

 
Solution, x, of ctrans(A)x = b
 1  (      -9.79,      11.23)
 2  (       2.96,      -3.13)
 3  (       1.85,       2.47)
 
                    LU factors of A
                           1                          2
1  (     -2.000,      1.000)  (     -4.000,      4.000)
2  (      0.600,      0.800)  (     -1.200,      1.400)
3  (      0.200,      0.600)  (     -1.118,      0.529)
 
                           3
1  (      5.000,      3.000)
2  (      2.200,      0.600)
3  (      4.824,      1.294)
 
           Pivot sequence
          1            2            3
          3            3            3

Example 3¶

This example computes the inverse of the 3 × 3 matrix A in the first example and also solves the linear system. The product matrix \(C = A^{-1}A\) is computed as a check. The approximate result is \(C = I\).

from numpy import *
from pyimsl.math.linSolGenComplex import linSolGenComplex
from pyimsl.math.matMulRectComplex import matMulRectComplex
from pyimsl.math.writeMatrixComplex import writeMatrixComplex

inv = []
a = array([[1.0 + 1.0j, 2.0 + 3.0j, 3.0 - 3.0j],
           [2.0 + 1.0j, 5.0 + 3.0j, 7.0 - 5.0j],
           [-2.0 + 1.0j, -4.0 + 4.0j, 5.0 + 3.0j]])
b = array([3.0 + 5.0j, 22.0 + 10.0j, -10.0 + 4.0j])

# Solve Ax=b for x
x = linSolGenComplex(a, b, inverse=inv)
writeMatrixComplex("Solution, x, of Ax = b", x, column=True)
writeMatrixComplex("Input A", a)
writeMatrixComplex("Inverse of A", inv)
c = matMulRectComplex("A*B",
                      aMatrix=inv,
                      bMatrix=a)
writeMatrixComplex("Product, inv(A)*A", c)

Output¶

 
   Solution, x, of Ax = b
1  (          1,         -1)
2  (          2,          4)
3  (          3,          0)
 
                        Input A
                           1                          2
1  (          1,          1)  (          2,          3)
2  (          2,          1)  (          5,          3)
3  (         -2,          1)  (         -4,          4)
 
                           3
1  (          3,         -3)
2  (          7,         -5)
3  (          5,          3)
 
                     Inverse of A
                           1                          2
1  (      1.330,      0.594)  (     -0.151,      0.028)
2  (     -0.632,     -0.538)  (      0.160,      0.189)
3  (     -0.189,      0.160)  (      0.193,     -0.052)
 
                           3
1  (     -0.604,      0.613)
2  (      0.142,     -0.245)
3  (      0.024,      0.042)
 
                   Product, inv(A)*A
                           1                          2
1  (          1,         -0)  (         -0,          0)
2  (          0,          0)  (          1,          0)
3  (         -0,         -0)  (         -0,         -0)
 
                           3
1  (          0,          0)
2  (         -0,         -0)
3  (          1,          0)

Warning Errors¶

IMSL_ILL_CONDITIONED The input matrix is too ill-conditioned. An estimate of the reciprocal of the \(L_1\) condition number is “rcond” = #. The solution might not be accurate.

Fatal Errors¶

IMSL_SINGULAR_MATRIX The input matrix is singular.