pounce_sensitivity/algorithm_backsolver.rs
1//! `PdSensBacksolver` — `SensBacksolver` adapter over the converged
2//! `PdFullSpaceSolver` from `pounce-algorithm`.
3//!
4//! This is the Phase B.2 piece tracked in
5//! [pounce#16](https://github.com/jkitchin/pounce/issues/16): it lets
6//! `pounce-sensitivity` drive backsolves against the real converged
7//! KKT factor, replacing the synthetic [`crate::DenseLuBacksolver`]
8//! used by Phase B.1 unit tests.
9//!
10//! # Use
11//!
12//! 1. Register an `on_converged` callback on `IpoptApplication` via
13//! [`pounce_algorithm::application::IpoptApplication::set_on_converged`].
14//! 2. Inside the callback, build a `PdSensBacksolver` from the four
15//! handles passed in (`data`, `cq`, `nlp`, `&mut pd_solver`).
16//! 3. Hand it to [`crate::SensApplication`] / a `SensStepCalc` /
17//! [`crate::compute_reduced_hessian`] like any other
18//! [`SensBacksolver`].
19//!
20//! Upstream `SensSimpleBacksolver`
21//! ([`ref/Ipopt/contrib/sIPOPT/src/SensSimpleBacksolver.cpp`](../../../ref/Ipopt/contrib/sIPOPT/src/SensSimpleBacksolver.cpp))
22//! is the analogous wrapper around `IpoptCalculatedQuantities` +
23//! `PDSystemSolver` upstream.
24//!
25//! # Flat-slice ↔ `IteratesVector` mapping
26//!
27//! The full primal-dual state of pounce's IPM is the eight-block
28//! compound `(x, s, λ_c, λ_d, z_l, z_u, v_l, v_u)` (see
29//! [`pounce_algorithm::iterates_vector::IteratesVector`]). This
30//! adapter packs / unpacks the flat slices that
31//! [`crate::SensBacksolver`] takes as the concatenation
32//! `x || s || λ_c || λ_d || z_l || z_u || v_l || v_u`, mirroring
33//! upstream's `CompoundVector` layout (`IpCompoundVector.hpp`).
34//!
35//! # Reference
36//!
37//! Pirnay, H.; López-Negrete, R.; Biegler, L. T. (2012). *Optimal
38//! sensitivity based on IPOPT*. Mathematical Programming Computation,
39//! **4**(4), 307–331. DOI:
40//! [10.1007/s12532-012-0043-2](https://doi.org/10.1007/s12532-012-0043-2).
41//! Verified via Crossref on 2026-05-13.
42
43use std::cell::RefCell;
44use std::rc::Rc;
45
46use pounce_algorithm::ipopt_cq::IpoptCqHandle;
47use pounce_algorithm::ipopt_data::IpoptDataHandle;
48use pounce_algorithm::iterates_vector::{IteratesVector, IteratesVectorMut};
49use pounce_algorithm::kkt::pd_full_space_solver::PdFullSpaceSolver;
50use pounce_common::types::Number;
51use pounce_linalg::dense_vector::DenseVector;
52use pounce_nlp::ipopt_nlp::IpoptNlp;
53
54use crate::backsolver::SensBacksolver;
55
56/// Adapter from `PdFullSpaceSolver` to [`SensBacksolver`]. Holds
57/// owning clones of the four pieces of the algorithm's converged
58/// state, plus the 8-block iterate template used to allocate fresh
59/// RHS / LHS vectors.
60///
61/// The PD solver lives behind an `Rc<RefCell<…>>` because
62/// [`SensBacksolver::solve`] is `&self` but the upstream signature
63/// for `PdFullSpaceSolver::solve` is `&mut self` (it caches the
64/// last-solve dependency tags and the augsys-improved flag). The
65/// `RefCell` is single-thread-only, single-borrow, exactly matching
66/// the call pattern from `pounce-sensitivity`'s pipeline.
67///
68/// Owning (rather than borrowing) the four handles is what lets a
69/// `PdSensBacksolver` outlive the `on_converged` callback frame —
70/// required by the public `Solver` session API in `pounce-algorithm`,
71/// which retains the backsolver for repeated `parametric_step` /
72/// `kkt_solve` / `compute_reduced_hessian` calls after the IPM has
73/// returned. The data, cq, and nlp handles are already
74/// `Rc<RefCell<…>>` cheap-clone handles upstream, so this carries no
75/// allocation overhead.
76#[derive(Clone)]
77pub struct PdSensBacksolver {
78 /// Shared, interior-mutable handle to the converged PD solver.
79 /// Cloned from `PdSearchDirCalc::pd_solver_rc()` at construction.
80 pd: Rc<RefCell<PdFullSpaceSolver>>,
81 data: IpoptDataHandle,
82 cq: IpoptCqHandle,
83 nlp: Rc<RefCell<dyn IpoptNlp>>,
84 /// Block dimensions in `(x, s, y_c, y_d, z_l, z_u, v_l, v_u)` order.
85 dims: [usize; 8],
86 /// 8-block prototype used to mint fresh vectors with the same
87 /// `VectorSpace`s as the converged iterate; cloned from
88 /// `data.borrow().curr`.
89 template: IteratesVector,
90}
91
92impl PdSensBacksolver {
93 /// Construct from the four handles handed in by the `on_converged`
94 /// callback. Returns `Err(())` if `data` has no `curr` (i.e. the
95 /// algorithm never reached an iterate — should not happen on
96 /// `SolveSucceeded`).
97 pub fn new(
98 data: &IpoptDataHandle,
99 cq: &IpoptCqHandle,
100 nlp: &Rc<RefCell<dyn IpoptNlp>>,
101 pd: Rc<RefCell<PdFullSpaceSolver>>,
102 ) -> Result<Self, ()> {
103 let curr = data.borrow().curr.clone().ok_or(())?;
104 let dims = [
105 curr.x.dim() as usize,
106 curr.s.dim() as usize,
107 curr.y_c.dim() as usize,
108 curr.y_d.dim() as usize,
109 curr.z_l.dim() as usize,
110 curr.z_u.dim() as usize,
111 curr.v_l.dim() as usize,
112 curr.v_u.dim() as usize,
113 ];
114 Ok(Self {
115 pd,
116 data: Rc::clone(data),
117 cq: Rc::clone(cq),
118 nlp: Rc::clone(nlp),
119 dims,
120 template: curr,
121 })
122 }
123
124 /// Block dimensions of the compound KKT vector at convergence, in
125 /// `(x, s, y_c, y_d, z_l, z_u, v_l, v_u)` order. Sum equals
126 /// [`SensBacksolver::dim`]. Useful when a caller needs to compute
127 /// the flat offset of a non-x block (e.g. `n_x + n_s` for the
128 /// start of the equality-multiplier `y_c` block).
129 pub fn block_dims(&self) -> [usize; 8] {
130 self.dims
131 }
132
133 /// Cumulative block offsets: `offset(i)` is the start index of
134 /// block `i` in the flat slice.
135 fn offsets(&self) -> [usize; 9] {
136 let mut o = [0usize; 9];
137 for i in 0..8 {
138 o[i + 1] = o[i] + self.dims[i];
139 }
140 o
141 }
142
143 /// Pack a flat slice into a freshly-allocated `IteratesVectorMut`
144 /// shaped like the converged iterate.
145 fn pack(&self, flat: &[Number]) -> Result<IteratesVectorMut, ()> {
146 let mut out = self.template.make_new_zeroed();
147 let off = self.offsets();
148 let blocks: [&mut Box<dyn pounce_linalg::vector::Vector>; 8] = [
149 &mut out.x,
150 &mut out.s,
151 &mut out.y_c,
152 &mut out.y_d,
153 &mut out.z_l,
154 &mut out.z_u,
155 &mut out.v_l,
156 &mut out.v_u,
157 ];
158 for (i, blk) in blocks.into_iter().enumerate() {
159 let slice = &flat[off[i]..off[i + 1]];
160 let dv = blk.as_any_mut().downcast_mut::<DenseVector>().ok_or(())?;
161 dv.set_values(slice);
162 }
163 Ok(out)
164 }
165
166 /// Read an `IteratesVectorMut` into a flat slice. Uses
167 /// [`DenseVector::expanded_values`] rather than `values()` so
168 /// blocks that the IPM left in homogeneous-scalar form (typical
169 /// for empty z_l/z_u/v_l/v_u when the TNLP has no bounds) are
170 /// materialized rather than panicking.
171 fn unpack(&self, iv: &IteratesVectorMut, out: &mut [Number]) -> Result<(), ()> {
172 let off = self.offsets();
173 let blocks: [&Box<dyn pounce_linalg::vector::Vector>; 8] = [
174 &iv.x, &iv.s, &iv.y_c, &iv.y_d, &iv.z_l, &iv.z_u, &iv.v_l, &iv.v_u,
175 ];
176 for (i, blk) in blocks.into_iter().enumerate() {
177 let dst = &mut out[off[i]..off[i + 1]];
178 if dst.is_empty() {
179 continue;
180 }
181 let dv = (**blk).as_any().downcast_ref::<DenseVector>().ok_or(())?;
182 let ev = dv.expanded_values();
183 dst.copy_from_slice(&ev);
184 }
185 Ok(())
186 }
187}
188
189impl SensBacksolver for PdSensBacksolver {
190 fn dim(&self) -> usize {
191 self.dims.iter().sum()
192 }
193
194 fn solve(&self, rhs: &[Number], lhs: &mut [Number]) -> bool {
195 let total = self.dim();
196 if rhs.len() != total || lhs.len() != total {
197 return false;
198 }
199 // Pack rhs into block form.
200 let rhs_mut = match self.pack(rhs) {
201 Ok(v) => v,
202 Err(()) => return false,
203 };
204 let rhs_iv = rhs_mut.freeze();
205 // Fresh result slot, zeroed.
206 let mut res_iv = self.template.make_new_zeroed();
207
208 // K · lhs = rhs ⇒ solve(α=1, β=0, rhs, res) writes
209 // res = K⁻¹ · rhs.
210 let ok = {
211 let mut pd_ref = self.pd.borrow_mut();
212 pd_ref.solve(
213 &self.data,
214 &self.cq,
215 &self.nlp,
216 1.0,
217 0.0,
218 &rhs_iv,
219 &mut res_iv,
220 false,
221 false,
222 )
223 };
224 if !ok {
225 return false;
226 }
227 self.unpack(&res_iv, lhs).is_ok()
228 }
229}