pounce_linsol/sparse_sym_iface.rs
1//! Low-level sparse-symmetric backend interface — port of
2//! `IpSparseSymLinearSolverInterface.hpp`.
3//!
4//! Concrete implementors:
5//! * `pounce_hsl::Ma57SolverInterface` (v1.0).
6//! * Future: MUMPS, FERAL.
7
8use crate::status::ESymSolverStatus;
9use pounce_common::types::{Index, Number};
10
11/// Snapshot of the most recent LDLᵀ factor's sparsity pattern (and
12/// optionally values) plus the fill-reducing permutation. Backends
13/// produce this on demand from [`SparseSymLinearSolverInterface::factor_pattern`]
14/// — it is purely diagnostic and is not part of the solve / refine
15/// hot path.
16///
17/// All `irn` / `jcn` indices are **1-based** in *permuted* coordinates
18/// (i.e. they reference the matrix `Pᵀ K P` that the backend actually
19/// factored, not the original-variable ordering). The `perm` array
20/// closes the loop: `perm[k] = original_row` for the k-th permuted
21/// row, so a consumer can render the L pattern in either coordinate
22/// system. `perm` is **0-based** to keep the array directly indexable.
23///
24/// Only the **strict lower triangle** of L is populated — the unit
25/// diagonal is implicit (`L_ii = 1`).
26#[derive(Debug, Clone)]
27pub struct FactorPattern {
28 /// Matrix dimension (rows = cols).
29 pub n: usize,
30 /// Fill-reducing permutation, 0-based, length `n`. `perm[k]` is
31 /// the original-variable row that landed at permuted-row `k`.
32 pub perm: Vec<usize>,
33 /// Row indices of L's strict-lower nonzeros, 1-based, permuted
34 /// coordinates.
35 pub l_irn: Vec<Index>,
36 /// Column indices of L's strict-lower nonzeros, 1-based, permuted
37 /// coordinates. Same length as `l_irn`.
38 pub l_jcn: Vec<Index>,
39 /// Optional numerical values aligned with `l_irn` / `l_jcn`. `None`
40 /// when only the pattern was requested.
41 pub l_vals: Option<Vec<Number>>,
42}
43
44/// Sparse matrix format that a backend wants its triplet/CSR data in.
45/// Mirrors `SparseSymLinearSolverInterface::EMatrixFormat`.
46#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
47pub enum EMatrixFormat {
48 /// Triplet (COO) of the lower triangle, 1-based indices
49 /// (MA27 / MA57 / MUMPS convention).
50 TripletFormat,
51 /// CSR of the upper triangle, 0-based indices.
52 CsrFormat0Offset,
53 /// CSR of the upper triangle, 1-based indices.
54 CsrFormat1Offset,
55 /// Full CSR (lower + upper), 0-based indices.
56 CsrFullFormat0Offset,
57 /// Full CSR (lower + upper), 1-based indices.
58 CsrFullFormat1Offset,
59}
60
61/// Backend-side trait. The lifecycle mirrors upstream's narrative
62/// comment in `IpSparseSymLinearSolverInterface.hpp`:
63///
64/// 1. caller asks [`Self::matrix_format`].
65/// 2. caller calls [`Self::initialize_structure`] once with `(ia, ja)`.
66/// 3. caller takes the values pointer from
67/// [`Self::values_array_mut`], fills it.
68/// 4. caller calls [`Self::multi_solve`] with `new_matrix=true` for
69/// each new value pattern.
70/// 5. caller may query [`Self::number_of_neg_evals`] /
71/// [`Self::increase_quality`] between solves.
72///
73/// `new_matrix=false` requests a back-substitution against the
74/// existing factorization.
75pub trait SparseSymLinearSolverInterface {
76 /// Initialize backend internal structures for a matrix of given
77 /// dimension and pattern.
78 fn initialize_structure(
79 &mut self,
80 dim: Index,
81 nonzeros: Index,
82 ia: &[Index],
83 ja: &[Index],
84 ) -> ESymSolverStatus;
85
86 /// Slice into which the caller writes the matrix nonzeros (in the
87 /// same order as `ja` from [`Self::initialize_structure`]).
88 fn values_array_mut(&mut self) -> &mut [Number];
89
90 /// Factor (if `new_matrix`) and back-substitute against `nrhs`
91 /// right-hand sides packed in `rhs_vals` (length `nrhs * dim`).
92 /// Solutions overwrite `rhs_vals`.
93 #[allow(clippy::too_many_arguments)]
94 fn multi_solve(
95 &mut self,
96 new_matrix: bool,
97 ia: &[Index],
98 ja: &[Index],
99 nrhs: Index,
100 rhs_vals: &mut [Number],
101 check_neg_evals: bool,
102 number_of_neg_evals: Index,
103 ) -> ESymSolverStatus;
104
105 /// Number of negative eigenvalues found in the most recent
106 /// factorization. Caller must check [`Self::provides_inertia`]
107 /// first.
108 fn number_of_neg_evals(&self) -> Index;
109
110 /// Ask the backend to use a more accurate (but slower) pivot
111 /// strategy on the next solve. Returns `false` if the maximum
112 /// quality is already reached.
113 fn increase_quality(&mut self) -> bool;
114
115 /// Whether this backend reports the number of negative
116 /// eigenvalues post-factor.
117 fn provides_inertia(&self) -> bool;
118
119 /// Whether a blocked `multi_solve` of `nrhs` columns returns
120 /// **bit-identical** results to `nrhs` separate `nrhs = 1` calls
121 /// against the same factor.
122 ///
123 /// Backends that block the triangular substitution across columns
124 /// reassociate the floating-point sums, so their batched answer is
125 /// tolerance-equal but not bit-equal to the per-column one. That is
126 /// fine for a caller that only wants *a* solution, and not fine for
127 /// a caller batching purely to save time inside an iteration whose
128 /// trajectory must not move: on a nonconvex problem the perturbation
129 /// can select a different local optimum. `pooling_rt2stp` under MA57
130 /// does exactly that (gh#729), landing on an objective 25% worse
131 /// while still reporting `Optimal Solution Found`.
132 ///
133 /// Defaults to `false` — the conservative answer, so a new backend
134 /// has to opt in deliberately rather than inherit a trajectory
135 /// change by omission. This gates only opportunistic batching;
136 /// callers that batch for their own reasons (`pounce-sensitivity`'s
137 /// `jacrev` backward, where each cotangent is an independent
138 /// question) do not consult it.
139 /// The answer is allowed to depend on `nrhs`: a backend may run a
140 /// bit-identical rank-1 cascade for narrow blocks and switch to a
141 /// reassociating BLAS-3 panel kernel once the block is wide enough
142 /// to pay for it. feral does exactly that.
143 fn multi_solve_matches_single_solve(&self, _nrhs: usize) -> bool {
144 false
145 }
146
147 /// Required matrix layout. Caller marshals data into this format.
148 fn matrix_format(&self) -> EMatrixFormat;
149
150 /// Downcast seam: the concrete backend behind the trait object, for
151 /// callers that need to inspect how it was configured.
152 ///
153 /// Defaults to `None`, so a backend is invisible here until it opts
154 /// in; nothing in the solve path consults this.
155 ///
156 /// It exists because gh#825 had no observable symptom. A factory
157 /// built its MA57 backend with `Ma57SolverInterface::new()` —
158 /// hard-coded defaults — so all nine `ma57_*` options were accepted
159 /// and discarded, and every arm of every solve came out identical to
160 /// all seventeen digits. Configuration reaching a backend is not
161 /// checkable through the rest of this trait: two differently
162 /// configured factorizations are both just "a solution". This gives
163 /// a test somewhere to stand. See
164 /// `pounce-algorithm/tests/ma57_options_reach_the_backend.rs`.
165 fn as_any(&self) -> Option<&dyn std::any::Any> {
166 None
167 }
168
169 /// Whether [`Self::determine_dependent_rows`] is supported.
170 fn provides_degeneracy_detection(&self) -> bool {
171 false
172 }
173
174 /// Find the linearly dependent rows of a constraint Jacobian `J`
175 /// (the Ipopt-style degeneracy probe). `J` is `n_rows × n_cols`,
176 /// supplied as a **1-based triplet** `(irn, jcn, vals)`; on
177 /// success `c_deps` is filled with the **0-based** indices of a
178 /// set of rows whose removal leaves `J` full row rank (each
179 /// dropped row is a linear combination of the retained ones).
180 ///
181 /// Callers must check [`Self::provides_degeneracy_detection`]
182 /// first. The default returns `FatalError`, matching upstream's
183 /// "not supported" default; backends that implement this set
184 /// `provides_degeneracy_detection() -> true`.
185 fn determine_dependent_rows(
186 &mut self,
187 _n_rows: Index,
188 _n_cols: Index,
189 _irn: &[Index],
190 _jcn: &[Index],
191 _vals: &[Number],
192 _c_deps: &mut Vec<Index>,
193 ) -> ESymSolverStatus {
194 ESymSolverStatus::FatalError
195 }
196
197 /// Snapshot of the most recent factor's L pattern and permutation.
198 /// Backends that expose their factor data structures (e.g. feral)
199 /// return `Some(_)`; backends that don't (e.g. MA57, which keeps
200 /// its factors inside opaque Fortran work arrays) return `None`.
201 /// Diagnostic-only — consumed by the `--dump kkt:*+L` path. Set
202 /// `want_values=true` to populate [`FactorPattern::l_vals`].
203 fn factor_pattern(&self, _want_values: bool) -> Option<FactorPattern> {
204 None
205 }
206}