sparseCholeskySmp¶

Computes the Cholesky factorization of a sparse real symmetric positive definite matrix A by an OpenMP parallelized supernodal algorithm and solves the sparse real positive definite system of linear equations Ax = b.

Synopsis¶

sparseCholeskySmp (a, b)

Required Arguments¶

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

Return Value¶

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

Optional Arguments¶

cscFormat, int colPtr[], int rowInd[], float values[] (Input)

Accept the coefficient matrix in compressed sparse column (CSC) format, as described in the Compressed Sparse Column (CSC) Format section of the “Introduction” chapter of this manual.

preordering, int (Input)

The variant of the Minimum Degree Ordering (MDO) algorithm used in the preordering of matrix A:

`preordering`	Method
0	George and Liu’s Quotient Minimum Degree algorithm.
1	A variant of George and Liu’s Quotient Minimum Degree algorithm using a preprocessing phase and external degrees.

Default: preordering = 0.

returnSymbolicFactor (Output)

A structure containing, on return, the supernodal symbolic factorization of the input matrix.

supplySymbolicFactor, structure (Input)

A structure. This structure contains the symbolic factorization of the input matrix computed by sparseCholeskySmp with the returnSymbolicFactor option.

symbolicFactorOnly, (Input)

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 (Input)

A structure. This structure contains the numeric factorization of the input matrix computed by sparseCholeskySmp with the returnNumericFactor option.

numericFactorOnly, (Input)

Compute the numeric factorization of the input matrix and return. The argument b is ignored.

solveOnly, (Input)

Solve \(Ax = b\) given the numeric or symbolic factorization of A. This option requires the use of either supplyNumericFactor or supplySymbolicFactor.

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 sparseCholeskySmp.

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 sparseCholeskySmp.

numNonzerosInFactor (Output)

A scalar containing the total number of nonzeros in the factor.

Description¶

The function sparseCholeskySmp solves a system of linear algebraic equations having a sparse symmetric positive definite coefficient matrix A. In this function’s default usage, a symbolic factorization of a permutation of the coefficient matrix is computed first. Then a numerical factorization exploiting OpenMP parallelism 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 supernodal 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 symmetric positive definite coefficient matrices that all have the same sparsity pattern, then by using returnSymbolicFactor and supplySymbolicFactor, the symbolic factorization needs 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^T\]

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

The numerical factorization is an implementation of a parallel supernodal algorithm that combines a left-looking and a right-looking column computation scheme. This algorithm is described in detail in Rauber et al. (1999).

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^Ty_2 = y_1\]

\[x = P^Ty_2\]

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

Examples¶

Example 1¶

Consider the 5 × 5 coefficient matrix A,

\[\begin{split}A = \begin{bmatrix} 10 & 0 & 1 & 0 & 2 \\ 0 & 20 & 0 & 0 & 3 \\ 1 & 0 & 30 & 4 & 0 \\ 0 & 0 & 4 & 40 & 5 \\ 2 & 3 & 0 & 5 & 50 \\ \end{bmatrix}\end{split}\]

The number of nonzeros in the lower triangle of A is nz = 10. We construct the solution \(x^T = (5, 4, 3, 2, 1)\) to the system \(Ax = b\) by setting \(b := Ax = (55, 83, 103, 97, 82)^T\). The solution is computed and printed.

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

a = [[0, 0, 10.0],
     [1, 1, 20.0],
     [2, 0, 1.0],
     [2, 2, 30.0],
     [3, 2, 4.0],
     [3, 3, 40.0],
     [4, 0, 2.0],
     [4, 1, 3.0],
     [4, 3, 5.0],
     [4, 4, 50.0]]
b = [55.0, 83.0, 103.0, 97.0, 82.0]
x = sparseCholeskySmp(a, b)
writeMatrix("solution", x)

Output¶

 
                           solution
          1            2            3            4            5
          5            4            3            2            1

Example 2¶

This example shows how a symbolic factor can be re-used. At first, the system \(Ax = b\) with \(A = E(2500, 50)\) is solved and the symbolic factorization of A is returned. Then, the system \(Cy = d\) with \(C = A+2*I\), I the identity matrix, is solved using the symbolic factorization just computed. This is possible because A and C have the same nonzero structure and therefore also the same symbolic factorization. The solution errors are printed.

from __future__ import print_function
from numpy import *
from pyimsl.math.ctime import ctime
from pyimsl.math.generateTestCoordinate import generateTestCoordinate
from pyimsl.math.freeSnodalSymbolicFactor import freeSnodalSymbolicFactor
from pyimsl.math.vectorNorm import vectorNorm
from pyimsl.math.matMulRectCoordinate import matMulRectCoordinate
from pyimsl.math.sparseCholeskySmp import sparseCholeskySmp

ic = 50
n = ic * ic

# Build coefficient matrix a
a = generateTestCoordinate(n, ic, symmetricStorage=True)

# Build coefficient matrix C
c = a
for i in range(n):
    c[i] = [a[i][0], a[i][1], 6.0]

# Form right hand side b
mod_vector = empty((n), dtype=double)
for i in range(0, n, 1):
    mod_vector[i] = (float)(i % 5)

b = matMulRectCoordinate("A*x",
                         aMatrix={'nrowa': n, 'ncola': n, 'a': a},
                         xVector=mod_vector,
                         symmetricStorage=True)

# Form right hand side d
d = matMulRectCoordinate("A*x",
                         aMatrix={'nrowa': n, 'ncola': n, 'a': c},
                         xVector=mod_vector,
                         symmetricStorage=True)

# Solve Ax = b and return the symbolic factorization
symbolic_factor = []
x = sparseCholeskySmp(a, b,
                      returnSymbolicFactor=symbolic_factor)

# Compute solution error |x - mod_vector|
index = []
error_1 = vectorNorm(x, secondVector=mod_vector, infNorm=index)

# Solve Cy = d given the symbolic factorization
y = sparseCholeskySmp(c, d,
                      supplySymbolicFactor=symbolic_factor[0])

# Compute solution error |y - mod_vector|
error_2 = vectorNorm(y, secondVector=mod_vector, infNorm=index)

print("Solution error |x - mod_vector| = %e" % error_1)
print("Solution error |y - mod_vector| = %e" % error_2)

freeSnodalSymbolicFactor(symbolic_factor[0])

Output¶

Solution error |x - mod_vector| = 8.437695e-15
Solution error |y - mod_vector| = 8.437695e-15

Example 3¶

In this example, set \(A = E(2500, 50)\). First 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. The ratio of execution times is printed. Be aware that timing results are highly machine dependent.

from __future__ import print_function
from numpy import *
from pyimsl.math.ctime import ctime
from pyimsl.math.generateTestCoordinate import generateTestCoordinate
from pyimsl.math.randomSeedSet import randomSeedSet
from pyimsl.math.randomUniform import randomUniform
from pyimsl.math.sparseCholeskySmp import sparseCholeskySmp
from pyimsl.math.freeNumericFactor import freeNumericFactor

ic = 50
n = ic * ic
# Generate two right hand sides
seed = 123457
randomSeedSet(seed)
b_1 = randomUniform(n)
randomSeedSet(seed)
b_2 = randomUniform(n)

# Build coefficient matrix a
a = generateTestCoordinate(n, ic, symmetricStorage=True)

# Now solve Ax_1 = b_1 and return the numeric factorization
time_1 = ctime()
numeric_factor = []
x_1 = sparseCholeskySmp(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 = sparseCholeskySmp(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\n" % ratio)

freeNumericFactor(numeric_factor[0])

Output¶

time_2/time_1 = 0.461442

Fatal Errors¶

IMSL_BAD_SQUARE_ROOT A zero or negative square root has occurred during the factorization. The coefficient matrix is not positive definite.