pounce_linsol/scaling.rs
1//! Symmetric-matrix scaling for triplet inputs.
2//!
3//! Port of `Algorithm/LinearSolvers/IpTSymScalingMethod.hpp`. A scaling
4//! method takes the matrix in triplet form `(airn, ajcn, a)` and writes
5//! a per-row scaling factor `s[i]` to `scaling_factors`. The
6//! `TSymLinearSolver` wrapper then applies the symmetric scaling
7//! `A' = diag(s) · A · diag(s)` (and the inverse to the RHS / solution)
8//! before / after delegating to the backend.
9//!
10//! Variants registered upstream:
11//!
12//! * `none` — no scaling, default in many problem classes
13//! ([`IdentityScalingMethod`]).
14//! * `mc19` — HSL MC19 row/column scaling. Bit-equivalence-default;
15//! implemented as `pounce_hsl::Mc19TSymScalingMethod` (FFI to
16//! `libcoinhsl.dylib`'s `mc19ad_`).
17//! * `slack-based` — slack-aware scaling driven by the current
18//! barrier slacks. Implemented as
19//! `pounce_algorithm::kkt::SlackBasedTSymScalingMethod`; lives in
20//! the algorithm crate because it reads `IpoptData::curr` /
21//! `IpoptCq::curr_slack_*`, which would otherwise create a
22//! circular dependency.
23
24use pounce_common::types::{Index, Number};
25
26/// Backend-agnostic scaling method.
27///
28/// Returns `true` on success. On `false` the caller must skip scaling
29/// (mirrors upstream's `ComputeSymTScalingFactors` contract).
30pub trait TSymScalingMethod {
31 fn compute_sym_t_scaling_factors(
32 &mut self,
33 n: Index,
34 nnz: Index,
35 airn: &[Index],
36 ajcn: &[Index],
37 a: &[Number],
38 scaling_factors: &mut [Number],
39 ) -> bool;
40
41 /// Hand the method the per-iterate data it needs, ahead of the
42 /// factorization. Default no-op: matrix-only methods (Ruiz, MC19,
43 /// identity) derive everything from the triplets they are given and
44 /// ignore this.
45 ///
46 /// Exists because upstream's slack-based method is an algorithm
47 /// strategy object that reads `IpCq()` and `IpNLP()` directly, and
48 /// this crate is below the algorithm and cannot. The iterate-derived
49 /// part is computed by the caller and pushed in here; the block
50 /// layout and the constant blocks stay in the method, where the
51 /// upstream algorithm keeps them.
52 ///
53 /// `nx` is the primal dimension and `s_scale` the `s`-block factors;
54 /// the remaining blocks are 1. Called at most once per iteration,
55 /// not once per solve — the quantity is a function of the iterate,
56 /// and several augmented solves share one iterate.
57 fn set_slack_scaling(&mut self, _nx: Index, _s_scale: &[Number]) {}
58}
59
60/// `linear_system_scaling=none` — write identity scaling factors. The
61/// `TSymLinearSolver` wrapper detects this case and skips the symmetric
62/// scaling pass entirely; this implementation exists so that callers
63/// who hand a scaling method unconditionally get a sensible default.
64#[derive(Debug, Default, Clone, Copy)]
65pub struct IdentityScalingMethod;
66
67impl TSymScalingMethod for IdentityScalingMethod {
68 fn compute_sym_t_scaling_factors(
69 &mut self,
70 n: Index,
71 _nnz: Index,
72 _airn: &[Index],
73 _ajcn: &[Index],
74 _a: &[Number],
75 scaling_factors: &mut [Number],
76 ) -> bool {
77 debug_assert_eq!(scaling_factors.len(), n as usize);
78 for s in scaling_factors.iter_mut() {
79 *s = 1.0;
80 }
81 true
82 }
83}
84
85/// `linear_system_scaling=slack-based` — port of
86/// `IpSlackBasedTSymScalingMethod.cpp:ComputeSymTScalingFactors`.
87///
88/// The augmented system is ordered `[x | s | y_c | y_d]`, and upstream
89/// writes
90///
91/// ```text
92/// x block 1
93/// s block min(Pd_L · slack_s_L + Pd_U · slack_s_U, 1)
94/// y_c, y_d blocks 1
95/// ```
96///
97/// Only the `s` block depends on the iterate. It arrives through
98/// [`TSymScalingMethod::set_slack_scaling`], computed by
99/// `IpoptCq::curr_slack_based_s_scaling` — see that method for why the
100/// split lands here rather than inside this type.
101///
102/// Until the first `set_slack_scaling` this behaves as identity, so a
103/// factorization that happens before the algorithm has an iterate (the
104/// least-square multiplier estimate at initialization, for instance) is
105/// scaled the way `none` would scale it rather than by a stale or empty
106/// vector.
107#[derive(Debug, Default, Clone)]
108pub struct SlackBasedTSymScalingMethod {
109 nx: Index,
110 s_scale: Vec<Number>,
111}
112
113impl SlackBasedTSymScalingMethod {
114 pub fn new() -> Self {
115 Self::default()
116 }
117}
118
119impl TSymScalingMethod for SlackBasedTSymScalingMethod {
120 fn set_slack_scaling(&mut self, nx: Index, s_scale: &[Number]) {
121 self.nx = nx;
122 self.s_scale.clear();
123 self.s_scale.extend_from_slice(s_scale);
124 }
125
126 fn compute_sym_t_scaling_factors(
127 &mut self,
128 n: Index,
129 _nnz: Index,
130 _airn: &[Index],
131 _ajcn: &[Index],
132 _a: &[Number],
133 scaling_factors: &mut [Number],
134 ) -> bool {
135 debug_assert_eq!(scaling_factors.len(), n as usize);
136 for s in scaling_factors.iter_mut() {
137 *s = 1.0;
138 }
139 // No iterate yet, or an augmented system this method was not
140 // built for: identity is the honest answer, and it is what
141 // `none` would have produced anyway.
142 let nx = self.nx as usize;
143 let ns = self.s_scale.len();
144 if ns == 0 || nx + ns > n as usize {
145 return true;
146 }
147 scaling_factors[nx..nx + ns].copy_from_slice(&self.s_scale);
148 true
149 }
150}
151
152#[cfg(test)]
153mod tests {
154 use super::*;
155
156 #[test]
157 fn slack_based_is_identity_before_the_first_iterate() {
158 // The initialization-time least-square solve factorizes before
159 // any iterate exists. Scaling it with an empty vector would
160 // silently leave the s block at whatever `scaling_factors`
161 // happened to hold.
162 let mut m = SlackBasedTSymScalingMethod::new();
163 let mut f = vec![0.0; 6];
164 assert!(m.compute_sym_t_scaling_factors(6, 0, &[], &[], &[], &mut f));
165 assert_eq!(f, vec![1.0; 6]);
166 }
167
168 #[test]
169 fn slack_based_writes_only_the_s_block() {
170 // n = 6 laid out as [x x | s s | y y].
171 let mut m = SlackBasedTSymScalingMethod::new();
172 m.set_slack_scaling(2, &[0.25, 0.5]);
173 let mut f = vec![0.0; 6];
174 assert!(m.compute_sym_t_scaling_factors(6, 0, &[], &[], &[], &mut f));
175 assert_eq!(f, vec![1.0, 1.0, 0.25, 0.5, 1.0, 1.0]);
176 }
177
178 #[test]
179 fn slack_based_declines_a_system_it_does_not_fit() {
180 // A shorter system than the stored blocks describe means this
181 // method is being asked about something else (the restoration
182 // sub-IPM's augmented system, say). Writing the s block anyway
183 // would scale unrelated rows.
184 let mut m = SlackBasedTSymScalingMethod::new();
185 m.set_slack_scaling(4, &[0.25, 0.5]);
186 let mut f = vec![0.0; 5];
187 assert!(m.compute_sym_t_scaling_factors(5, 0, &[], &[], &[], &mut f));
188 assert_eq!(f, vec![1.0; 5], "must fall back to identity, not misplace");
189 }
190
191 #[test]
192 fn identity_writes_unit_factors() {
193 let mut method = IdentityScalingMethod;
194 let irn = [1, 2, 2];
195 let jcn = [1, 1, 2];
196 let vals = [2.0, 1.0, 3.0];
197 let mut s = vec![0.0; 2];
198 assert!(method.compute_sym_t_scaling_factors(2, 3, &irn, &jcn, &vals, &mut s));
199 assert_eq!(s, &[1.0, 1.0]);
200 }
201}