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 :doc:`kernel_authoring` (the **per-element**
``ke`` / ``me`` contract); this page is the **system-level**
``K x = b`` and ``(K - σ M) v = 0`` contract.

.. contents:: Page contents
    :local:
    :depth: 2

Where solvers live
------------------

Two parallel hierarchies under
``src/femorph_solver/solvers/``:

.. list-table::
    :header-rows: 1
    :widths: 26 26 48

    * - 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
:class:`femorph_solver.solvers.linear._base.LinearSolverBase`:

.. code-block:: python

    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
:class:`femorph_solver.solvers.eigen._base.EigenSolverBase`:

.. code-block:: python

    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:

#. Walk the registered backends in **priority order** (the order
   they're listed in the registry — preferred backends first).
#. For each, call ``available()``; if it imports cleanly, use it.
#. Fall back to the always-on backend (``superlu`` for linear,
   ``lapack_dense`` for eigen on small problems).
#. 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``:

.. code-block:: python

    """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 :doc:`provenance`).
* ``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``:

.. code-block:: python

    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
-----------------

.. list-table::
    :header-rows: 1
    :widths: 32 68

    * - 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)