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
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
use crate::error::FeralError;
/// Symmetric matrix stored as full n×n column-major. Only the lower triangle
/// is meaningful; the strict upper triangle is ignored on input.
/// Entry (i, j) is at index j*n + i. Size: n*n f64 values.
pub struct SymmetricMatrix {
pub n: usize,
pub data: Vec<f64>,
}
impl SymmetricMatrix {
/// Create a new n×n symmetric matrix initialized to zero.
pub fn zeros(n: usize) -> Self {
Self {
n,
data: vec![0.0; n * n],
}
}
/// Reuse a pooled buffer to construct an n×n `SymmetricMatrix`
/// with the lower triangle zeroed. The strict upper triangle is
/// left with whatever stale contents the buffer held — the
/// multifrontal kernel (`factor_frontal_blocked_in_place`) and
/// every accessor on this type read only the lower triangle, so
/// the upper-triangle bytes are dead memory and the memset
/// would be wasted bandwidth.
///
/// Halves memset traffic on pool reuse compared to
/// `buf.clear(); buf.resize(n*n, 0.0)`. Diagnosed as the
/// bandwidth bottleneck on c-big-class matrices in the
/// `solver_parallel_threadcount_sweep` run on 2026-05-12
/// (sp@8 plateaued at 1.10× before this change).
///
/// Safety contract: callers must not depend on the strict
/// upper triangle being zero. Audited readers as of 2026-05-12:
/// kernel only touches lower triangle (`dense/factor.rs:1138-1140`),
/// `extend_add` normalises to lower (`numeric/factorize.rs:2371`),
/// `symv` and `validate` iterate only `i >= j`.
pub fn from_pooled_buf(n: usize, mut buf: Vec<f64>) -> Self {
let needed = n * n;
// `resize` grows the tail with zeros or truncates. It does
// NOT re-zero entries `[0, min(old_len, needed))`, which is
// where the bandwidth saving comes from in steady state
// (capacity already reached after the first few large
// supernodes; subsequent calls only re-zero the lower
// triangle).
buf.resize(needed, 0.0);
// Zero only the lower triangle of the n×n column-major layout:
// for each column j, positions `[j*n + j, j*n + n)`.
for j in 0..n {
let col_base = j * n;
buf[col_base + j..col_base + n].fill(0.0);
}
Self { n, data: buf }
}
/// Create a symmetric matrix from a flat column-major vector.
/// The lower triangle is authoritative; the upper triangle is ignored.
pub fn from_column_major(n: usize, data: Vec<f64>) -> Result<Self, FeralError> {
if data.len() != n * n {
return Err(FeralError::InvalidInput(format!(
"matrix data length {} != expected {} for n={}",
data.len(),
n * n,
n
)));
}
Ok(Self { n, data })
}
/// Create a symmetric matrix from a dense 2D lower-triangular representation.
/// `entries` provides (i, j, value) triples where i >= j.
pub fn from_lower_triangle(n: usize, entries: &[(usize, usize, f64)]) -> Self {
let mut mat = Self::zeros(n);
for &(i, j, v) in entries {
mat.set(i, j, v);
}
mat
}
/// Get entry (i, j), reading from lower triangle.
/// For i >= j, returns data[j*n + i].
/// For i < j, returns data[i*n + j] (symmetric).
#[inline]
pub fn get(&self, i: usize, j: usize) -> f64 {
if i >= j {
self.data[j * self.n + i]
} else {
self.data[i * self.n + j]
}
}
/// Set entry (i, j) in the lower triangle.
/// Also sets (j, i) for symmetry in the stored data.
#[inline]
pub fn set(&mut self, i: usize, j: usize, val: f64) {
if i >= j {
self.data[j * self.n + i] = val;
} else {
self.data[i * self.n + j] = val;
}
}
/// Validate the matrix for factorization input.
/// Checks: n > 0, data length, no NaN/Inf in lower triangle.
pub fn validate(&self) -> Result<(), FeralError> {
if self.n == 0 {
return Err(FeralError::InvalidInput(
"matrix dimension is zero".to_string(),
));
}
if self.data.len() != self.n * self.n {
return Err(FeralError::InvalidInput(format!(
"matrix data length {} != expected {} for n={}",
self.data.len(),
self.n * self.n,
self.n
)));
}
// Check lower triangle for NaN/Inf
for j in 0..self.n {
for i in j..self.n {
let val = self.data[j * self.n + i];
if val.is_nan() || val.is_infinite() {
return Err(FeralError::InvalidInput(format!(
"matrix contains NaN or Inf at index ({},{})",
i, j
)));
}
}
}
Ok(())
}
/// Symmetric matrix-vector product: y = A * x.
/// Uses only the lower triangle.
pub fn symv(&self, x: &[f64], y: &mut [f64]) {
let n = self.n;
for yi in y.iter_mut().take(n) {
*yi = 0.0;
}
for j in 0..n {
// Diagonal
y[j] += self.data[j * n + j] * x[j];
// Off-diagonal (lower triangle)
for i in (j + 1)..n {
let a_ij = self.data[j * n + i];
y[i] += a_ij * x[j];
y[j] += a_ij * x[i];
}
}
}
}