Solver-backend authoring#

How to add a new linear or eigen solver backend. The solver family is pluggable by design: every backend is a class that satisfies a small protocol; Model.solve_static / Model.solve_modal dispatch to the registered backend without caring which one is in use.

Companion to Element-kernel authoring (the per-element ke / me contract); this page is the system-level K x = b and (K - σ M) v = 0 contract.

Where solvers live#

Two parallel hierarchies under src/femorph_solver/solvers/:

Family

Path

Backends shipped today

Linear (K x = b)

solvers/linear/

cholmod (SPD direct, default), mkl_pardiso, pardiso, mumps, superlu, umfpack, cg, gmres, pyamg

Eigen ((K - σ M) v = 0)

solvers/eigen/

arpack (default), lobpcg, primme, lapack_dense (small problems / fallback)

Each backend is a single .py module under the relevant folder. __init__.py exposes the registry and the get_linear_solver(name) / get_eigen_solver(name) lookup.

The linear-solver contract#

Every linear backend subclasses femorph_solver.solvers.linear._base.LinearSolverBase:

class LinearSolverBase:
    # Class attributes
    name: ClassVar[str]            # "cholmod", "mkl_pardiso", ...
    kind: ClassVar[str]            # "direct" or "iterative"
    spd_only: ClassVar[bool]       # True for cholmod, etc.
    install_hint: ClassVar[str]    # "pip install scikit-sparse"

    @classmethod
    def available(cls) -> bool:
        """True if the optional dependency is importable."""

    def _factor(self, A: sp.csr_matrix) -> None:
        """Factor (or pre-process) the system matrix."""

    def solve(self, b: np.ndarray, x: np.ndarray | None = None) -> np.ndarray:
        """Solve ``A x = b`` against the previously-factored ``A``."""

The factorisation / solve split lets Model.solve_static reuse a factor across multiple right-hand sides (load cases, modal- superposition harmonic, etc.).

Symmetric-positive-definite-only (spd_only=True) backends are chosen by Model for static / modal problems where K is SPD; indefinite-capable backends are needed for problems with constraints, hourglass-stabilised reduced integration, or sigma < 0 shift in modal solves.

Iterative backends additionally accept a preconditioner kwarg (M) and tolerance / max-iterations through their constructor; the dispatch passes those through.

The eigen-solver contract#

Every eigen backend subclasses femorph_solver.solvers.eigen._base.EigenSolverBase:

class EigenSolverBase:
    name: ClassVar[str]
    install_hint: ClassVar[str]

    @classmethod
    def available(cls) -> bool: ...

    def solve(
        self,
        K: sp.csr_matrix,
        M: sp.csr_matrix,
        n_modes: int,
        sigma: float = 0.0,
    ) -> tuple[np.ndarray, np.ndarray]:
        """Return (eigenvalues, eigenvectors) for ``K v = λ M v``.

        ``sigma`` is a shift; non-zero shifts factor ``K - σ M``
        via shift-invert and require an SPD-or-indefinite linear
        backend underneath.
        """

The dispatch protocol#

get_linear_solver(name) (in solvers/linear/__init__.py) and get_eigen_solver(name) (in solvers/eigen/__init__.py) are the lookup points. The selection rule when no name is passed:

  1. Walk the registered backends in priority order (the order they’re listed in the registry — preferred backends first).

  2. For each, call available(); if it imports cleanly, use it.

  3. Fall back to the always-on backend (superlu for linear, lapack_dense for eigen on small problems).

  4. Raise SolverUnavailableError only if every backend fails — this should never happen in a default install.

The user can override via Model.solve_static(linear="mkl_pardiso") or by setting an environment variable (see the user-guide solver page).

Adding a new linear backend#

Worked example: a hypothetical hypre backend (PETSc / Hypre BoomerAMG for distributed sparse direct).

Step 1 — write the module#

Drop src/femorph_solver/solvers/linear/_hypre.py:

"""Hypre BoomerAMG linear solver — distributed-memory direct.

References
----------
* Henson, V. E. and Yang, U. M., 2002.  *BoomerAMG: a parallel
  algebraic multigrid solver and preconditioner.*  Applied
  Numerical Mathematics, 41, 155–177.
"""

from __future__ import annotations
import numpy as np
import scipy.sparse as sp
from femorph_solver.solvers.linear._base import (
    LinearSolverBase, SolverUnavailableError,
)


class HypreSolver(LinearSolverBase):
    name = "hypre"
    kind = "direct"
    spd_only = False
    install_hint = "pip install pyamg-hypre  # plus a working MPI"

    @classmethod
    def available(cls) -> bool:
        try:
            import pyamg_hypre  # noqa: F401
            return True
        except ImportError:
            return False

    def _factor(self, A: sp.csr_matrix) -> None:
        import pyamg_hypre
        self._solver = pyamg_hypre.boomeramg(A)

    def solve(self, b, x=None):
        return self._solver.solve_static(b)

Conventions:

  • Top-of-file References block citing the published algorithm — provenance is non-negotiable (see Provenance inventory).

  • available() returns False on a clean ImportError, not on a runtime error. Runtime failures bubble through solve.

  • Don’t import the optional dep at module-import time. The rest of the codebase imports the solver registry eagerly; importing a missing dep would cascade.

Step 2 — register#

Add to src/femorph_solver/solvers/linear/__init__.py:

from femorph_solver.solvers.linear._hypre import HypreSolver

_BACKENDS = (
    # ... existing entries unchanged ...
    HypreSolver,
)

The order in _BACKENDS is the auto-selection priority. Insert at the right tier — before the always-on superlu fallback, after the high-priority defaults.

Step 3 — unit test#

tests/solvers/linear/test_hypre.py covers:

  • available() returns False cleanly when the dep is missing (mock the import).

  • When the dep IS available, the backend solves a small SPD reference problem (the standard 1-D Poisson tridiagonal is the convention).

  • The backend matches the always-on superlu to rtol=1e-9 on the same problem (cross-backend parity).

The eigen pattern#

Eigen backends mirror the same shape:

  • Subclass EigenSolverBase, set name / install_hint, implement available() and solve(K, M, n_modes, sigma).

  • Drop into solvers/eigen/_<name>.py.

  • Register in solvers/eigen/__init__.py.

  • Unit test against a small analytical eigenproblem (a 2-DOF spring-mass system is the convention).

Shift-invert eigen backends#

Backends that support shift-invert (arpack at sigma != 0, primme) need an inner linear solver to factor K - σ M. They take a linear_solver kwarg in their constructor and dispatch that linear solve internally; arpack._arpack.py is the reference implementation.

Common pitfalls#

  • Importing the optional dep at module-import time. import pyamg_hypre at top-of-file means a missing Hypre install breaks the registry import, which breaks Model.solve_static for every user. Always lazy-import inside available() / _factor() / solve().

  • Returning a non-ndarray from solve(). The contract is np.ndarray (1-D for linear, 2-D for eigen vectors). Some backends return a wrapper type that acts like an ndarray but breaks downstream code that does np.asarray(result, dtype=np.float64). Always np.asarray(...) before returning.

  • Mutating the input A / b. The dispatch caller may reuse the matrix across solves; in-place modifications cascade.

  • Not handling spd_only=False at the modal layer. Model.solve_modal(sigma=-1.0) factors K + |σ| M which is SPD when M is positive definite; that’s still the spd_only=True path. sigma > 0 factors K - σ M which is generally indefinite and needs a non-SPD backend.

  • Forgetting the install_hint. The hint is what the SolverUnavailableError message tells the user to do; an empty string makes the error unactionable.

Where things live#

Concern

Path

Linear-backend module

src/femorph_solver/solvers/linear/_<name>.py

Eigen-backend module

src/femorph_solver/solvers/eigen/_<name>.py

Linear registry / dispatch

src/femorph_solver/solvers/linear/__init__.py

Eigen registry / dispatch

src/femorph_solver/solvers/eigen/__init__.py

Linear base class

src/femorph_solver/solvers/linear/_base.py

Eigen base class

src/femorph_solver/solvers/eigen/_base.py

Per-backend unit tests

tests/solvers/{linear,eigen}/test_<name>.py

User-guide solver page

doc/source/user-guide/solving/ (linked from this page)