Source code for femorph_solver.solvers.linear._base

"""Common base class for pluggable linear solvers."""

from __future__ import annotations

from typing import ClassVar

import numpy as np
import scipy.sparse as sp

# Re-export ``SolverUnavailableError`` from the canonical exceptions
# module so callers can ``except SolverUnavailableError`` from either
# import path; new code should prefer ``from femorph_solver.exceptions
# import SolverUnavailableError``.
from femorph_solver.exceptions import SolverUnavailableError  # noqa: F401




class LinearSolverBase:
    """Protocol every linear-solver backend implements.

    Concrete subclasses set the class attributes (name, kind, spd_only,
    install_hint), override :meth:`available`, and implement
    :meth:`_factor` and :meth:`solve`.
    """

    #: Short identifier used by ``get_linear_solver("...")``.
    name: ClassVar[str] = ""

    #: ``"direct"`` (factor + forward/back solve) or ``"iterative"``.
    kind: ClassVar[str] = "direct"

    #: ``True`` if the solver requires a symmetric positive-definite matrix.
    spd_only: ClassVar[bool] = False

    #: User-facing hint printed when the backend is unavailable.
    install_hint: ClassVar[str] = ""

    def __init__(
        self,
        A: sp.spmatrix | sp.sparray,
        *,
        assume_spd: bool = False,
        assume_hermitian: bool = False,
        mixed_precision: bool | None = False,
    ) -> None:
        """``mixed_precision``:

        * ``False`` (default) — factor and solve in double precision.
        * ``True`` — factor in single precision and refine in double
          (for backends that support it; no-op elsewhere).  Halves the
          factor footprint when the backend honours the request.
        * ``None`` — let the backend choose the default based on the
          problem (``assume_spd`` + size heuristic); the
          :class:`PardisoSolver` is currently the only backend that
          opts into MP here, and only when the input is large enough
          for the memory saving to matter.

        ``assume_hermitian``:

        * ``False`` (default) — treat ``A`` as general (real
          unsymmetric or complex unsymmetric).
        * ``True`` — caller asserts ``A`` is Hermitian (real
          symmetric or complex Hermitian).  Backends that expose a
          Hermitian fast path (Pardiso ``mtype=4`` for complex
          Hermitian; ``mtype=-2`` for real symmetric indefinite)
          may pick it up.  Real SPD callers should keep using
          ``assume_spd=True`` instead — that's a stricter assertion
          (positive-definite implies Hermitian) and unlocks the
          Cholesky path.
        """
        self._n = A.shape[0]
        self._assume_spd = assume_spd
        self._assume_hermitian = assume_hermitian
        self._mixed_precision = mixed_precision
        self._factor(A)



    @staticmethod
    def available() -> bool:
        """Return ``True`` if the solver can be constructed on this install."""
        return True


    def _factor(self, A: sp.spmatrix | sp.sparray) -> None:
        """Precompute any factorization / preconditioner.

        Subclasses override.  Called once by :meth:`__init__`.
        """
        raise NotImplementedError



    def solve(self, b: np.ndarray) -> np.ndarray:
        """Solve ``A x = b`` for a single RHS."""
        raise NotImplementedError




    def refactor(self, A: sp.spmatrix | sp.sparray) -> None:
        """Update the numeric factor for a new ``A`` with the same sparsity.

        Contract: ``A`` must have **identical sparsity pattern**
        (``indptr`` + ``indices``) to the matrix this solver was
        originally constructed on; only ``A.data`` differs.  Used by
        sweep-style callers (e.g. cyclic-symmetry harmonics, parametric
        modal sweeps) to amortise the symbolic factor across many
        numeric refactorisations.

        Default implementation re-runs :meth:`_factor` from scratch —
        backends that expose a numeric-only refactor hook (Pardiso
        phase 22, UMFPACK ``numeric()`` reusing ``symbolic()``)
        override and skip the analysis cost.

        Callers must not mutate ``A``'s pattern; backends that override
        won't notice and will produce wrong results silently.
        """
        self._factor(A)