Optional parallel sparse solver for FD (PARDISO/MUMPS via scikit-umfpack or pyMUMPS)

[awickert]

Summary

The FD solver currently uses UMFpack via scipy.sparse.linalg.spsolve(..., use_umfpack=True), which is single-threaded. For large grids the sparse LU factorization dominates runtime (e.g., ~5 s at 400×400). Parallel sparse direct solvers — PARDISO (Intel MKL) or MUMPS — can factorize in near-linear time with core count.

Options

PARDISO (via pypardiso):

• Wraps Intel MKL PARDISO; drop-in replacement (from pypardiso import spsolve)
• Multithreaded factorization and solve; typically 4–8× faster than UMFpack on 8+ cores
• Intel hardware only: optimized for Intel silicon and may underperform or behave unpredictably on AMD CPUs due to MKL's CPU dispatch. Not recommended for AMD-based machines.
• Requires Intel MKL at runtime (available in most conda environments; less reliable on pure pip installs)

CHOLMOD (via scikit-sparse, SuiteSparse):

• Hardware-agnostic — works correctly on AMD, ARM, and any other CPU
• Cholesky factorization: faster than LU but requires the coefficient matrix to be symmetric positive definite (SPD)
• The gFlex biharmonic FD matrix (∇²(D∇²·)) should be SPD for most standard BCs, but this needs verification across all BC combinations before relying on it
• conda install -c conda-forge scikit-sparse

MUMPS (via python-mumps):

• Open-source, MPI-parallel, hardware-agnostic
• Handles non-symmetric matrices — no SPD requirement
• More complex installation; better suited to HPC/cluster scenarios than single-node use
• conda install -c conda-forge python-mumps

Suggested approach

• Keep UMFpack as the default (no new required dependencies)
• Check for an available solver at runtime and use it, or expose flex.sparse_solver = "pardiso" / "cholmod" / "mumps" as an opt-in
• Document in configuration reference

Expected benefit

Near-linear speedup in core count for the factorization step on large grids. Most impactful for 200×200 and above where factorization takes >0.5 s.

Cost

• Optional dependency only — no change for users without it
• CHOLMOD requires verifying SPD property across all gFlex BC combinations before use
• pypardiso is Intel-specific; MUMPS adds MPI complexity
• All options need benchmarking and correctness testing across BC combinations before shipping

[awickert]

A few implementation notes to flesh out the work spec:

Cache integration complexity

pypardiso does not expose the same callable interface as scipy.sparse.linalg.factorized(). Instead it uses a PyPardisoSolver object with separate .factorize(A) and .solve(b) calls. The cache_factorization = True / "no_check" paths in fd_solve() currently rely on the factorized() callable; those paths would need separate handling for pypardiso (store the solver object rather than a callable).

Suggested implementation order

1. Verify SPD — run a quick test across all BC combinations (clamped, free, pinned, mirror, periodic, no_outside_loads) to confirm the coefficient matrix is symmetric positive-definite. This is a prerequisite for CHOLMOD and takes an hour.
2. CHOLMOD first (via scikit-sparse) — hardware-agnostic, no Intel dependency, faster than UMFpack for SPD matrices. Implement if SPD is confirmed.
3. pypardiso second — genuine speedup on Intel hardware; document the AMD caveat clearly.

UMFpack should remain the default throughout; all alternatives are opt-in via flex.sparse_solver.

🤖 Generated with Claude Code