Solving Systems Using Direct Methods

Use direct methods to solve systems of equations where the coefficient matrix is sparse.


Direct methods offer high-precision solving with a simple API when compared to iterative methods (see Implementing Iterative Methods for a detailed comparison). The code discussed here uses sparse Cholesky factorization to solve the following equation:

Equation to solve using sparse Cholesky factorization.

In this illustration, A refers to the 4 x 4 matrix, and b to the right-hand-side vector. The code in this article solves the equation Ax = b by finding x.

Note that A is sparse: Some entries (those that are blank) are zero. For small matrices such as this, there is little gain in exploiting this structure. However it is essential for larger systems that would not otherwise fit in memory, even on a large machine.

The code below performs a sparse Cholesky factorization, equivalent to calling the LAPACK function DPOTRF on a dense matrix. The main requirement for sparse Cholesky factorization is that the matrix is symmetric positive-definite (that is, A=Aᵀ), and xAx 0 for all x (that is, all eigenvalues are greater than zero). A sufficient (but not necessary) condition is that the matrix is diagonally dominant (the sum of the absolute values of the off-diagonal entries in each row or column is less than the value of the diagonal). This is the case for the above matrix.

Create the Matrix Structure

Use the code below to define the matrix structure. Recall from Creating Sparse Matrices that because A is symmetric, it stores only half of the data.

int rowIndices[]    = {0, 1, 3, 1, 2, 3, 2, 3}; 
double values[] = {10.0, 1.0, 2.5, 12.0, -0.3, 1.1, 9.5, 6.0};
long columnStarts[] = { 0, 3, 6, 7, 8}; 
SparseMatrixStructure structure = {
    .rowCount     = 4,
    .columnCount  = 4,
    .columnStarts = columnStarts,
    .rowIndices   = rowIndices,
    .attributes = {
        .kind = SparseSymmetric,
        .triangle = SparseLowerTriangle
    .blockSize = 1

Create and Factorize the Matrix

The SparseFactor function performs the actual Cholesky factorization, finding L such that A = LL:

var values = [10.0, 1.0, 2.5,           // Column 0
              12.0, -0.3, 1.1,          // Column 1
              9.5,                      // Column 2
              6.0]                      // Column 3

SparseMatrix_Double A = {
    .structure = structure,
    .data = values

__auto_type LLT = SparseFactor(SparseFactorizationCholesky, A);

If there is an error, the above code prints an error message and terminates. You may instead want to capture the error by using the optional SparseSymbolicFactorOptions parameter with the reportError parameter set to a user-supplied error-handling routine. The returned SparseOpaqueFactorization_Double structure reflects the error.

Solve the Equation

Once you have the factorization, you can solve the equation. The right-hand-side and solution vectors are arrays wrapped in flexible data types. The actual values of x do not matter, because they are overwritten by the solution:

double bValues[] = { 2.20, 2.85, 2.79, 2.87 };
DenseVector_Double b = {
    .count = 4,     
    .data = bValues 
double xValues[] = { 0.00, 0.00, 0.00, 0.00 };
DenseVector_Double x = {
    .count = 4,    
    .data = xValues 

The solve call takes the factorization A = LL and solves the system Ax = b as LLx = b by solving the two triangular systems:

  • Ly = b

  • Lᵀx = y

However, you need only to supply the right-hand-side vector and the factorization.

SparseSolve solves the equation, populating x with the solution:

SparseSolve(LLT, b, x);

Again, if there is an error, the code prints an error message and terminates, unless you set reportError on the initial call to SparseFactor.

You can then iterate over the solution vector, x, and print the solution.

printf("x = "); for(int i=0; i<x.count; i++) printf(" %.2f",[i]); printf("\n");

You should get this output: x = 0.10 0.20 0.30 0.40.

See Also

Sparse Matrices

Creating Sparse Matrices

Create sparse matrices for factorization and solving systems.

Implementing Iterative Methods

Use iterative methods to solve large problems faster and with a lower memory overhead than with direct methods.

Solving Systems Using Iterative Methods

Use iterative methods to solve systems of equations where the coefficient matrix is sparse.

Creating a Sparse Matrix from Coordinate Format Arrays

Use separate coordinate format arrays to create sparse matrices.

Sparse Solvers

Solve systems of equations where the coefficient matrix is sparse.