Add selinv_extract: read the selected inverse at a sparse pattern

[timweiland]

Summary

Adds selinv_extract / selinv_extract! / selinv_extract_setup — a primitive that reads the selected-inverse values at a given sparse pattern straight from the supernodal blocks, without materializing the full sparse(selinv(F).Z).

It is the read analogue of the existing dot(S::SupernodalMatrix, B::SparseMatrixCSC): the same block traversal, but instead of accumulating the scalar tr(ΣB), it stores Σ's values at B's nonzero positions.

Motivation (GaussianMarkovRandomFields.jl #168): the workspace selinv path materializes the entire sparse(selinv.Z) on every refactorization, but the real consumer — diag(A Σ Aᵀ) for predictor marginals — only needs Σ at the observation-local pattern (≈ pattern(AᵀA)), a tiny subset of the factor fill. Streaming the extract (and reusing a precomputed plan) avoids that materialization and restores parallel scaling.

API

selinv_extract(S, B)                      # allocating: SparseMatrixCSC with B's pattern
selinv_extract!(dest, S, B)               # in-place into a dest carrying B's pattern
plan = selinv_extract_setup(S, B)         # precompute reuse plan (depermute baked in)
selinv_extract!(dest, S, plan)            # O(nnz), zero allocation, reusable

Semantics. Z = selinv_extract(S, B) has exactly B's sparsity pattern, and Z[i,j] == sparse(S)[i,j] for every (i,j) ∈ pattern(B) — including positions outside the stored selinv pattern, which are kept as 0.0. A generic AbstractMatrix fallback handles the simplicial case (where selinv(F).Z is a Symmetric/SparseMatrixCSC), so callers never branch on factor type.

The supernodal structure is invariant across refactorizations of the same symbolic factor (only S.vals changes), so selinv_extract_setup builds a plan once and selinv_extract!(dest, S, plan) reuses it with zero allocation.

Benchmark

benchmark/bench_extract.jl, 2D Laplacian factor n = 14400, observation-local pattern (B/fill ≈ 0.18), depermuted supernodal .Z:

method	min time	allocated
1. sparse(S) + mask	22.9 ms	97 MiB
2. per-entry getindex over S	14.8 ms	30 MiB
3. selinv_extract (allocating)	9.8 ms	12 MiB
4. selinv_extract! + plan	0.27 ms	0 B

The plan-based fill is ~85× faster than materialize-and-mask with zero allocation; even the one-shot allocating selinv_extract is ~2.3× faster while bit-identical.

Tests

test/test_extract.jl (wired into runtests.jl) covers:

• supernodal and simplicial factors,
• depermute = true / false,
• random sprand patterns and obs-style AᵀA patterns,
• exact pattern preservation including off-pattern (zero) entries,
• the allocation-free reuse path (@allocated == 0).

Full suite (including Aqua) passes locally.

Notes

• Project.toml version bumped 0.2.0 → 0.2.1 so GaussianMarkovRandomFields can depend on the new release.

🤖 Generated with Claude Code