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

Framework

- Accelerate

## Overview

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**ᵀ**A**x **≥** 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.

### Factorize the Matrix

The `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:

Again, if there is an error, the code prints an error message and terminates, unless you set `report`

on the initial call to `Sparse`

.

You can then iterate over the solution vector, `x`

, and print the solution.

You should get this output: `x = 0`

.