dune_core/linalg/
mod.rs

1//! Linear solver abstraction (D-0005). Devices never see the backend;
2//! nothing outside `linalg::backend` may call `faer` or
3//! `xingo-sparse` directly.
4//!
5//! The split mirrors the 3f5's SMP layer: `analyze` fixes the
6//! sparsity pattern (`SMPpreOrder`), `factor` refactors with new
7//! values and the diagonal gmin of `CKTop`'s gmin stepping
8//! (`SMPluFac` + `LoadGmin`), `solve_in_place` back-substitutes
9//! (`SMPsolve`).
10
11pub mod backend;
12
13/// Matrix sparsity pattern: `n`×`n`, ground already removed
14/// (0-based). `entries` are unique positions; values passed to
15/// [`LinearSolver::factor`] are parallel to it.
16#[derive(Clone, Debug)]
17pub struct Pattern {
18    pub n: usize,
19    pub entries: Vec<(usize, usize)>,
20}
21
22#[derive(Debug, thiserror::Error)]
23pub enum LinalgError {
24    /// Structural or numerical singularity; `row` is the 0-based
25    /// matrix row (caller maps it back to a node name).
26    #[error("singular matrix at row {row}")]
27    Singular { row: usize },
28}
29
30pub trait LinearSolver {
31    /// Fix the pattern for subsequent factorizations.
32    fn analyze(&mut self, pattern: &Pattern);
33
34    /// Factor with `values` (parallel to the pattern's entries).
35    /// `diag_gmin` is added to every diagonal position present in
36    /// the pattern, exactly like the 3f5's `LoadGmin`.
37    fn factor(&mut self, values: &[f64], diag_gmin: f64) -> Result<(), LinalgError>;
38
39    /// Solve `A x = rhs`, overwriting `rhs` with `x`.
40    fn solve_in_place(&mut self, rhs: &mut [f64]) -> Result<(), LinalgError>;
41
42    /// Factor the complex matrix `G + jB` over the same pattern
43    /// (`SMPcLUfac`); `re`/`im` are parallel to the pattern entries.
44    fn factor_complex(&mut self, re: &[f64], im: &[f64]) -> Result<(), LinalgError>;
45
46    /// Solve the complex system in place (`SMPcSolve`).
47    fn solve_complex_in_place(
48        &mut self,
49        rhs_re: &mut [f64],
50        rhs_im: &mut [f64],
51    ) -> Result<(), LinalgError>;
52
53    /// Solve the **transposed** complex system in place
54    /// (`SMPcaSolve`, the noise adjoint). Requires a prior
55    /// [`LinearSolver::factor_complex`].
56    fn solve_complex_transposed_in_place(
57        &mut self,
58        rhs_re: &mut [f64],
59        rhs_im: &mut [f64],
60    ) -> Result<(), LinalgError>;
61
62    /// Determinant of the complex matrix `G + jB` as a base-2
63    /// mantissa/exponent pair (`SMPcLUfac` + `SMPcDProd`): the
64    /// mantissa is normalized so `max(|re|,|im|) ∈ (0.5, 1]` like
65    /// the 3f5's `C_NORM`. A singular matrix returns a zero
66    /// mantissa — the pole-zero search treats that as a root.
67    fn det_complex(&mut self, re: &[f64], im: &[f64]) -> (num_complex::Complex64, i32);
68}
69
70/// The backend for this build: the dense faer LU by default; the
71/// `solver-xingo` feature switches every engine to the xingo-sparse LU
72/// (D-0005; the choice of default is P-3, still open — enabling
73/// `solver-xingo` is an explicit opt-in, not a new default).
74pub fn default_solver() -> Box<dyn LinearSolver> {
75    #[cfg(feature = "solver-xingo")]
76    {
77        Box::new(backend::xingo_sparse::XingoSparse::default())
78    }
79    #[cfg(all(feature = "solver-faer", not(feature = "solver-xingo")))]
80    {
81        Box::new(backend::faer_dense::FaerDense::default())
82    }
83    #[cfg(not(any(feature = "solver-faer", feature = "solver-xingo")))]
84    {
85        compile_error!(
86            "dune-core needs a linear solver backend (enable the `solver-faer` or the `solver-xingo` feature)"
87        );
88    }
89}
dune_core/linalg/mod.rs

dune_core/linalg/
mod.rs
