Sparse Matrices And Solvers

Use this guide for Eigen/src/SparseCore, the sparse decompositions in SparseCholesky, SparseLU, SparseQR, and OrderingMethods, the external *Support backend wrappers, and their tests. IterativeLinearSolvers shares the solver contract below; numerics.md governs accuracy expectations for all of them. The checked-out headers are authoritative:

Compressed And Uncompressed Storage

A SparseMatrix is in one of two storage modes, and most bugs in this module come from a path that silently assumes one. m_innerNonZeros == nullptr means compressed; when it is non-null, inner vector j occupies [outerIndexPtr()[j], outerIndexPtr()[j] + innerNonZeroPtr()[j]) rather than running to outerIndexPtr()[j + 1].

  • insert() and coeffRef() turn a compressed matrix into uncompressed mode when they add an entry. A function that takes a SparseMatrix& and inserts is therefore free to change the caller's storage mode; makeCompressed() restores it.
  • Do not derive an entry count or an iteration bound from consecutive outerIndexPtr() differences. Use nonZeros(), innerNonZeroPtr(), or the uniform loop that SparsityPatternRef.h documents, which is correct in both modes.
  • resize() zeroes the matrix, drops to compressed mode, and keeps the allocation; conservativeResize() preserves contents. Neither is a way to change storage mode deliberately.
  • Ref<SparseMatrix> accepts an uncompressed argument unless it is declared with StandardCompressedFormat. With that option a writable Ref asserts isCompressed(), while a Ref<const SparseMatrix, StandardCompressedFormat> silently materializes a compressed copy instead of failing. A Ref parameter is consequently not proof that no copy happened; state which form a new API takes and why.
  • InnerIterator and every raw pointer obtained from the matrix are invalidated by an insertion. Finish iterating, or collect the coordinates first and mutate afterwards.
  • A consumer may require compressed input outright rather than adapting: SparseQR::analyzePattern opens with an eigen_assert(mat.isCompressed()), while SparseLU branches on isCompressed() and falls back to copying the outer index array. Because that check is an eigen_assert, a release build does not diagnose the violation at all — call makeCompressed() before handing an assembled matrix to a direct solver rather than relying on the assertion.

Sorted Inner Indices

Sorted inner indices within each inner vector are an invariant most of the API maintains and several parts require. coeff() resolves an entry by binary search over the inner range, and the sparse-sparse coefficient-wise binary evaluator merges two inner streams by advancing whichever index is smaller. Both return wrong values on unsorted input rather than failing.

The invariant is not universal, and the exception is a trap. SparseQR::matrixR() hands back a reference to a stored factor built with insertBackByOuterInnerUnordered, so it is compressed but not sorted, and the rank-deficient path preserves that: right-multiplying a column-major matrix by the pivot permutation takes the outer-permutation branch, which relocates whole inner vectors without reordering within them. Compressed is not sorted, and a factor returned by reference has had no opportunity to become either.

A storage-order round-trip is the fix, and it works for a reason worth knowing: a cross-order sparse assignment is a counting transpose that walks the source in outer order, so each destination inner vector receives its entries in ascending order no matter how the source was ordered. Eigen relies on this internally — SparseQR::_sort_matrix_Q() sorts the stored reflectors exactly that way before matrixQ() is materialized into a sparse destination, which is why Q is not a hazard in the way R is. Verify with innerIndicesAreSorted() rather than assuming a decomposition's output factor satisfies the invariant.

  • setFromTriplets() accepts unsorted input with duplicates and produces a sorted, compressed matrix with duplicates summed. It destroys the previous contents and does not resize — construct or resize() the matrix first, since the dimensions are not inferred from the triplets.
  • setFromSortedTriplets(), insertFromTriplets(), and insertFromSortedTriplets() are the other three corners: the Sorted variants promise pre-sorted input, and the insertFrom variants merge into existing entries rather than replacing them. All four take an optional duplicate functor; the default sums.
  • insert() requires that the entry not already exist. Use coeffRef() when it may, and reserve(const SizesType&) before random-order insertion — the sequential fast path only holds for increasing outer indices.
  • The sparse-sparse product selectors preserve sortedness deliberately, choosing between a sorted insertion and an unsorted pass followed by a transpose round-trip that sorts as a side effect (ConservativeSparseSparseProduct.h). A new product, permutation, or assembly path must restore the invariant; sortInnerIndices() and innerIndicesAreSorted() on SparseCompressedBase are the tools, and the latter belongs in a test rather than only in reasoning.

Products

A * B on two sparse operands uses the conservative product; (A * B).pruned() selects the pruning product in SparseSparseProductWithPruning.h instead. The distinction is semantic before it is a matter of performance: the conservative path stores every structurally generated entry, including one whose accumulation cancels to exactly zero, while the pruning path drops completed values at or below its tolerance. The two therefore produce different patterns from the same operands, and a test or benchmark written against one does not transfer to the other. Neither reserves the exact result size up front — the conservative path starts from the heuristic nonZerosEstimate() sum documented at its definition and grows or over-allocates from there.

Threaded SpMV is opt-in: Eigen/SparseCore includes ThreadedSparseProduct.h and Eigen/ThreadPool only under EIGEN_USE_THREADS. Coverage lives in test/sparse_threaded_product.cpp, and tensor-threadpool.md applies to its threading.

Solver Contract

Direct sparse solvers split pattern analysis from numerical work: analyzePattern(), then factorize(), with compute() doing both. Re-solving with the same pattern and new values must reuse the analysis; a change that forces a re-analysis is a performance regression even when results match.

  • Sparse solvers report failure through info(), not exceptions, and info() lives on each concrete solver and on IterativeSolverBase — not on SparseSolverBase. Check it after compute()/factorize() and again after solve() where the solver documents doing so. A test that ignores info() can pass on a matrix the solver rejected.
  • solve() asserts that the solver was initialized, so a missing compute() surfaces only in a debug build.
  • Reordering is part of the result: a solver's permutation affects fill-in and the achievable accuracy, so an ordering change needs the fill-in or timing evidence benchmarking.md asks for, not only a residual check.

Testing Sparse Changes

initSparse() in test/sparse.h fills a dense reference and a sparse matrix together, with ForceNonZeroDiag, MakeLowerTriangular, MakeUpperTriangular, and ForceRealDiag for the shapes solvers require. test/sparse_solver.h provides the check_sparse_solving, check_sparse_spd_solving, check_sparse_nonhermitian_solving, and determinant harnesses; prefer them to a hand-rolled solve so a new solver inherits the established coverage.

Scale coverage to the axes this module actually branches on: both storage orders, both storage modes, a non-default StorageIndex width, complex scalars where conjugation is not a no-op, and a matrix with an empty inner vector. Comparing against a dense reference computed by Eigen is the standard technique; keep the tolerance a named epsilon multiple scaled by dimension or conditioning as testing.md requires.

External backend tests are registered conditionally in test/CMakeLists.txt, so a green local run says nothing about any of them. The full set is cholmod_support, umfpack_support, klu_support, superlu_support, pastix_support, spqr_support, accelerate_support, metis_support — an ordering backend rather than a solver — and pardiso_support. Most are gated on a find_package result and register in EIGEN_MISSING_BACKENDS when absent, several additionally require EIGEN_BUILD_BLAS or EIGEN_BUILD_LAPACK, and metis_support and pastix_support depend on variables the PaStiX search sets when the METIS component is requested.

pardiso_support is the exception worth knowing: the tree contains no find_package(PARDISO) and no EIGEN_MISSING_BACKENDS entry for it, so it is registered only when PARDISO_FOUND arrives from outside the project and its absence is silent even in the missing-backend summary. Report which sparse backends were unavailable rather than implying full coverage.