linSolPosdefCoordinateComplex¶

Solves a sparse Hermitian positive definite system of linear equations \(Ax = b\). Using optional arguments, any of several related computations can be performed. These extra tasks include returning the symbolic factorization of A, returning the numeric factorization of A, and computing the solution of \(Ax = b\) given either the symbolic or numeric factorizations.

Synopsis¶

linSolPosdefCoordinateComplex (a, b)

Required Arguments¶

a[[]] (Input): Vector of length nz containing the location and value of each nonzero entry in lower triangle of the matrix.
complex b (Input): Vector of length n containing the right-hand side.

Return Value¶

The solution x of the sparse Hermitian positive definite linear system \(Ax = b\). If no solution was computed, then None is returned.

Optional Arguments¶

returnSymbolicFactor (Output): A structure containing, on return, the symbolic factorization of the input matrix.
supplySymbolicFactor (Input): A structure. This structure contains the symbolic factorization of the input matrix computed by linSolPosdefCoordinateComplex with the returnSymbolicFactor option. The structure is described in the returnSymbolicFactor optional argument description.
symbolicFactorOnly,: Compute the symbolic factorization of the input matrix and return. The argument b is ignored.
returnNumericFactor (Output): A structure containing, on return, the numeric factorization of the input matrix.
supplyNumericFactor, structure (Input): A structure. This structure contains the numeric factorization of the input matrix computed by linSolPosdefCoordinateComplex with the returnNumericFactor option.
numericFactorOnly,: Compute the numeric factorization of the input matrix and return. The argument b is ignored.
solveOnly,: Solve \(Ax = b\) given the numeric or symbolic factorization of A. This option requires the use of either supplyNumericFactor or supplySymbolicFactor.
multifrontalFactorization,: Perform the numeric factorization using a multifrontal technique. By default a standard factorization is computed based on a sparse compressed storage scheme.
smallestDiagonalElement(Output): A scalar containing the smallest diagonal element that occurred during the numeric factorization. This option is valid only if the numeric factorization is computed during this call to linSolPosdefCoordinateComplex.
largestDiagonalElement (Output): A scalar containing the largest diagonal element that occurred during the numeric factorization. This option is valid only if the numeric factorization is computed during this call to linSolPosdefCoordinateComplex.
numNonzerosInFactor (Output): A scalar containing the total number of nonzeros in the factor.
cscFormat, int colPtr, int rowInd, float values (Input): Accept the coefficient matrix in Compressed Sparse Column (CSC) Format. See the “Matrix Storage Modes” section of the “Introduction” at the beginning of this manual for a discussion of this storage scheme.

Description¶

The function linSolPosdefCoordinateComplex solves a system of linear algebraic equations having a sparse Hermitian positive definite coefficient matrix A. In this function’s default use, a symbolic factorization of a permutation of the coefficient matrix is computed first. Then a numerical factorization is performed. The solution of the linear system is then found using the numeric factor.

The symbolic factorization step of the computation consists of determining a minimum degree ordering and then setting up a sparse data structure for the Cholesky factor, L. This step only requires the “pattern” of the sparse coefficient matrix, i.e., the locations of the nonzeros elements but not any of the elements themselves. Thus, the value field in the sparse element array is ignored. If an application generates different sparse Hermitian positive definite coefficient matrices that all have the same sparsity pattern, then by using returnSymbolicFactor and supplySymbolicFactor, the symbolic factorization need only be computed once.

Given the sparse data structure for the Cholesky factor L, as supplied by the symbolic factor, the numeric factorization produces the entries in L so that

\[PAP^T = LL^H\]

Here P is the permutation matrix determined by the minimum degree ordering.

The numerical factorization can be carried out in one of two ways. By default, the standard factorization is performed based on a sparse compressed storage scheme. This is fully described in George and Liu (1981). Optionally, a multifrontal technique can be used. The multifrontal method requires more storage but will be faster in certain cases. The multifrontal factorization is based on the routines in Liu (1987). For a detailed description of this method, see Liu (1990), also Duff and Reid (1983, 1984), Ashcraft (1987), Ashcraft et al. (1987), and Liu (1986, 1989).

If an application requires that several linear systems be solved where the coefficient matrix is the same but the right-hand sides change, the options returnNumericFactor and supplyNumericFactor can be used to precompute the Cholesky factor. Then the solveOnly option can be used to efficiently solve all subsequent systems.

Given the numeric factorization, the solution x is obtained by the following calculations:

\[Ly_1 = Pb\]

\[L^Hy_2 = y_1\]

\[x = P^Ty_2\]

The permutation information, P, is carried in the numeric factor structure.

Examples¶

Example 1¶

As a simple example of default use, consider the following Hermitian positive definite matrix

\[\begin{split}A = \begin{bmatrix} 2 & -1+i & 0 \\ -1-i & 4 & 1+2i \\ 0 & 1-2i & 10 \\ \end{bmatrix}\end{split}\]

Let \(x^T = (1 + i, 2 + 2i, 3 + 3i)\) so that \(Ax = (-2 + 2i, 5 +15i, 36 + 28i)^T\). The number of nonzeros in the lower triangle is nz = 5.

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

a = [[0, 0, (2.0 + 0.0j)],
     [1, 1, (4.0 + 0.0j)],
     [2, 2, (10.0 + 0.0j)],
     [1, 0, (-1.0 - 1.0j)],
     [2, 1, (1.0 - 2.0j)]]
b = [-2.0 + 2.0j, 5.0 + 15.0j, 36.0 + 28.0j]
n = 3
nz = 5
x = linSolPosdefCoordinateComplex(a, b)
writeMatrixComplex("Solution, x, of Ax=b", x, column=True)

Output¶

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

Example 2¶

Set \(A = E(2500, 50)\). Then solve the system \(Ax = b_1\) and return the numeric factorization resulting from that call. Then solve the system \(Ax = b_2\) using the numeric factorization just computed. Absolute errors and execution time are printed.

from __future__ import print_function

from numpy import *
from pyimsl.math.ctime import ctime
from pyimsl.math.generateTestCoordinateComplex import generateTestCoordinateComplex
from pyimsl.math.linSolPosdefCoordinateComplex import linSolPosdefCoordinateComplex
from pyimsl.math.randomSeedSet import randomSeedSet
from pyimsl.math.randomUniform import randomUniform

ic = 50
n = ic * ic
index = 0
b_1 = empty(n, dtype=complex)
b_2 = empty(n, dtype=complex)

# Generate two right hand sides
seed = 123457
randomSeedSet(seed)
index = 0
rand_vec = randomUniform(4 * n)
for i in range(0, n):
    temp_re = rand_vec[index]
    index = index + 1
    temp_im = rand_vec[index]
    b_1[i] = complex(temp_re, temp_im)
    index = index + 1
    temp_re = rand_vec[index]
    index = index + 1
    temp_im = rand_vec[index]
    b_2[i] = complex(temp_re, temp_im)
    index = index + 1

# Build coefficient matrix a
numeric_factor = []
a = generateTestCoordinateComplex(n, ic, symmetricStorage=True)

# Now solve Ax_1 = b_1 and return the numeric factorization
time_1 = ctime()
x_1 = linSolPosdefCoordinateComplex(a, b_1,
                                    returnNumericFactor=numeric_factor)
time_1 = ctime() - time_1

# Now solve Ax_2 = b_2 given the numeric factorization
time_2 = ctime()
x_2 = linSolPosdefCoordinateComplex(a, b_2,
                                    supplyNumericFactor=numeric_factor[0],
                                    solveOnly=True)
time_2 = ctime() - time_2
try:
    ratio = time_2 / time_1
except ZeroDivisionError:
    ratio = 1.0
print("time_2/time_1 = %lf" % ratio)

Output¶

time_2/time_1 = 0.751358