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

Use the code below to define the matrix A. 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
SparseMatrix_Double A = {
    .structure = structure,
    .data = values

Factorize the Matrix

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

__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

Direct Sparse Solving Methods

Sparse Factor Functions

Compute the factorization of a matrix.

Sparse Direct Solving Functions (Matrix RHS)

Solve a system using a factored matrix.

Advanced Functions for Sparse Solving

Work with symbolic factorizations, refactoring, and subfactors.