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:
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 xᵀAx ≥ 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.
Create and Factorize the Matrix
Sparse function performs the actual Cholesky factorization, finding L such that A = LLᵀ:
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
Sparse parameter with the
report parameter set to a user-supplied error-handling routine. The returned
Sparse 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:
The solve call takes the factorization A = LLᵀ and solves the system Ax = b as LLᵀx = 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.
Sparse solves the equation, populating
x with the solution:
You can then iterate over the solution vector,
x, and print the solution.
You should get this output:
x = 0.