MeanField_SERiF/docs/sparse_direct.md

# Sparse direct preconditioner backend

`strategies::inverse::SparseDirect` factors an assembled preconditioner block once per Newton state with a
pivoted sparse LU and applies it with one exact solve. `make_primal_preconditioner(layout, SparseDirect{})`
uses it for the gravitational potential; `complete_sandbox` and the profiling drivers use it by default
when it is available (`--potential-inverse auto`).

The factorization library is chosen when SERiF is configured:

```sh
meson setup build -Dsparse_direct=auto        # default
meson setup build -Dsparse_direct=accelerate  # Apple Accelerate (macOS)
meson setup build -Dsparse_direct=umfpack     # SuiteSparse UMFPACK (any platform)
meson setup build -Dsparse_direct=none        # no backend
```

| Value | Library | Notes |
|---|---|---|
| `auto` | Accelerate on macOS when the SDK provides `SparseFactorizationLU`, else UMFPACK when `pkg-config` finds it, else `none`. `feature_profile=minimal` selects `none`. | The configure summary prints the result as `SparseDirect backend`. |
| `accelerate` | Apple Accelerate sparse LU, METIS ordering | A system framework: nothing to bundle. Needs a macOS 15.5 SDK at build time and macOS 15.5 at run time (checked at run time; older systems report it unavailable). Its headers are not valid C++26, so only `inverse/direct_lu_accelerate.cpp` includes them, compiled as C++20. |
| `umfpack` | SuiteSparse UMFPACK (`dependency('UMFPACK')`), METIS ordering when the SuiteSparse build has it, otherwise its default | Portable. Linux: `libsuitesparse-dev`; macOS: `brew install suite-sparse`. A wheel must bundle `libumfpack` and its SuiteSparse and BLAS dependencies (auditwheel/delocate do this). |
| `none` | — | The strategy still exists; preparing it fails with `PRECONDITIONER_INCOMPATIBLE_STRATEGY`. `--potential-inverse auto` then uses AMG. |

An explicit `accelerate` or `umfpack` that cannot be satisfied stops the configure step instead of falling
back. At run time, `serif::solver::preconditioner::backend::mfem::sparse_direct_backend()` returns the
compiled backend name and `sparse_direct_available()` says whether it works on this machine.

## Wheels

Pass the backend per platform through meson-python, for example with cibuildwheel:

```toml
[tool.cibuildwheel.macos]
config-settings = { "setup-args" = "-Dsparse_direct=accelerate" }

[tool.cibuildwheel.linux]
before-all = "yum install -y suitesparse-devel || apt-get install -y libsuitesparse-dev"
config-settings = { "setup-args" = "-Dsparse_direct=umfpack" }
```

or `pip wheel . -Csetup-args=-Dsparse_direct=umfpack`. Setting it explicitly keeps a wheel from silently
picking up whatever the build machine happens to have; use `none` for a minimal wheel without extra
native dependencies.

## Reuse across Newton steps and batched solves

- **Symbolic reuse.** Each sparse-direct block keeps the symbolic analysis (fill-reducing ordering) of the last
  matrix it factored (`SparseDirectCache`, owned by the preconditioner's inverse computation). While the sparsity
  pattern is unchanged, which is every step after the first in the stellar solves measured, only the numeric
  factorization is repeated: 0.033 instead of 0.258 s with Accelerate, 0.17 instead of 0.22 s with UMFPACK. The
  result is bitwise the same as a fresh factorization for both backends (checked on potential blocks from two
  different Newton states), so it never changes a solve.
- **Batched solves.** The bordered preconditioner probes its border columns in chunks of 32 and passes each chunk
  to the field inverse as one `ArrayMult`; the sparse-direct inverse then solves the chunk in one Accelerate call,
  reading the factor once (32 columns: 0.08 s for 171 columns instead of 0.26 s). These columns agree with single
  solves to about 4e-15 relative, not bit for bit, so preconditioned Newton iterates change at that level.
  UMFPACK has no multiple right-hand-side solve and solves the chunk column by column (bitwise as before).

## Measured on the n=2, J=0.05 benchmark (16,189-row potential block, M4 Max)

| Backend | factor | one solve | factorization footprint |
|---|---:|---:|---:|
| Accelerate, METIS | 0.25 s | 1.5 ms | +271 MB |
| UMFPACK, METIS | 0.22 s | 10 ms | +348 MB |

Accelerate's triangular solves are multithreaded, UMFPACK's are not, so with UMFPACK each Newton step's 83
border solves cost more: the benchmark solve takes 21.4 s instead of 18.4 s, with the same 7 Newton steps and
the same converged star. On macOS the UMFPACK build's peak footprint also rose by about 200 MB per Newton step
(2.1 GiB after 7 steps, no leaks reported by `leaks`): freed factorization memory that the macOS allocator
did not return. Measure again on Linux before relying on that number.