pounce_sensitivity/backsolver.rs
1//! `SensBacksolver` trait — abstract backsolver against a converged KKT factor.
2//!
3//! Mirrors upstream
4//! [`SensBacksolver.hpp`](../../../ref/Ipopt/contrib/sIPOPT/src/SensBacksolver.hpp).
5//!
6//! # What this is
7//!
8//! After pounce's IPM converges, the KKT factor `K` is stored on the
9//! algorithm side (`pounce-algorithm::kkt::PdFullSpaceSolver`).
10//! sIPOPT runs a small number of backsolves against that factor to
11//! build the sensitivity matrix `P = K⁻¹ A` (see
12//! [`crate::p_calculator::PCalculator`]). The trait surface is just
13//! "solve `K · lhs = rhs`"; concrete impls plug in either the real
14//! `AugSystemSolver` (Phase B.2) or a synthetic dense-LU backsolver
15//! (Phase B.1, this file, for unit testing the sensitivity math).
16//!
17//! Upstream `SensSimpleBacksolver` is the analogous wrapper around
18//! Ipopt's `AugSystemSolver`
19//! ([`SensSimpleBacksolver.{hpp,cpp}`](../../../ref/Ipopt/contrib/sIPOPT/src/SensSimpleBacksolver.hpp)).
20
21use pounce_common::types::Number;
22
23/// Solve `K · lhs = rhs` against the converged KKT factor. Returns
24/// `false` on failure (e.g. backend reports `Singular`).
25///
26/// Mirrors `Ipopt::SensBacksolver::Solve`
27/// ([`SensBacksolver.hpp:28-31`](../../../ref/Ipopt/contrib/sIPOPT/src/SensBacksolver.hpp)).
28/// Pounce takes flat `&[Number]` / `&mut [Number]` rather than
29/// upstream's block-structured `IteratesVector` because the
30/// sensitivity-side data is naturally flat; if the algorithm-side
31/// wrapper in Phase B.2 needs the block layout it converts before
32/// calling.
33pub trait SensBacksolver {
34 /// Size of the linear system in entries (length of `lhs` and
35 /// `rhs`). The backsolver's notion of "full state" — pounce's
36 /// IPM uses the compound `(x, s, λ_c, λ_d, z_l, z_u, v_l, v_u)`
37 /// concatenation here.
38 fn dim(&self) -> usize;
39
40 /// Solve `K · lhs = rhs`. The implementation may use `rhs` as
41 /// scratch; callers should treat it as moved-from on return.
42 /// `lhs` must have length `self.dim()`.
43 fn solve(&self, rhs: &[Number], lhs: &mut [Number]) -> bool;
44}
45
46/// Synthetic dense-LU backsolver. Used in this crate's tests to
47/// validate the sensitivity math against known-good linear-algebra
48/// answers without standing up the full pounce IPM. Phase B.2 ships
49/// a real `pounce-algorithm`-backed implementation; this stays for
50/// regression tests and as a reference for the trait contract.
51///
52/// Stores a row-major `n × n` matrix and reuses an in-place
53/// Gaussian-elimination LU factor with partial pivoting. Numerical
54/// stability is fine for the small problem sizes the unit tests
55/// exercise (n ≤ 16).
56#[derive(Debug, Clone)]
57pub struct DenseLuBacksolver {
58 n: usize,
59 /// Factored `K` in row-major order: contains `L` (unit lower)
60 /// and `U` (upper) packed in-place per Doolittle.
61 lu: Vec<Number>,
62 /// Row-permutation order: row `piv[i]` of the original `K`
63 /// ended up in row `i` of the factor.
64 piv: Vec<usize>,
65}
66
67impl DenseLuBacksolver {
68 /// Build from a row-major `n × n` matrix. Returns `Err(())` if
69 /// the matrix is exactly singular at LU time (zero pivot after
70 /// pivoting).
71 pub fn from_dense(n: usize, a_row_major: &[Number]) -> Result<Self, ()> {
72 if a_row_major.len() != n * n {
73 return Err(());
74 }
75 let mut lu = a_row_major.to_vec();
76 let mut piv: Vec<usize> = (0..n).collect();
77 for k in 0..n {
78 // Partial pivot: find row with largest |lu[r,k]| for r >= k.
79 let mut best = k;
80 let mut best_mag = lu[k * n + k].abs();
81 for r in (k + 1)..n {
82 let mag = lu[r * n + k].abs();
83 if mag > best_mag {
84 best = r;
85 best_mag = mag;
86 }
87 }
88 if best_mag == 0.0 {
89 return Err(());
90 }
91 if best != k {
92 piv.swap(k, best);
93 for j in 0..n {
94 let tmp = lu[k * n + j];
95 lu[k * n + j] = lu[best * n + j];
96 lu[best * n + j] = tmp;
97 }
98 }
99 // Eliminate below pivot.
100 let inv_p = 1.0 / lu[k * n + k];
101 for r in (k + 1)..n {
102 let m = lu[r * n + k] * inv_p;
103 lu[r * n + k] = m;
104 for j in (k + 1)..n {
105 let upd = lu[k * n + j];
106 lu[r * n + j] -= m * upd;
107 }
108 }
109 }
110 Ok(Self { n, lu, piv })
111 }
112}
113
114impl SensBacksolver for DenseLuBacksolver {
115 fn dim(&self) -> usize {
116 self.n
117 }
118
119 fn solve(&self, rhs: &[Number], lhs: &mut [Number]) -> bool {
120 if rhs.len() != self.n || lhs.len() != self.n {
121 return false;
122 }
123 // Apply row permutation first.
124 for i in 0..self.n {
125 lhs[i] = rhs[self.piv[i]];
126 }
127 // Forward solve: L · y = P · rhs (L unit lower in `lu`).
128 for i in 0..self.n {
129 let mut s = lhs[i];
130 for j in 0..i {
131 s -= self.lu[i * self.n + j] * lhs[j];
132 }
133 lhs[i] = s;
134 }
135 // Back solve: U · x = y.
136 for i in (0..self.n).rev() {
137 let mut s = lhs[i];
138 for j in (i + 1)..self.n {
139 s -= self.lu[i * self.n + j] * lhs[j];
140 }
141 let p = self.lu[i * self.n + i];
142 if p == 0.0 {
143 return false;
144 }
145 lhs[i] = s / p;
146 }
147 true
148 }
149}
150
151#[cfg(test)]
152mod tests {
153 use super::*;
154
155 /// Solve the small symmetric system `A x = b`:
156 /// 2 -1 0 | 1
157 /// -1 2 -1 | 0
158 /// 0 -1 2 | 0
159 /// Closed-form: x = (3/4, 1/2, 1/4).
160 #[test]
161 fn dense_lu_solves_3x3_symmetric() {
162 #[rustfmt::skip]
163 let a = vec![
164 2.0, -1.0, 0.0,
165 -1.0, 2.0, -1.0,
166 0.0, -1.0, 2.0,
167 ];
168 let solver = DenseLuBacksolver::from_dense(3, &a).expect("factor");
169 let b = [1.0, 0.0, 0.0];
170 let mut x = [0.0; 3];
171 assert!(solver.solve(&b, &mut x));
172 assert!((x[0] - 0.75).abs() < 1e-12, "x[0] = {}", x[0]);
173 assert!((x[1] - 0.50).abs() < 1e-12, "x[1] = {}", x[1]);
174 assert!((x[2] - 0.25).abs() < 1e-12, "x[2] = {}", x[2]);
175 }
176
177 /// Pivoting check: a matrix whose first pivot is zero (must swap
178 /// rows). Direct closed-form: A x = b for
179 /// 0 1 | 2 → x = (2, 1) after swap
180 /// 1 0 | 1
181 #[test]
182 fn dense_lu_handles_zero_first_pivot() {
183 let a = vec![0.0, 1.0, 1.0, 0.0];
184 let solver = DenseLuBacksolver::from_dense(2, &a).expect("factor");
185 let b = [2.0, 1.0];
186 let mut x = [0.0; 2];
187 assert!(solver.solve(&b, &mut x));
188 assert!((x[0] - 1.0).abs() < 1e-12, "x[0] = {}", x[0]);
189 assert!((x[1] - 2.0).abs() < 1e-12, "x[1] = {}", x[1]);
190 }
191
192 #[test]
193 fn dense_lu_rejects_singular_matrix() {
194 // Rank-deficient: rows are linearly dependent.
195 let a = vec![1.0, 2.0, 2.0, 4.0];
196 assert!(DenseLuBacksolver::from_dense(2, &a).is_err());
197 }
198
199 #[test]
200 fn solve_rejects_wrong_dim() {
201 let a = vec![1.0, 0.0, 0.0, 1.0];
202 let s = DenseLuBacksolver::from_dense(2, &a).expect("ok");
203 let b = [1.0];
204 let mut x = [0.0; 2];
205 assert!(!s.solve(&b, &mut x));
206 }
207}