1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
//! Integration test: feed a real triplet matrix through
//! [`pounce_linsol::TSymLinearSolver`] backed by MA57 (via
//! `pounce-hsl`'s `Ma57SolverInterface`), confirm the solve hits
//! machine-precision. Gated on the `ma57` cargo feature; without
//! HSL linked in, MA57 is not available.
//!
//! Lives in `pounce-algorithm` because the algorithm crate is the
//! consumer that wires Phase-4 (linsol + hsl) into Phase-6 (KKT
//! solvers). The wrapper itself is canonical at
//! `pounce_linsol::TSymLinearSolver`; the duplicate local copy under
//! `pounce-algorithm/src/kkt/` was removed when `StdAugSystemSolver`
//! was rewired in Phase 6.
#![cfg(feature = "ma57")]
use pounce_common::types::{Index, Number};
use pounce_hsl::Ma57SolverInterface;
use pounce_linsol::{ESymSolverStatus, TSymLinearSolver};
/// Solve `A x = b` for the symmetric indefinite 3×3:
///
/// ```text
/// A = [ 4 -1 0 ] b = [ 4 ] x = ...
/// [-1 4 -1 ] [ 4 ]
/// [ 0 -1 4 ] [ 6 ]
/// ```
///
/// Lower-triangular triplet (1-based for MA57):
/// (1,1)=4, (2,1)=-1, (2,2)=4, (3,2)=-1, (3,3)=4
#[test]
fn ma57_solves_3x3_spd_via_t_sym_wrapper() {
let irn: [Index; 5] = [1, 2, 2, 3, 3];
let jcn: [Index; 5] = [1, 1, 2, 2, 3];
let vals: [Number; 5] = [4.0, -1.0, 4.0, -1.0, 4.0];
let mut tsym = TSymLinearSolver::new(Box::new(Ma57SolverInterface::new()), None, false);
assert_eq!(
tsym.initialize_structure(3, &irn, &jcn),
ESymSolverStatus::Success
);
// Hand-solved x:
// x3 = 110/56 = 55/28
// x2 = 4 x3 - 6 = 13/7
// x1 = (4 + x2)/4 = 41/28
let x1: Number = 41.0 / 28.0;
let x2: Number = 13.0 / 7.0;
let x3: Number = 55.0 / 28.0;
let mut rhs: [Number; 3] = [4.0, 4.0, 6.0];
let status = tsym.multi_solve(&vals, true, 1, &mut rhs, true, 0);
assert_eq!(status, ESymSolverStatus::Success);
assert!(
(rhs[0] - x1).abs() < 1e-12,
"x1: got {} want {}",
rhs[0],
x1
);
assert!(
(rhs[1] - x2).abs() < 1e-12,
"x2: got {} want {}",
rhs[1],
x2
);
assert!(
(rhs[2] - x3).abs() < 1e-12,
"x3: got {} want {}",
rhs[2],
x3
);
use pounce_linsol::SymLinearSolver;
assert_eq!(
<TSymLinearSolver as SymLinearSolver>::number_of_neg_evals(&tsym),
0
);
}
/// Indefinite 2×2: `[[1, 2], [2, 1]]` has eigenvalues `3, -1`, i.e. one
/// negative.
#[test]
fn ma57_reports_correct_inertia_on_indefinite_matrix() {
let irn: [Index; 3] = [1, 2, 2];
let jcn: [Index; 3] = [1, 1, 2];
let vals: [Number; 3] = [1.0, 2.0, 1.0];
let mut tsym = TSymLinearSolver::new(Box::new(Ma57SolverInterface::new()), None, false);
assert_eq!(
tsym.initialize_structure(2, &irn, &jcn),
ESymSolverStatus::Success
);
let mut rhs: [Number; 2] = [3.0, 3.0];
let status = tsym.multi_solve(&vals, true, 1, &mut rhs, false, 0);
assert_eq!(status, ESymSolverStatus::Success);
use pounce_linsol::SymLinearSolver;
assert_eq!(
<TSymLinearSolver as SymLinearSolver>::number_of_neg_evals(&tsym),
1
);
// [[1,2],[2,1]] · [1,1]^T = [3,3] → x = [1,1].
assert!((rhs[0] - 1.0_f64).abs() < 1e-12);
assert!((rhs[1] - 1.0_f64).abs() < 1e-12);
}