Module sparse

Sparse matrices and the preconditioned conjugate-gradient solve.

Mesh operators (graph and cotangent Laplacians, their squares, mass matrices) have a handful of non-zeros per row, so a dense Matrix of a few thousand mesh vertices already costs more memory than the mesh itself and an O(n^3) factorisation more time than the rest of the algorithm. SparseMatrix stores only the non-zeros in compressed rows, and conjugate_gradient solves a symmetric positive definite system with it in O(iterations * non-zeros).

The solve follows Hestenes and Stiefel (1952) with a Jacobi (diagonal) preconditioner, the textbook form in Saad, Iterative Methods for Sparse Linear Systems, algorithm 9.1. Its stopping test is on the true residual b - A x, recomputed from the matrix when the recurrence says the tolerance is met, so drift in the recurred residual cannot report a convergence that did not happen.

Structs

ConjugateGradientOptions

Stopping rule and budget for conjugate_gradient.

ConjugateGradientSolution

The result of conjugate_gradient.

SparseMatrix

A sparse matrix in compressed sparse row form.

Functions

conjugate_gradient

Solve A x = b for a symmetric positive definite sparse A by the Jacobi-preconditioned conjugate-gradient method, starting from initial (zero when None).