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    /// Required matrix layout. Caller marshals data into this format.
120    fn matrix_format(&self) -> EMatrixFormat;
121
122    /// Whether [`Self::determine_dependent_rows`] is supported.
123    fn provides_degeneracy_detection(&self) -> bool {
124        false
125    }
126
127    /// Find the linearly dependent rows of a constraint Jacobian `J`
128    /// (the Ipopt-style degeneracy probe). `J` is `n_rows × n_cols`,
129    /// supplied as a **1-based triplet** `(irn, jcn, vals)`; on
130    /// success `c_deps` is filled with the **0-based** indices of a
131    /// set of rows whose removal leaves `J` full row rank (each
132    /// dropped row is a linear combination of the retained ones).
133    ///
134    /// Callers must check [`Self::provides_degeneracy_detection`]
135    /// first. The default returns `FatalError`, matching upstream's
136    /// "not supported" default; backends that implement this set
137    /// `provides_degeneracy_detection() -> true`.
138    fn determine_dependent_rows(
139        &mut self,
140        _n_rows: Index,
141        _n_cols: Index,
142        _irn: &[Index],
143        _jcn: &[Index],
144        _vals: &[Number],
145        _c_deps: &mut Vec<Index>,
146    ) -> ESymSolverStatus {
147        ESymSolverStatus::FatalError
148    }
149
150    /// Snapshot of the most recent factor's L pattern and permutation.
151    /// Backends that expose their factor data structures (e.g. feral)
152    /// return `Some(_)`; backends that don't (e.g. MA57, which keeps
153    /// its factors inside opaque Fortran work arrays) return `None`.
154    /// Diagnostic-only — consumed by the `--dump kkt:*+L` path. Set
155    /// `want_values=true` to populate [`FactorPattern::l_vals`].
156    fn factor_pattern(&self, _want_values: bool) -> Option<FactorPattern> {
157        None
158    }
159}
pounce_linsol/sparse_sym_iface.rs

pounce_linsol/
sparse_sym_iface.rs
