la_stack/matrix.rs
1//! Fixed-size, stack-allocated square matrices.
2
3use core::hint::cold_path;
4
5use crate::ldlt::Ldlt;
6use crate::lu::Lu;
7use crate::{ERR_COEFF_2, ERR_COEFF_3, ERR_COEFF_4, LDLT_SYMMETRY_REL_TOL, LaError, Tolerance};
8
9/// Fixed-size square matrix `D×D`, stored inline.
10#[must_use]
11#[derive(Clone, Copy, Debug, PartialEq)]
12pub struct Matrix<const D: usize> {
13 pub(crate) rows: [[f64; D]; D],
14}
15
16/// Fixed-size square matrix whose stored entries are all finite.
17///
18/// This proof-bearing wrapper makes the finite invariant explicit for internal
19/// algorithms that should not repeatedly check stored entries for NaN or
20/// infinity.
21#[must_use]
22#[derive(Clone, Copy, Debug, PartialEq)]
23#[allow(clippy::redundant_pub_crate)]
24pub(crate) struct FiniteMatrix<const D: usize> {
25 matrix: Matrix<D>,
26}
27
28impl<const D: usize> FiniteMatrix<D> {
29 /// Construct a finite matrix without checking the invariant.
30 ///
31 /// This crate-internal escape hatch is only for paths with a local proof
32 /// that every stored entry is finite.
33 #[inline]
34 pub(crate) const fn new_unchecked(matrix: Matrix<D>) -> Self {
35 Self { matrix }
36 }
37
38 /// Validate a matrix and wrap it for algorithms that carry the finite
39 /// invariant explicitly.
40 ///
41 /// # Errors
42 /// Returns [`LaError::NonFinite`] with matrix coordinates for the first
43 /// stored NaN or infinity in row-major order.
44 #[inline]
45 pub(crate) const fn new(matrix: Matrix<D>) -> Result<Self, LaError> {
46 if let Some((row, col)) = Matrix::<D>::first_non_finite_cell_in(&matrix.rows) {
47 Err(LaError::non_finite_cell(row, col))
48 } else {
49 Ok(Self::new_unchecked(matrix))
50 }
51 }
52
53 /// Validate raw row-major storage and construct a finite matrix.
54 /// # Errors
55 /// Returns [`LaError::NonFinite`] with matrix coordinates for the first
56 /// offending entry in row-major order when `rows` contains NaN or infinity.
57 #[inline]
58 pub(crate) const fn from_rows(rows: [[f64; D]; D]) -> Result<Self, LaError> {
59 match Matrix::try_from_rows(rows) {
60 Ok(matrix) => Ok(Self::new_unchecked(matrix)),
61 Err(err) => Err(err),
62 }
63 }
64
65 /// All-zeros finite matrix.
66 #[inline]
67 pub(crate) const fn zero() -> Self {
68 Self::new_unchecked(Matrix::zero())
69 }
70
71 /// Consume the wrapper and return the underlying raw matrix.
72 #[inline]
73 pub(crate) const fn into_matrix(self) -> Matrix<D> {
74 self.matrix
75 }
76
77 /// Borrow the underlying raw matrix without revalidating stored entries.
78 #[cfg(feature = "exact")]
79 #[inline]
80 pub(crate) const fn as_matrix(&self) -> &Matrix<D> {
81 &self.matrix
82 }
83
84 /// Infinity norm (maximum absolute row sum) for a finite matrix.
85 ///
86 /// Stored entries are known finite, so this path only checks whether a row
87 /// sum overflows to NaN or infinity.
88 ///
89 /// # Errors
90 /// Returns [`LaError::NonFinite`] when a row sum overflows to NaN or infinity.
91 #[inline]
92 pub(crate) const fn inf_norm(&self) -> Result<f64, LaError> {
93 let mut max_row_sum: f64 = 0.0;
94
95 let mut r = 0;
96 while r < D {
97 let row = &self.matrix.rows[r];
98 let mut row_sum: f64 = 0.0;
99 let mut c = 0;
100 while c < D {
101 row_sum += row[c].abs();
102 c += 1;
103 }
104 if !row_sum.is_finite() {
105 cold_path();
106 return Err(LaError::non_finite_at(r));
107 }
108 if row_sum > max_row_sum {
109 max_row_sum = row_sum;
110 }
111 r += 1;
112 }
113
114 Ok(max_row_sum)
115 }
116
117 /// Returns `true` if the finite matrix is symmetric within a relative tolerance.
118 ///
119 /// # Errors
120 /// Returns [`LaError::NonFinite`] when computing the scaled symmetry tolerance
121 /// overflows to NaN or infinity.
122 #[inline]
123 pub(crate) fn is_symmetric(&self, rel_tol: Tolerance) -> Result<bool, LaError> {
124 Ok(self.first_asymmetry(rel_tol)?.is_none())
125 }
126
127 /// Returns the first asymmetric off-diagonal pair, if any.
128 ///
129 /// # Errors
130 /// Returns [`LaError::NonFinite`] when computing the scaled symmetry tolerance
131 /// overflows to NaN or infinity.
132 #[inline]
133 pub(crate) fn first_asymmetry(
134 &self,
135 rel_tol: Tolerance,
136 ) -> Result<Option<(usize, usize)>, LaError> {
137 let eps = self.symmetry_epsilon(rel_tol)?;
138 for r in 0..D {
139 for c in (r + 1)..D {
140 let upper = self.matrix.rows[r][c];
141 let lower = self.matrix.rows[c][r];
142
143 let diff = (upper - lower).abs();
144 if !diff.is_finite() || diff > eps {
145 cold_path();
146 return Ok(Some((r, c)));
147 }
148 }
149 }
150 Ok(None)
151 }
152
153 /// Compute the symmetry tolerance scale for a finite matrix.
154 fn symmetry_epsilon(&self, rel_tol: Tolerance) -> Result<f64, LaError> {
155 let rel_tol = rel_tol.get();
156 let mut eps = rel_tol;
157
158 for r in 0..D {
159 let mut row_eps = 0.0;
160 for c in 0..D {
161 row_eps = rel_tol.mul_add(self.matrix.rows[r][c].abs(), row_eps);
162 if !row_eps.is_finite() {
163 cold_path();
164 return Err(LaError::non_finite_at(c));
165 }
166 }
167 if row_eps > eps {
168 eps = row_eps;
169 }
170 }
171
172 Ok(eps)
173 }
174
175 /// Compute an LU decomposition with partial pivoting.
176 ///
177 /// # Errors
178 /// Returns [`LaError::Singular`] for an unusable pivot, or
179 /// [`LaError::NonFinite`] if elimination computes a non-finite intermediate.
180 #[inline]
181 pub(crate) fn lu(self, tol: Tolerance) -> Result<Lu<D>, LaError> {
182 Lu::factor_finite(self, tol)
183 }
184
185 /// Compute an LDLT factorization (`A = L D Lᵀ`) without pivoting.
186 ///
187 /// # Errors
188 /// Returns [`LaError::Asymmetric`] if the matrix is not symmetric,
189 /// [`LaError::NotPositiveSemidefinite`] for a negative LDLT diagonal,
190 /// [`LaError::Singular`] for a zero or too-small non-negative diagonal, or
191 /// [`LaError::NonFinite`] if factorization computes a non-finite intermediate.
192 #[inline]
193 pub(crate) fn ldlt(self, tol: Tolerance) -> Result<Ldlt<D>, LaError> {
194 Ldlt::factor_symmetric(SymmetricMatrix::try_new(self)?, tol)
195 }
196
197 /// Closed-form determinant for dimensions 0–4, bypassing LU factorization.
198 ///
199 /// # Errors
200 /// Returns [`LaError::NonFinite`] when the closed-form determinant overflows
201 /// to NaN or infinity.
202 #[inline]
203 pub(crate) const fn det_direct(&self) -> Result<Option<f64>, LaError> {
204 match D {
205 0 => Ok(Some(1.0)),
206 1 => Self::computed_scalar_result(Some(self.matrix.rows[0][0])),
207 2 => Self::computed_scalar_result(Some(self.matrix.rows[0][0].mul_add(
208 self.matrix.rows[1][1],
209 -(self.matrix.rows[0][1] * self.matrix.rows[1][0]),
210 ))),
211 3 => {
212 let m00 = self.matrix.rows[1][1].mul_add(
213 self.matrix.rows[2][2],
214 -(self.matrix.rows[1][2] * self.matrix.rows[2][1]),
215 );
216 let m01 = self.matrix.rows[1][0].mul_add(
217 self.matrix.rows[2][2],
218 -(self.matrix.rows[1][2] * self.matrix.rows[2][0]),
219 );
220 let m02 = self.matrix.rows[1][0].mul_add(
221 self.matrix.rows[2][1],
222 -(self.matrix.rows[1][1] * self.matrix.rows[2][0]),
223 );
224 Self::computed_scalar_result(Some(self.matrix.rows[0][0].mul_add(
225 m00,
226 (-self.matrix.rows[0][1]).mul_add(m01, self.matrix.rows[0][2] * m02),
227 )))
228 }
229 4 => {
230 let r = &self.matrix.rows;
231
232 let s23 = r[2][2].mul_add(r[3][3], -(r[2][3] * r[3][2]));
233 let s13 = r[2][1].mul_add(r[3][3], -(r[2][3] * r[3][1]));
234 let s12 = r[2][1].mul_add(r[3][2], -(r[2][2] * r[3][1]));
235 let s03 = r[2][0].mul_add(r[3][3], -(r[2][3] * r[3][0]));
236 let s02 = r[2][0].mul_add(r[3][2], -(r[2][2] * r[3][0]));
237 let s01 = r[2][0].mul_add(r[3][1], -(r[2][1] * r[3][0]));
238
239 let c00 = r[1][1].mul_add(s23, (-r[1][2]).mul_add(s13, r[1][3] * s12));
240 let c01 = r[1][0].mul_add(s23, (-r[1][2]).mul_add(s03, r[1][3] * s02));
241 let c02 = r[1][0].mul_add(s13, (-r[1][1]).mul_add(s03, r[1][3] * s01));
242 let c03 = r[1][0].mul_add(s12, (-r[1][1]).mul_add(s02, r[1][2] * s01));
243
244 Self::computed_scalar_result(Some(r[0][0].mul_add(
245 c00,
246 (-r[0][1]).mul_add(c01, r[0][2].mul_add(c02, -(r[0][3] * c03))),
247 )))
248 }
249 _ => {
250 cold_path();
251 Ok(None)
252 }
253 }
254 }
255
256 /// Floating-point determinant for a finite matrix.
257 ///
258 /// # Errors
259 /// Returns [`LaError::NonFinite`] if a computed determinant or LU fallback
260 /// intermediate overflows to NaN or infinity.
261 #[inline]
262 pub(crate) fn det(self) -> Result<f64, LaError> {
263 if let Some(d) = self.det_direct()? {
264 return Ok(d);
265 }
266 match self.lu(Tolerance::new_unchecked(0.0)) {
267 Ok(lu) => lu.det(),
268 Err(LaError::Singular { .. }) => Ok(0.0),
269 Err(err) => Err(err),
270 }
271 }
272
273 /// Conservative absolute error bound for [`det_direct`](Self::det_direct).
274 ///
275 /// # Errors
276 /// Returns [`LaError::NonFinite`] when the bound computation overflows to NaN
277 /// or infinity.
278 #[inline]
279 pub(crate) const fn det_errbound(&self) -> Result<Option<f64>, LaError> {
280 match D {
281 0 | 1 => Self::computed_scalar_result(Some(0.0)),
282 2 => {
283 let r = &self.matrix.rows;
284 let permanent = (r[0][0] * r[1][1]).abs() + (r[0][1] * r[1][0]).abs();
285 Self::computed_scalar_result(Some(ERR_COEFF_2 * permanent))
286 }
287 3 => {
288 let r = &self.matrix.rows;
289 let pm00 = (r[1][1] * r[2][2]).abs() + (r[1][2] * r[2][1]).abs();
290 let pm01 = (r[1][0] * r[2][2]).abs() + (r[1][2] * r[2][0]).abs();
291 let pm02 = (r[1][0] * r[2][1]).abs() + (r[1][1] * r[2][0]).abs();
292 let permanent = r[0][2]
293 .abs()
294 .mul_add(pm02, r[0][1].abs().mul_add(pm01, r[0][0].abs() * pm00));
295 Self::computed_scalar_result(Some(ERR_COEFF_3 * permanent))
296 }
297 4 => {
298 let r = &self.matrix.rows;
299 let sp23 = (r[2][2] * r[3][3]).abs() + (r[2][3] * r[3][2]).abs();
300 let sp13 = (r[2][1] * r[3][3]).abs() + (r[2][3] * r[3][1]).abs();
301 let sp12 = (r[2][1] * r[3][2]).abs() + (r[2][2] * r[3][1]).abs();
302 let sp03 = (r[2][0] * r[3][3]).abs() + (r[2][3] * r[3][0]).abs();
303 let sp02 = (r[2][0] * r[3][2]).abs() + (r[2][2] * r[3][0]).abs();
304 let sp01 = (r[2][0] * r[3][1]).abs() + (r[2][1] * r[3][0]).abs();
305 let pc0 = r[1][3]
306 .abs()
307 .mul_add(sp12, r[1][2].abs().mul_add(sp13, r[1][1].abs() * sp23));
308 let pc1 = r[1][3]
309 .abs()
310 .mul_add(sp02, r[1][2].abs().mul_add(sp03, r[1][0].abs() * sp23));
311 let pc2 = r[1][3]
312 .abs()
313 .mul_add(sp01, r[1][1].abs().mul_add(sp03, r[1][0].abs() * sp13));
314 let pc3 = r[1][2]
315 .abs()
316 .mul_add(sp01, r[1][1].abs().mul_add(sp02, r[1][0].abs() * sp12));
317 let permanent = r[0][3].abs().mul_add(
318 pc3,
319 r[0][2]
320 .abs()
321 .mul_add(pc2, r[0][1].abs().mul_add(pc1, r[0][0].abs() * pc0)),
322 );
323 Self::computed_scalar_result(Some(ERR_COEFF_4 * permanent))
324 }
325 _ => {
326 cold_path();
327 Ok(None)
328 }
329 }
330 }
331
332 /// Return a computed scalar result for a matrix with finite stored entries.
333 const fn computed_scalar_result(value: Option<f64>) -> Result<Option<f64>, LaError> {
334 match value {
335 Some(value) if value.is_finite() => Ok(Some(value)),
336 Some(_) => Err(LaError::non_finite_at(0)),
337 None => Ok(None),
338 }
339 }
340}
341
342impl<const D: usize> Default for FiniteMatrix<D> {
343 #[inline]
344 fn default() -> Self {
345 Self::zero()
346 }
347}
348
349impl<const D: usize> From<FiniteMatrix<D>> for Matrix<D> {
350 #[inline]
351 fn from(value: FiniteMatrix<D>) -> Self {
352 value.into_matrix()
353 }
354}
355
356impl<const D: usize> TryFrom<Matrix<D>> for FiniteMatrix<D> {
357 type Error = LaError;
358
359 #[inline]
360 fn try_from(value: Matrix<D>) -> Result<Self, Self::Error> {
361 Self::new(value)
362 }
363}
364
365impl<const D: usize> TryFrom<[[f64; D]; D]> for FiniteMatrix<D> {
366 type Error = LaError;
367
368 #[inline]
369 fn try_from(value: [[f64; D]; D]) -> Result<Self, Self::Error> {
370 Self::from_rows(value)
371 }
372}
373
374/// Matrix proven finite and symmetric under the crate's LDLT symmetry tolerance.
375#[must_use]
376#[derive(Clone, Copy, Debug, PartialEq)]
377#[allow(clippy::redundant_pub_crate)]
378pub(crate) struct SymmetricMatrix<const D: usize> {
379 matrix: Matrix<D>,
380}
381
382impl<const D: usize> SymmetricMatrix<D> {
383 /// Construct a symmetric matrix proof without checking the invariant.
384 ///
385 /// This constructor is only for paths that have already validated finite
386 /// entries and LDLT symmetry with the same predicate as
387 /// [`try_new`](Self::try_new).
388 #[inline]
389 pub(crate) const fn new_unchecked(matrix: Matrix<D>) -> Self {
390 Self { matrix }
391 }
392
393 /// Validate that a finite matrix is symmetric under the LDLT symmetry tolerance.
394 ///
395 /// The predicate is the same one used by [`Matrix::ldlt`]:
396 /// `|A[i][j] - A[j][i]| <= 1e-12 * max(1, inf_norm(A))`, with scaling that
397 /// preserves strict tolerances when an unscaled row sum would overflow.
398 ///
399 /// # Errors
400 /// Returns [`LaError::Asymmetric`] when the first off-diagonal pair violates
401 /// the LDLT symmetry predicate.
402 ///
403 /// Returns [`LaError::NonFinite`] when computing the scaled symmetry
404 /// tolerance overflows to NaN or infinity.
405 #[inline]
406 pub(crate) fn try_new(matrix: FiniteMatrix<D>) -> Result<Self, LaError> {
407 if let Some((row, col)) = matrix.first_asymmetry(LDLT_SYMMETRY_REL_TOL)? {
408 cold_path();
409 Err(LaError::asymmetric(row, col, D))
410 } else {
411 Ok(Self::new_unchecked(matrix.into_matrix()))
412 }
413 }
414
415 /// Consume the wrapper and return the underlying matrix.
416 #[inline]
417 pub(crate) const fn into_matrix(self) -> Matrix<D> {
418 self.matrix
419 }
420}
421
422impl<const D: usize> Matrix<D> {
423 /// Test-only infallible constructor for finite literal fixtures.
424 #[cfg(test)]
425 #[inline]
426 pub(crate) const fn from_rows(rows: [[f64; D]; D]) -> Self {
427 match Self::try_from_rows(rows) {
428 Ok(matrix) => matrix,
429 Err(_) => panic!("Matrix::from_rows requires finite entries"),
430 }
431 }
432
433 /// Try to create a finite matrix from row-major storage.
434 ///
435 /// This is the public raw-storage boundary for matrices. Public compute
436 /// methods parse stored rows into crate-internal proof-bearing types before
437 /// arithmetic, including when crate-internal unchecked storage exists.
438 ///
439 /// # Examples
440 /// ```
441 /// use la_stack::prelude::*;
442 ///
443 /// # fn main() -> Result<(), LaError> {
444 /// let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
445 /// assert_eq!(m.get(0, 1), Some(2.0));
446 /// # Ok(())
447 /// # }
448 /// ```
449 ///
450 /// # Errors
451 /// Returns [`LaError::NonFinite`] with matrix coordinates for the first
452 /// offending entry in row-major order when `rows` contains NaN or infinity.
453 #[inline]
454 pub const fn try_from_rows(rows: [[f64; D]; D]) -> Result<Self, LaError> {
455 if let Some((row, col)) = Self::first_non_finite_cell_in(&rows) {
456 Err(LaError::non_finite_cell(row, col))
457 } else {
458 Ok(Self::from_rows_unchecked(rows))
459 }
460 }
461
462 /// Construct a matrix without checking that entries are finite.
463 ///
464 /// This crate-internal escape hatch is reserved for literals and algorithm
465 /// outputs whose finite invariant is visible at the call site.
466 #[inline]
467 pub(crate) const fn from_rows_unchecked(rows: [[f64; D]; D]) -> Self {
468 Self { rows }
469 }
470
471 /// All-zeros matrix.
472 ///
473 /// # Examples
474 /// ```
475 /// use la_stack::prelude::*;
476 ///
477 /// let z = Matrix::<2>::zero();
478 /// assert_eq!(z.get(1, 1), Some(0.0));
479 /// ```
480 #[inline]
481 pub const fn zero() -> Self {
482 Self::from_rows_unchecked([[0.0; D]; D])
483 }
484
485 /// Identity matrix.
486 ///
487 /// # Examples
488 /// ```
489 /// use la_stack::prelude::*;
490 ///
491 /// let i = Matrix::<3>::identity();
492 /// assert_eq!(i.get(0, 0), Some(1.0));
493 /// assert_eq!(i.get(0, 1), Some(0.0));
494 /// assert_eq!(i.get(2, 2), Some(1.0));
495 /// ```
496 #[inline]
497 pub const fn identity() -> Self {
498 let mut m = Self::zero();
499
500 let mut i = 0;
501 while i < D {
502 m.rows[i][i] = 1.0;
503 i += 1;
504 }
505
506 m
507 }
508
509 /// Get an element with bounds checking.
510 ///
511 /// # Examples
512 /// ```
513 /// use la_stack::prelude::*;
514 ///
515 /// # fn main() -> Result<(), LaError> {
516 /// let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
517 /// assert_eq!(m.get(1, 0), Some(3.0));
518 /// assert_eq!(m.get(2, 0), None);
519 /// # Ok(())
520 /// # }
521 /// ```
522 #[inline]
523 #[must_use]
524 pub const fn get(&self, r: usize, c: usize) -> Option<f64> {
525 if r < D && c < D {
526 Some(self.rows[r][c])
527 } else {
528 None
529 }
530 }
531
532 /// Get an element, preserving index context on failure.
533 ///
534 /// Prefer [`get`](Self::get) for const or hot paths that only need
535 /// `Option`-style absence. Use this method at public runtime boundaries
536 /// where row, column, and dimension context should survive in a typed error.
537 ///
538 /// # Examples
539 /// ```
540 /// use la_stack::prelude::*;
541 ///
542 /// # fn main() -> Result<(), LaError> {
543 /// let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
544 /// assert_eq!(m.get_checked(1, 0)?, 3.0);
545 /// assert_eq!(
546 /// m.get_checked(2, 0),
547 /// Err(LaError::IndexOutOfBounds {
548 /// row: 2,
549 /// col: 0,
550 /// dim: 2,
551 /// })
552 /// );
553 /// # Ok(())
554 /// # }
555 /// ```
556 ///
557 /// # Errors
558 /// Returns [`LaError::IndexOutOfBounds`] when either index is not `< D`.
559 #[inline]
560 pub const fn get_checked(&self, row: usize, col: usize) -> Result<f64, LaError> {
561 if row < D && col < D {
562 Ok(self.rows[row][col])
563 } else {
564 Err(LaError::index_out_of_bounds(row, col, D))
565 }
566 }
567
568 /// Set a finite element with bounds checking.
569 ///
570 /// # Examples
571 /// ```
572 /// use la_stack::prelude::*;
573 ///
574 /// # fn main() -> Result<(), LaError> {
575 /// let mut m = Matrix::<2>::zero();
576 /// assert_eq!(m.set(0, 1, 2.5), Ok(()));
577 /// assert_eq!(m.get(0, 1), Some(2.5));
578 /// assert_eq!(
579 /// m.set(10, 0, 1.0),
580 /// Err(LaError::IndexOutOfBounds {
581 /// row: 10,
582 /// col: 0,
583 /// dim: 2,
584 /// })
585 /// );
586 /// # Ok(())
587 /// # }
588 /// ```
589 ///
590 /// # Errors
591 /// Returns [`LaError::IndexOutOfBounds`] when either index is not `< D`.
592 /// Returns [`LaError::NonFinite`] when `value` is NaN or infinity.
593 #[inline]
594 pub const fn set(&mut self, row: usize, col: usize, value: f64) -> Result<(), LaError> {
595 self.set_checked(row, col, value)
596 }
597
598 /// Set an element, preserving index context on failure.
599 ///
600 /// The matrix is mutated only when `(row, col)` is in bounds and `value` is
601 /// finite.
602 ///
603 /// # Examples
604 /// ```
605 /// use la_stack::prelude::*;
606 ///
607 /// # fn main() -> Result<(), LaError> {
608 /// let mut m = Matrix::<2>::zero();
609 /// m.set_checked(0, 1, 2.5)?;
610 /// assert_eq!(m.get_checked(0, 1)?, 2.5);
611 ///
612 /// assert_eq!(
613 /// m.set_checked(10, 0, 1.0),
614 /// Err(LaError::IndexOutOfBounds {
615 /// row: 10,
616 /// col: 0,
617 /// dim: 2,
618 /// })
619 /// );
620 /// # Ok(())
621 /// # }
622 /// ```
623 ///
624 /// # Errors
625 /// Returns [`LaError::IndexOutOfBounds`] when either index is not `< D`.
626 /// Returns [`LaError::NonFinite`] when `value` is NaN or infinity.
627 #[inline]
628 pub const fn set_checked(&mut self, row: usize, col: usize, value: f64) -> Result<(), LaError> {
629 if row >= D || col >= D {
630 return Err(LaError::index_out_of_bounds(row, col, D));
631 }
632 if !value.is_finite() {
633 return Err(LaError::non_finite_cell(row, col));
634 }
635 self.rows[row][col] = value;
636 Ok(())
637 }
638
639 /// Infinity norm (maximum absolute row sum).
640 ///
641 /// # Non-finite handling
642 /// Public constructors and setters reject raw non-finite entries, but
643 /// crate-internal unchecked storage can still contain NaN or infinity.
644 /// `inf_norm` returns [`LaError::NonFinite`] if it encounters stored NaN/∞
645 /// or if a row sum overflows to a non-finite value.
646 ///
647 /// Row sums are accumulated in `f64` with ordinary addition. This method
648 /// checks for overflowed accumulators, but it does not provide a certified
649 /// absolute rounding bound for the returned norm.
650 ///
651 /// # Examples
652 /// ```
653 /// use la_stack::prelude::*;
654 ///
655 /// # fn main() -> Result<(), LaError> {
656 /// let m = Matrix::<2>::try_from_rows([[1.0, -2.0], [3.0, 4.0]])?;
657 /// assert!((m.inf_norm()? - 7.0).abs() <= 1e-12);
658 ///
659 /// // Raw NaN entries are rejected with coordinates.
660 /// assert_eq!(
661 /// Matrix::<2>::try_from_rows([[f64::NAN, 1.0], [2.0, 3.0]]),
662 /// Err(LaError::NonFinite {
663 /// row: Some(0),
664 /// col: 0,
665 /// })
666 /// );
667 /// # Ok(())
668 /// # }
669 /// ```
670 ///
671 /// # Errors
672 /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or a
673 /// row sum overflows to NaN or infinity.
674 #[inline]
675 pub const fn inf_norm(&self) -> Result<f64, LaError> {
676 match FiniteMatrix::new(*self) {
677 Ok(matrix) => matrix.inf_norm(),
678 Err(err) => Err(err),
679 }
680 }
681
682 /// Returns `true` if the matrix is symmetric within a relative tolerance.
683 ///
684 /// Two entries `self[r][c]` and `self[c][r]` are considered equal (for the
685 /// purposes of symmetry) when
686 /// `|self[r][c] - self[c][r]| <= rel_tol * max(1.0, inf_norm(self))`.
687 /// This mirrors the predicate used internally by [`ldlt`](Self::ldlt), so
688 /// callers can pre-validate matrices that may come from untrusted sources.
689 ///
690 /// Use [`first_asymmetry`](Self::first_asymmetry) to locate the first
691 /// offending pair when this returns `Ok(false)`.
692 ///
693 /// The `rel_tol` argument is a [`Tolerance`], so raw caller input must be
694 /// finite and non-negative before it can reach this predicate. Use
695 /// [`Tolerance::new`] or [`LaError::validate_tolerance`] when accepting a
696 /// raw `f64`; negative, NaN, and infinite tolerances return
697 /// [`LaError::InvalidTolerance`].
698 ///
699 /// # Overflow handling
700 /// A finite matrix can return [`LaError::NonFinite`] if computing the scaled
701 /// symmetry tolerance overflows to NaN or infinity. If both stored entries
702 /// are finite but their difference overflows to ±∞, the pair is reported as
703 /// asymmetric.
704 ///
705 /// # Examples
706 /// ```
707 /// use la_stack::prelude::*;
708 ///
709 /// # fn main() -> Result<(), LaError> {
710 /// let a = Matrix::<2>::try_from_rows([[4.0, 2.0], [2.0, 3.0]])?;
711 /// let tol = Tolerance::new(1e-12)?;
712 /// assert!(a.is_symmetric(tol)?);
713 ///
714 /// let b = Matrix::<2>::try_from_rows([[4.0, 2.0], [3.0, 3.0]])?;
715 /// assert!(!b.is_symmetric(tol)?);
716 /// # Ok(())
717 /// # }
718 /// ```
719 ///
720 /// # Errors
721 /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or
722 /// computing the scaled symmetry tolerance overflows to NaN or infinity.
723 #[inline]
724 pub fn is_symmetric(&self, rel_tol: Tolerance) -> Result<bool, LaError> {
725 FiniteMatrix::new(*self)?.is_symmetric(rel_tol)
726 }
727
728 /// Returns the indices `(r, c)` (with `r < c`) of the first off-diagonal
729 /// pair that violates symmetry, or `None` if the matrix is symmetric
730 /// within `rel_tol`.
731 ///
732 /// Iteration order is row-major over the strict upper triangle, so the
733 /// returned indices are the lexicographically smallest such pair. The
734 /// predicate is the same as [`is_symmetric`](Self::is_symmetric):
735 /// `|self[r][c] - self[c][r]| <= rel_tol * max(1.0, inf_norm(self))`.
736 ///
737 /// A finite matrix can return [`LaError::NonFinite`] if computing the scaled
738 /// symmetry tolerance overflows to NaN or infinity. If both stored entries
739 /// are finite but their difference overflows to ±∞, the pair is reported as
740 /// asymmetric.
741 ///
742 /// The `rel_tol` argument is a [`Tolerance`], so raw caller input must be
743 /// finite and non-negative before it can reach this predicate. Use
744 /// [`Tolerance::new`] or [`LaError::validate_tolerance`] when accepting a
745 /// raw `f64`; negative, NaN, and infinite tolerances return
746 /// [`LaError::InvalidTolerance`].
747 ///
748 /// # Examples
749 /// ```
750 /// use la_stack::prelude::*;
751 ///
752 /// # fn main() -> Result<(), LaError> {
753 /// let a = Matrix::<3>::try_from_rows([
754 /// [1.0, 2.0, 0.0],
755 /// [2.0, 4.0, 5.0],
756 /// [0.0, 6.0, 9.0], // 6.0 breaks symmetry with a[1][2] = 5.0
757 /// ])?;
758 /// let tol = Tolerance::new(1e-12)?;
759 /// assert_eq!(a.first_asymmetry(tol)?, Some((1, 2)));
760 /// assert_eq!(Matrix::<3>::identity().first_asymmetry(tol)?, None);
761 /// # Ok(())
762 /// # }
763 /// ```
764 ///
765 /// # Errors
766 /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or
767 /// computing the scaled symmetry tolerance overflows to NaN or infinity.
768 #[inline]
769 pub fn first_asymmetry(&self, rel_tol: Tolerance) -> Result<Option<(usize, usize)>, LaError> {
770 FiniteMatrix::new(*self)?.first_asymmetry(rel_tol)
771 }
772
773 /// Compute an LU decomposition with partial pivoting.
774 ///
775 /// # Examples
776 /// ```
777 /// use la_stack::prelude::*;
778 ///
779 /// # fn main() -> Result<(), LaError> {
780 /// let a = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
781 /// let lu = a.lu(DEFAULT_PIVOT_TOL)?;
782 ///
783 /// let b = Vector::<2>::try_new([5.0, 11.0])?;
784 /// let x = lu.solve_vec(b)?.into_array();
785 ///
786 /// assert!((x[0] - 1.0).abs() <= 1e-12);
787 /// assert!((x[1] - 2.0).abs() <= 1e-12);
788 /// # Ok(())
789 /// # }
790 /// ```
791 ///
792 /// The `tol` argument is a [`Tolerance`], so raw caller input must be
793 /// finite and non-negative before it can reach factorization. Use
794 /// [`Tolerance::new`] or [`LaError::validate_tolerance`] when accepting a
795 /// raw `f64`; negative, NaN, and infinite tolerances return
796 /// [`LaError::InvalidTolerance`].
797 ///
798 /// # Errors
799 /// Returns [`LaError::Singular`] if, for some column `k`, the largest-magnitude candidate pivot
800 /// in that column satisfies `|pivot| <= tol` (so no numerically usable pivot exists).
801 /// Returns [`LaError::NonFinite`] if stored entries are NaN/infinity or an
802 /// elimination intermediate overflows to NaN/∞ before it can be stored in
803 /// the returned [`Lu`].
804 #[inline]
805 pub fn lu(self, tol: Tolerance) -> Result<Lu<D>, LaError> {
806 FiniteMatrix::new(self)?.lu(tol)
807 }
808
809 /// Compute an LDLT factorization (`A = L D Lᵀ`) without pivoting.
810 ///
811 /// This is intended for symmetric positive definite (SPD) and positive semi-definite (PSD)
812 /// matrices such as Gram matrices.
813 ///
814 /// # Symmetry validation
815 /// The input matrix `self` must be symmetric — that is,
816 /// `self[i][j] == self[j][i]` within the crate's LDLT symmetry tolerance
817 /// (`1e-12`, scaled like [`is_symmetric`](Self::is_symmetric)). This is a
818 /// correctness invariant, not merely a performance hint, so asymmetric inputs return
819 /// [`LaError::Asymmetric`] before factorization starts. If you need a
820 /// general-purpose factorization that tolerates non-symmetric inputs, use
821 /// [`lu`](Self::lu) instead.
822 ///
823 /// The `tol` argument is a [`Tolerance`], so raw caller input must be
824 /// finite and non-negative before it can reach factorization. Use
825 /// [`Tolerance::new`] or [`LaError::validate_tolerance`] when accepting a
826 /// raw `f64`; negative, NaN, and infinite tolerances return
827 /// [`LaError::InvalidTolerance`].
828 ///
829 /// # Examples
830 /// ```
831 /// use la_stack::prelude::*;
832 ///
833 /// # fn main() -> Result<(), LaError> {
834 /// // Note the symmetric layout: a[0][1] == a[1][0] == 2.0.
835 /// let a = Matrix::<2>::try_from_rows([[4.0, 2.0], [2.0, 3.0]])?;
836 /// let ldlt = a.ldlt(DEFAULT_SINGULAR_TOL)?;
837 ///
838 /// // det(A) = 8
839 /// assert!((ldlt.det()? - 8.0).abs() <= 1e-12);
840 ///
841 /// // Solve A x = b
842 /// let b = Vector::<2>::try_new([1.0, 2.0])?;
843 /// let x = ldlt.solve_vec(b)?.into_array();
844 /// assert!((x[0] - (-0.125)).abs() <= 1e-12);
845 /// assert!((x[1] - 0.75).abs() <= 1e-12);
846 /// # Ok(())
847 /// # }
848 /// ```
849 ///
850 /// # Errors
851 /// Returns [`LaError::NotPositiveSemidefinite`] if, for some step `k`, the required
852 /// diagonal entry `d = D[k,k]` is negative.
853 /// Returns [`LaError::Singular`] if `0 <= d <= tol`, treating PSD degeneracy
854 /// as singular/degenerate.
855 /// Returns [`LaError::NonFinite`] if stored entries are NaN/infinity or
856 /// factorization computes a non-finite intermediate.
857 /// Returns [`LaError::Asymmetric`] if the input matrix is not symmetric.
858 #[inline]
859 pub fn ldlt(self, tol: Tolerance) -> Result<Ldlt<D>, LaError> {
860 FiniteMatrix::new(self)?.ldlt(tol)
861 }
862
863 /// Return the first non-finite stored cell in row-major order.
864 const fn first_non_finite_cell_in(rows: &[[f64; D]; D]) -> Option<(usize, usize)> {
865 let mut r = 0;
866 while r < D {
867 let mut c = 0;
868 while c < D {
869 if !rows[r][c].is_finite() {
870 return Some((r, c));
871 }
872 c += 1;
873 }
874 r += 1;
875 }
876 None
877 }
878
879 /// Closed-form determinant for dimensions 0–4, bypassing LU factorization.
880 ///
881 /// Returns `Ok(Some(det))` for `D` ∈ {0, 1, 2, 3, 4}, `Ok(None)` for D ≥ 5.
882 /// `D = 0` returns `Ok(Some(1.0))` (empty product).
883 /// This is a `const fn` (Rust 1.94+) and uses fused multiply-add (`mul_add`)
884 /// for improved accuracy and performance.
885 ///
886 /// For a determinant that works for any dimension (falling back to LU for D ≥ 5),
887 /// use [`det`](Self::det).
888 ///
889 /// # Examples
890 /// ```
891 /// use la_stack::prelude::*;
892 ///
893 /// # fn main() -> Result<(), LaError> {
894 /// let m = Matrix::<2>::try_from_rows([[1.0, 2.0], [3.0, 4.0]])?;
895 /// assert_eq!(m.det_direct()?, Some(-2.0));
896 ///
897 /// // D = 0 is the empty product.
898 /// assert_eq!(Matrix::<0>::zero().det_direct()?, Some(1.0));
899 ///
900 /// // D ≥ 5 returns None.
901 /// assert!(Matrix::<5>::identity().det_direct()?.is_none());
902 /// # Ok(())
903 /// # }
904 /// ```
905 ///
906 /// # Errors
907 /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or
908 /// the closed-form determinant overflows to NaN or infinity.
909 #[inline]
910 pub const fn det_direct(&self) -> Result<Option<f64>, LaError> {
911 match FiniteMatrix::new(*self) {
912 Ok(matrix) => matrix.det_direct(),
913 Err(err) => Err(err),
914 }
915 }
916
917 /// Floating-point determinant, using closed-form formulas for D ≤ 4 and
918 /// LU decomposition for D ≥ 5.
919 ///
920 /// For D ∈ {1, 2, 3, 4}, this bypasses LU factorization entirely for a significant
921 /// speedup (see [`det_direct`](Self::det_direct)).
922 ///
923 /// Finite inputs return a floating-point determinant estimate in every dimension;
924 /// this method does not surface [`LaError::Singular`]. Because it mixes
925 /// closed-form paths from [`det_direct`](Self::det_direct) with an LU fallback,
926 /// the returned value has no certified absolute error bound. Use
927 /// [`det_errbound`](Self::det_errbound) for D ≤ 4 bounds, or the exact
928 /// determinant APIs when exact singularity classification or certified values
929 /// matter. For D ≥ 5, the LU fallback only maps an exactly zero pivot to
930 /// `Ok(0.0)`. Use [`lu`](Self::lu) directly when you need tolerance-aware
931 /// singularity detection or the pivot column.
932 ///
933 /// # Examples
934 /// ```
935 /// use la_stack::prelude::*;
936 ///
937 /// # fn main() -> Result<(), LaError> {
938 /// let det = Matrix::<3>::identity().det()?;
939 /// assert!((det - 1.0).abs() <= 1e-12);
940 /// # Ok(())
941 /// # }
942 /// ```
943 ///
944 /// # Errors
945 /// Returns [`LaError::NonFinite`] if stored entries are NaN/infinity, the
946 /// LU fallback computes a non-finite factorization cell, or the determinant
947 /// product overflows to NaN or infinity.
948 #[inline]
949 pub fn det(self) -> Result<f64, LaError> {
950 FiniteMatrix::new(self)?.det()
951 }
952
953 /// Conservative absolute error bound for `det_direct()`.
954 ///
955 /// Returns `Ok(Some(bound))` such that `|det_direct() - det_exact| ≤ bound`,
956 /// or `Ok(None)` for D ≥ 5 where no fast bound is available.
957 ///
958 /// For D ≤ 4, the bound is derived from the absolute Leibniz sum using
959 /// Shewchuk-style error analysis (see `REFERENCES.md` \[8\] and the
960 /// per-constant docs on [`ERR_COEFF_2`], [`ERR_COEFF_3`], and
961 /// [`ERR_COEFF_4`]). For D = 0 or 1, returns
962 /// `Some(0.0)` since the determinant computation is exact (no
963 /// arithmetic).
964 ///
965 /// This method does NOT require the `exact` feature — the bounds use
966 /// pure f64 arithmetic and are useful for custom adaptive-precision logic.
967 ///
968 /// # When to use
969 ///
970 /// Use this to build adaptive-precision logic: if `|det_direct()| > bound`,
971 /// the f64 sign is provably correct. Otherwise fall back to exact arithmetic.
972 ///
973 /// # Examples
974 /// ```
975 /// use la_stack::prelude::*;
976 ///
977 /// # fn main() -> Result<(), LaError> {
978 /// let m = Matrix::<3>::try_from_rows([
979 /// [1.0, 2.0, 3.0],
980 /// [4.0, 5.0, 6.0],
981 /// [7.0, 8.0, 9.0],
982 /// ])?;
983 /// if let (Some(bound), Some(det_approx)) = (m.det_errbound()?, m.det_direct()?) {
984 /// // If |det_approx| > bound, the sign is guaranteed correct.
985 /// let sign_is_certified = det_approx.abs() > bound;
986 /// assert!(!sign_is_certified);
987 /// }
988 /// # Ok(())
989 /// # }
990 /// ```
991 ///
992 /// # Adaptive precision pattern (requires `exact` feature)
993 /// ```ignore
994 /// use la_stack::prelude::*;
995 ///
996 /// let m = Matrix::<3>::identity();
997 /// if let Some(bound) = m.det_errbound()? {
998 /// if let Some(det) = m.det_direct()? {
999 /// if det.abs() > bound {
1000 /// // f64 sign is guaranteed correct
1001 /// let sign = det.signum() as i8;
1002 /// } else {
1003 /// // Fall back to exact arithmetic (requires `exact` feature)
1004 /// let sign = m.det_sign_exact()?;
1005 /// }
1006 /// }
1007 /// } else {
1008 /// // D ≥ 5: no fast filter, use exact directly
1009 /// let sign = m.det_sign_exact()?;
1010 /// }
1011 /// ```
1012 ///
1013 /// # Errors
1014 /// Returns [`LaError::NonFinite`] when stored entries are NaN/infinity or
1015 /// the bound computation overflows to NaN or infinity.
1016 #[inline]
1017 pub const fn det_errbound(&self) -> Result<Option<f64>, LaError> {
1018 match FiniteMatrix::new(*self) {
1019 Ok(matrix) => matrix.det_errbound(),
1020 Err(err) => Err(err),
1021 }
1022 }
1023}
1024
1025impl<const D: usize> Default for Matrix<D> {
1026 #[inline]
1027 fn default() -> Self {
1028 Self::zero()
1029 }
1030}
1031
1032#[cfg(test)]
1033mod tests {
1034 use super::*;
1035 use crate::DEFAULT_PIVOT_TOL;
1036 use crate::vector::{FiniteVector, Vector};
1037
1038 use approx::assert_abs_diff_eq;
1039 use core::assert_matches;
1040 use pastey::paste;
1041 use std::hint::black_box;
1042
1043 fn assert_matrix_abs_eq<const D: usize>(actual: &Matrix<D>, expected: &Matrix<D>) {
1044 for r in 0..D {
1045 for c in 0..D {
1046 assert_abs_diff_eq!(
1047 actual.get(r, c).unwrap(),
1048 expected.get(r, c).unwrap(),
1049 epsilon = 0.0
1050 );
1051 }
1052 }
1053 }
1054
1055 fn assert_array_abs_eq<const D: usize>(actual: &[f64; D], expected: &[f64; D]) {
1056 for i in 0..D {
1057 assert_abs_diff_eq!(actual[i], expected[i], epsilon = 0.0);
1058 }
1059 }
1060
1061 macro_rules! gen_public_api_matrix_tests {
1062 ($d:literal) => {
1063 paste! {
1064 #[test]
1065 fn [<public_api_matrix_from_rows_get_set_bounds_checked_ $d d>]() {
1066 let mut rows = [[0.0f64; $d]; $d];
1067 rows[0][0] = 1.0;
1068 rows[$d - 1][$d - 1] = -2.0;
1069
1070 let mut m = Matrix::<$d>::from_rows(rows);
1071
1072 assert_eq!(m.get(0, 0), Some(1.0));
1073 assert_eq!(m.get($d - 1, $d - 1), Some(-2.0));
1074 assert_eq!(m.get_checked(0, 0), Ok(1.0));
1075 assert_eq!(m.get_checked($d - 1, $d - 1), Ok(-2.0));
1076
1077 // Out-of-bounds is None.
1078 assert_eq!(m.get($d, 0), None);
1079 assert_eq!(
1080 m.get_checked($d, 0),
1081 Err(LaError::IndexOutOfBounds {
1082 row: $d,
1083 col: 0,
1084 dim: $d,
1085 })
1086 );
1087
1088 // Out-of-bounds set fails.
1089 let before_failed_set = m;
1090 assert_eq!(
1091 m.set($d, 0, 3.0),
1092 Err(LaError::IndexOutOfBounds {
1093 row: $d,
1094 col: 0,
1095 dim: $d,
1096 })
1097 );
1098 assert_eq!(m, before_failed_set);
1099 assert_eq!(
1100 m.set_checked($d, 0, 3.0),
1101 Err(LaError::IndexOutOfBounds {
1102 row: $d,
1103 col: 0,
1104 dim: $d,
1105 })
1106 );
1107 assert_eq!(m, before_failed_set);
1108 assert_eq!(
1109 m.set_checked(0, $d, 3.0),
1110 Err(LaError::IndexOutOfBounds {
1111 row: 0,
1112 col: $d,
1113 dim: $d,
1114 })
1115 );
1116 assert_eq!(m, before_failed_set);
1117 assert_eq!(m.get(0, 0), Some(1.0));
1118
1119 // In-bounds set works.
1120 assert_eq!(m.set(0, $d - 1, 3.0), Ok(()));
1121 assert_eq!(m.get(0, $d - 1), Some(3.0));
1122 assert_eq!(m.set_checked($d - 1, 0, 4.0), Ok(()));
1123 assert_eq!(m.get_checked($d - 1, 0), Ok(4.0));
1124 }
1125
1126 #[test]
1127 fn [<public_api_matrix_zero_and_default_are_zero_ $d d>]() {
1128 let z = Matrix::<$d>::zero();
1129 assert_abs_diff_eq!(z.inf_norm().unwrap(), 0.0, epsilon = 0.0);
1130
1131 let d = Matrix::<$d>::default();
1132 assert_abs_diff_eq!(d.inf_norm().unwrap(), 0.0, epsilon = 0.0);
1133 }
1134
1135 #[test]
1136 fn [<public_api_matrix_inf_norm_max_row_sum_ $d d>]() {
1137 let mut rows = [[0.0f64; $d]; $d];
1138
1139 // Row 0 has absolute row sum = D.
1140 for c in 0..$d {
1141 rows[0][c] = -1.0;
1142 }
1143
1144 // Row 1 has smaller absolute row sum.
1145 for c in 0..$d {
1146 rows[1][c] = 0.5;
1147 }
1148
1149 let m = Matrix::<$d>::from_rows(rows);
1150 assert_abs_diff_eq!(m.inf_norm().unwrap(), f64::from($d), epsilon = 0.0);
1151 }
1152
1153 #[test]
1154 fn [<public_api_matrix_identity_lu_det_solve_vec_ $d d>]() {
1155 let m = Matrix::<$d>::identity();
1156
1157 // Identity has ones on diag and zeros off diag.
1158 for r in 0..$d {
1159 for c in 0..$d {
1160 let expected = if r == c { 1.0 } else { 0.0 };
1161 assert_abs_diff_eq!(m.get(r, c).unwrap(), expected, epsilon = 0.0);
1162 }
1163 }
1164
1165 // Determinant is 1.
1166 let det = m.det().unwrap();
1167 assert_abs_diff_eq!(det, 1.0, epsilon = 1e-12);
1168
1169 // LU solve on identity returns the RHS.
1170 let lu = m.lu(DEFAULT_PIVOT_TOL).unwrap();
1171
1172 let b_arr = {
1173 let mut arr = [0.0f64; $d];
1174 let values = [1.0f64, 2.0, 3.0, 4.0, 5.0];
1175 for (dst, src) in arr.iter_mut().zip(values.iter()) {
1176 *dst = *src;
1177 }
1178 arr
1179 };
1180
1181 let b = Vector::<$d>::new(b_arr);
1182 let x = lu.solve_vec(b).unwrap().into_array();
1183
1184 for (x_i, b_i) in x.iter().zip(b_arr.iter()) {
1185 assert_abs_diff_eq!(*x_i, *b_i, epsilon = 1e-12);
1186 }
1187 }
1188
1189 #[test]
1190 fn [<finite_matrix_accepts_and_roundtrips_ $d d>]() {
1191 let mut rows = [[0.0f64; $d]; $d];
1192 for r in 0..$d {
1193 rows[r][r] = 1.0;
1194 }
1195
1196 let finite = FiniteMatrix::<$d>::from_rows(rows).unwrap();
1197
1198 assert_matrix_abs_eq(&finite.into_matrix(), &Matrix::<$d>::from_rows(rows));
1199 assert_matrix_abs_eq(&FiniteMatrix::<$d>::zero().into_matrix(), &Matrix::<$d>::zero());
1200 assert_matrix_abs_eq(&FiniteMatrix::<$d>::default().into_matrix(), &Matrix::<$d>::zero());
1201 assert_eq!(finite.into_matrix().get(0, 0), Some(1.0));
1202 assert_eq!(finite.into_matrix().get($d, 0), None);
1203 assert_eq!(
1204 finite.into_matrix().get_checked($d, 0),
1205 Err(LaError::IndexOutOfBounds {
1206 row: $d,
1207 col: 0,
1208 dim: $d,
1209 })
1210 );
1211 assert_matrix_abs_eq(&Matrix::from(finite), &Matrix::<$d>::from_rows(rows));
1212 assert_matrix_abs_eq(
1213 &FiniteMatrix::<$d>::try_from(Matrix::<$d>::from_rows(rows))
1214 .unwrap()
1215 .into_matrix(),
1216 &finite.into_matrix()
1217 );
1218 assert_matrix_abs_eq(
1219 &FiniteMatrix::<$d>::try_from(rows).unwrap().into_matrix(),
1220 &finite.into_matrix()
1221 );
1222 }
1223
1224 #[test]
1225 fn [<finite_matrix_rejects_nonfinite_with_coordinates_ $d d>]() {
1226 let mut rows = [[1.0f64; $d]; $d];
1227 rows[$d - 1][0] = f64::NAN;
1228
1229 assert_eq!(
1230 FiniteMatrix::<$d>::from_rows(rows),
1231 Err(LaError::NonFinite {
1232 row: Some($d - 1),
1233 col: 0,
1234 })
1235 );
1236 }
1237
1238 #[test]
1239 fn [<finite_matrix_try_from_raw_matrix_revalidates_entries_ $d d>]() {
1240 let mut rows = [[1.0f64; $d]; $d];
1241 rows[$d - 1][$d - 1] = f64::INFINITY;
1242 let raw = Matrix::<$d>::from_rows_unchecked(rows);
1243
1244 assert_eq!(
1245 FiniteMatrix::<$d>::try_from(raw),
1246 Err(LaError::NonFinite {
1247 row: Some($d - 1),
1248 col: $d - 1,
1249 })
1250 );
1251 }
1252
1253 #[test]
1254 fn [<finite_matrix_algorithms_match_raw_boundary_ $d d>]() {
1255 let mut rows = [[0.0f64; $d]; $d];
1256 let diag_values = [2.0f64, 3.0, 5.0, 7.0, 11.0];
1257 for i in 0..$d {
1258 rows[i][i] = diag_values[i];
1259 }
1260
1261 let raw = Matrix::<$d>::from_rows(rows);
1262 let finite = FiniteMatrix::<$d>::from_rows(rows).unwrap();
1263
1264 assert_abs_diff_eq!(finite.inf_norm().unwrap(), raw.inf_norm().unwrap(), epsilon = 0.0);
1265 assert_eq!(finite.det_direct(), raw.det_direct());
1266 assert_abs_diff_eq!(finite.det().unwrap(), raw.det().unwrap(), epsilon = 0.0);
1267 assert_eq!(finite.det_errbound(), raw.det_errbound());
1268 assert_eq!(
1269 finite.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(),
1270 raw.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap()
1271 );
1272 assert_eq!(
1273 finite.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap(),
1274 raw.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap()
1275 );
1276 }
1277
1278 #[test]
1279 fn [<finite_matrix_lu_and_ldlt_accept_finite_rhs_ $d d>]() {
1280 let finite = FiniteMatrix::<$d>::new(Matrix::<$d>::identity()).unwrap();
1281 let rhs = {
1282 let mut arr = [0.0f64; $d];
1283 let values = [1.0f64, 2.0, 3.0, 4.0, 5.0];
1284 for (dst, src) in arr.iter_mut().zip(values.iter()) {
1285 *dst = *src;
1286 }
1287 FiniteVector::<$d>::from_array(arr).unwrap()
1288 };
1289
1290 let lu = finite.lu(DEFAULT_PIVOT_TOL).unwrap();
1291 assert_array_abs_eq(
1292 &lu.solve_finite_vec(rhs).unwrap().into_array(),
1293 &rhs.into_array()
1294 );
1295
1296 let ldlt = finite.ldlt(DEFAULT_PIVOT_TOL).unwrap();
1297 assert_array_abs_eq(
1298 &ldlt.solve_finite_vec(rhs).unwrap().into_array(),
1299 &rhs.into_array()
1300 );
1301 }
1302 }
1303 };
1304 }
1305
1306 // Mirror delaunay-style multi-dimension tests.
1307 gen_public_api_matrix_tests!(2);
1308 gen_public_api_matrix_tests!(3);
1309 gen_public_api_matrix_tests!(4);
1310 gen_public_api_matrix_tests!(5);
1311
1312 macro_rules! gen_public_matrix_revalidation_tests {
1313 ($d:literal) => {
1314 paste! {
1315 #[test]
1316 fn [<public_matrix_compute_methods_revalidate_unchecked_storage_ $d d>]() {
1317 let mut rows = [[0.0f64; $d]; $d];
1318 for i in 0..$d {
1319 rows[i][i] = 1.0;
1320 }
1321 rows[0][$d - 1] = f64::NAN;
1322 let raw = Matrix::<$d>::from_rows_unchecked(rows);
1323 let expected = LaError::NonFinite {
1324 row: Some(0),
1325 col: $d - 1,
1326 };
1327
1328 assert_eq!(raw.inf_norm(), Err(expected));
1329 assert_eq!(
1330 raw.is_symmetric(Tolerance::new(0.0).unwrap()),
1331 Err(expected)
1332 );
1333 assert_eq!(
1334 raw.first_asymmetry(Tolerance::new(0.0).unwrap()),
1335 Err(expected)
1336 );
1337 assert_eq!(raw.lu(DEFAULT_PIVOT_TOL), Err(expected));
1338 assert_eq!(raw.ldlt(DEFAULT_PIVOT_TOL), Err(expected));
1339 assert_eq!(raw.det_direct(), Err(expected));
1340 assert_eq!(raw.det(), Err(expected));
1341 assert_eq!(raw.det_errbound(), Err(expected));
1342 }
1343 }
1344 };
1345 }
1346
1347 gen_public_matrix_revalidation_tests!(2);
1348 gen_public_matrix_revalidation_tests!(3);
1349 gen_public_matrix_revalidation_tests!(4);
1350 gen_public_matrix_revalidation_tests!(5);
1351
1352 #[test]
1353 fn public_matrix_fast_none_paths_revalidate_unchecked_storage() {
1354 let mut rows = [[1.0f64; 5]; 5];
1355 rows[4][1] = f64::INFINITY;
1356 let raw = Matrix::<5>::from_rows_unchecked(rows);
1357
1358 let expected = Err(LaError::NonFinite {
1359 row: Some(4),
1360 col: 1,
1361 });
1362
1363 assert_eq!(raw.det_direct(), expected);
1364 assert_eq!(raw.det_errbound(), expected);
1365 }
1366
1367 // === det_direct tests ===
1368
1369 #[test]
1370 fn det_direct_d0_is_one() {
1371 assert_eq!(Matrix::<0>::zero().det_direct(), Ok(Some(1.0)));
1372 }
1373
1374 #[test]
1375 fn det_direct_d1_returns_element() {
1376 let m = Matrix::<1>::from_rows([[42.0]]);
1377 assert_eq!(m.det_direct(), Ok(Some(42.0)));
1378 }
1379
1380 #[test]
1381 fn det_direct_d2_known_value() {
1382 // [[1,2],[3,4]] → det = 1*4 - 2*3 = -2
1383 // black_box prevents compile-time constant folding of the const fn.
1384 let m = black_box(Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]));
1385 assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), -2.0, epsilon = 1e-15);
1386 }
1387
1388 #[test]
1389 fn det_direct_d3_known_value() {
1390 // Classic 3×3: det = 0
1391 let m = black_box(Matrix::<3>::from_rows([
1392 [1.0, 2.0, 3.0],
1393 [4.0, 5.0, 6.0],
1394 [7.0, 8.0, 9.0],
1395 ]));
1396 assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), 0.0, epsilon = 1e-12);
1397 }
1398
1399 #[test]
1400 fn det_direct_d3_nonsingular() {
1401 // [[2,1,0],[0,3,1],[1,0,2]] → det = 2*(6-0) - 1*(0-1) + 0 = 13
1402 let m = black_box(Matrix::<3>::from_rows([
1403 [2.0, 1.0, 0.0],
1404 [0.0, 3.0, 1.0],
1405 [1.0, 0.0, 2.0],
1406 ]));
1407 assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), 13.0, epsilon = 1e-12);
1408 }
1409
1410 #[test]
1411 fn det_direct_d4_identity() {
1412 let m = black_box(Matrix::<4>::identity());
1413 assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), 1.0, epsilon = 1e-15);
1414 }
1415
1416 #[test]
1417 fn det_direct_d4_known_value() {
1418 // Diagonal matrix: det = product of diagonal entries.
1419 let mut rows = [[0.0f64; 4]; 4];
1420 rows[0][0] = 2.0;
1421 rows[1][1] = 3.0;
1422 rows[2][2] = 5.0;
1423 rows[3][3] = 7.0;
1424 let m = black_box(Matrix::<4>::from_rows(rows));
1425 assert_abs_diff_eq!(m.det_direct().unwrap().unwrap(), 210.0, epsilon = 1e-12);
1426 }
1427
1428 #[test]
1429 fn det_direct_d5_returns_none() {
1430 assert_eq!(Matrix::<5>::identity().det_direct(), Ok(None));
1431 }
1432
1433 #[test]
1434 fn det_direct_d5_rejects_nonfinite_before_returning_none() {
1435 let mut m = Matrix::<5>::identity();
1436 assert_eq!(
1437 m.set(3, 4, f64::NAN),
1438 Err(LaError::NonFinite {
1439 row: Some(3),
1440 col: 4,
1441 })
1442 );
1443 }
1444
1445 #[test]
1446 fn det_direct_d8_returns_none() {
1447 assert_eq!(Matrix::<8>::zero().det_direct(), Ok(None));
1448 }
1449
1450 #[test]
1451 fn det_direct_rejects_nonfinite_entry_with_coordinates() {
1452 assert_eq!(
1453 Matrix::<3>::try_from_rows([[1.0, 0.0, 0.0], [0.0, f64::NAN, 0.0], [0.0, 0.0, 1.0]]),
1454 Err(LaError::NonFinite {
1455 row: Some(1),
1456 col: 1,
1457 })
1458 );
1459 }
1460
1461 #[test]
1462 fn det_direct_rejects_computed_overflow() {
1463 let m = Matrix::<2>::from_rows([[1e300, 0.0], [0.0, 1e300]]);
1464 assert_eq!(
1465 m.det_direct(),
1466 Err(LaError::NonFinite { row: None, col: 0 })
1467 );
1468 }
1469
1470 #[test]
1471 fn det_d5_rejects_lu_product_overflow() {
1472 let m = Matrix::<5>::from_rows([
1473 [1.0e100, 0.0, 0.0, 0.0, 0.0],
1474 [0.0, 1.0e100, 0.0, 0.0, 0.0],
1475 [0.0, 0.0, 1.0e100, 0.0, 0.0],
1476 [0.0, 0.0, 0.0, 1.0e100, 0.0],
1477 [0.0, 0.0, 0.0, 0.0, 1.0e100],
1478 ]);
1479 assert_eq!(m.det(), Err(LaError::NonFinite { row: None, col: 3 }));
1480 }
1481
1482 #[test]
1483 fn det_d5_rejects_lu_trailing_update_overflow() {
1484 let m = Matrix::<5>::from_rows([
1485 [1.0, f64::MAX, 0.0, 0.0, 0.0],
1486 [-1.0, f64::MAX, 0.0, 0.0, 0.0],
1487 [0.0, 0.0, 1.0, 0.0, 0.0],
1488 [0.0, 0.0, 0.0, 1.0, 0.0],
1489 [0.0, 0.0, 0.0, 0.0, 1.0],
1490 ]);
1491
1492 assert_eq!(
1493 m.det(),
1494 Err(LaError::NonFinite {
1495 row: Some(1),
1496 col: 1,
1497 })
1498 );
1499 }
1500
1501 macro_rules! gen_det_direct_agrees_with_lu {
1502 ($d:literal) => {
1503 paste! {
1504 #[test]
1505 #[allow(clippy::cast_precision_loss)] // r, c, D are tiny integers
1506 fn [<det_direct_agrees_with_lu_ $d d>]() {
1507 // Well-conditioned matrix: diagonally dominant.
1508 let mut rows = [[0.0f64; $d]; $d];
1509 for r in 0..$d {
1510 for c in 0..$d {
1511 rows[r][c] = if r == c {
1512 (r as f64) + f64::from($d) + 1.0
1513 } else {
1514 0.1 / ((r + c + 1) as f64)
1515 };
1516 }
1517 }
1518 let m = Matrix::<$d>::from_rows(rows);
1519 let direct = m.det_direct().unwrap().unwrap();
1520 let lu_det = m.lu(DEFAULT_PIVOT_TOL).unwrap().det().unwrap();
1521 let eps = lu_det.abs().mul_add(1e-12, 1e-12);
1522 assert_abs_diff_eq!(direct, lu_det, epsilon = eps);
1523 }
1524 }
1525 };
1526 }
1527
1528 gen_det_direct_agrees_with_lu!(1);
1529 gen_det_direct_agrees_with_lu!(2);
1530 gen_det_direct_agrees_with_lu!(3);
1531 gen_det_direct_agrees_with_lu!(4);
1532
1533 #[test]
1534 fn det_direct_identity_all_dims() {
1535 assert_abs_diff_eq!(
1536 Matrix::<1>::identity().det_direct().unwrap().unwrap(),
1537 1.0,
1538 epsilon = 0.0
1539 );
1540 assert_abs_diff_eq!(
1541 Matrix::<2>::identity().det_direct().unwrap().unwrap(),
1542 1.0,
1543 epsilon = 0.0
1544 );
1545 assert_abs_diff_eq!(
1546 Matrix::<3>::identity().det_direct().unwrap().unwrap(),
1547 1.0,
1548 epsilon = 0.0
1549 );
1550 assert_abs_diff_eq!(
1551 Matrix::<4>::identity().det_direct().unwrap().unwrap(),
1552 1.0,
1553 epsilon = 0.0
1554 );
1555 }
1556
1557 #[test]
1558 fn det_direct_zero_matrix() {
1559 assert_abs_diff_eq!(
1560 Matrix::<2>::zero().det_direct().unwrap().unwrap(),
1561 0.0,
1562 epsilon = 0.0
1563 );
1564 assert_abs_diff_eq!(
1565 Matrix::<3>::zero().det_direct().unwrap().unwrap(),
1566 0.0,
1567 epsilon = 0.0
1568 );
1569 assert_abs_diff_eq!(
1570 Matrix::<4>::zero().det_direct().unwrap().unwrap(),
1571 0.0,
1572 epsilon = 0.0
1573 );
1574 }
1575
1576 macro_rules! gen_det_singular_zero_matrix_tests {
1577 ($d:literal) => {
1578 paste! {
1579 #[test]
1580 fn [<det_singular_zero_matrix_returns_zero_ $d d>]() {
1581 assert_abs_diff_eq!(
1582 Matrix::<$d>::zero().det().unwrap(),
1583 0.0,
1584 epsilon = 0.0
1585 );
1586 }
1587 }
1588 };
1589 }
1590
1591 gen_det_singular_zero_matrix_tests!(2);
1592 gen_det_singular_zero_matrix_tests!(3);
1593 gen_det_singular_zero_matrix_tests!(4);
1594 gen_det_singular_zero_matrix_tests!(5);
1595
1596 #[test]
1597 fn det_d5_ignores_pivot_tolerance_for_tiny_nonsingular_matrix() {
1598 // A small nonzero determinant is still a determinant. `det` must not
1599 // flatten the value to zero merely because the default LU tolerance
1600 // would reject a pivot this small.
1601 let m = Matrix::<5>::from_rows([
1602 [1e-13, 0.0, 0.0, 0.0, 0.0],
1603 [0.0, 1.0, 0.0, 0.0, 0.0],
1604 [0.0, 0.0, 1.0, 0.0, 0.0],
1605 [0.0, 0.0, 0.0, 1.0, 0.0],
1606 [0.0, 0.0, 0.0, 0.0, 1.0],
1607 ]);
1608
1609 assert_abs_diff_eq!(m.det().unwrap(), 1e-13, epsilon = 0.0);
1610 assert_eq!(
1611 m.lu(DEFAULT_PIVOT_TOL),
1612 Err(LaError::Singular { pivot_col: 0 })
1613 );
1614 }
1615
1616 #[test]
1617 fn det_returns_nonfinite_error_for_nan_d2() {
1618 assert_eq!(
1619 Matrix::<2>::try_from_rows([[f64::NAN, 1.0], [1.0, 1.0]]),
1620 Err(LaError::NonFinite {
1621 row: Some(0),
1622 col: 0
1623 })
1624 );
1625 }
1626
1627 #[test]
1628 fn det_returns_nonfinite_error_for_inf_d3() {
1629 assert_eq!(
1630 Matrix::<3>::try_from_rows([
1631 [f64::INFINITY, 0.0, 0.0],
1632 [0.0, 1.0, 0.0],
1633 [0.0, 0.0, 1.0]
1634 ]),
1635 Err(LaError::NonFinite {
1636 row: Some(0),
1637 col: 0
1638 })
1639 );
1640 }
1641
1642 #[test]
1643 fn det_returns_nonfinite_error_for_overflow_with_finite_entries() {
1644 // det_direct produces an overflowing f64 (1e300 * 1e300 = ∞) even
1645 // though every matrix entry is finite. The entry scan in `det`
1646 // falls through and returns NonFinite { row: None, col: 0 } to signal
1647 // a computed overflow rather than a NaN/∞ input.
1648 let m = Matrix::<2>::from_rows([[1e300, 0.0], [0.0, 1e300]]);
1649 assert_eq!(m.det(), Err(LaError::NonFinite { row: None, col: 0 }));
1650 }
1651
1652 // === det_direct const-evaluability tests (D = 2..=5) ===
1653 //
1654 // Every dimension hits a distinct arm of the `match D { … }` body inside
1655 // `det_direct`, so exercising each at compile time is the tightest
1656 // const-fn proof available.
1657
1658 macro_rules! gen_det_direct_const_eval_tests {
1659 ($d:literal) => {
1660 paste! {
1661 /// `Matrix::<D>::det_direct()` on the identity must const-evaluate
1662 /// to `Ok(Some(1.0))` for every closed-form dimension `D ∈ {1, 2, 3, 4}`.
1663 #[test]
1664 fn [<det_direct_const_eval_ $d d>]() {
1665 const DET: Result<Option<f64>, LaError> = Matrix::<$d>::identity().det_direct();
1666 assert_eq!(DET, Ok(Some(1.0)));
1667 }
1668 }
1669 };
1670 }
1671
1672 gen_det_direct_const_eval_tests!(2);
1673 gen_det_direct_const_eval_tests!(3);
1674 gen_det_direct_const_eval_tests!(4);
1675
1676 #[test]
1677 fn det_direct_const_eval_d5_is_none() {
1678 // D ≥ 5 has no closed-form arm; `det_direct` returns `Ok(None)`. Verify
1679 // that the wildcard arm is reachable in a `const { … }` context.
1680 const DET: Result<Option<f64>, LaError> = Matrix::<5>::identity().det_direct();
1681 assert_eq!(DET, Ok(None));
1682 }
1683
1684 // === det_errbound tests (no `exact` feature required) ===
1685
1686 #[test]
1687 fn det_errbound_available_without_exact_feature() {
1688 // Verify det_errbound is accessible without exact feature
1689 let m = Matrix::<3>::identity();
1690 let bound = m.det_errbound().unwrap();
1691 assert!(bound.is_some());
1692 assert!(bound.unwrap() > 0.0);
1693 }
1694
1695 #[test]
1696 fn det_errbound_matches_documented_coefficient_scale() {
1697 let m2 = Matrix::<2>::from_rows([[1.0, 2.0], [3.0, 4.0]]);
1698 let expected_2 = ERR_COEFF_2 * ((1.0_f64 * 4.0).abs() + (2.0_f64 * 3.0).abs());
1699 assert_abs_diff_eq!(
1700 m2.det_errbound().unwrap().unwrap(),
1701 expected_2,
1702 epsilon = 0.0
1703 );
1704
1705 assert_abs_diff_eq!(
1706 Matrix::<3>::identity().det_errbound().unwrap().unwrap(),
1707 ERR_COEFF_3,
1708 epsilon = 0.0
1709 );
1710 assert_abs_diff_eq!(
1711 Matrix::<4>::identity().det_errbound().unwrap().unwrap(),
1712 ERR_COEFF_4,
1713 epsilon = 0.0
1714 );
1715 }
1716
1717 #[test]
1718 fn det_errbound_d5_returns_none() {
1719 // D=5 has no fast filter
1720 assert_eq!(Matrix::<5>::identity().det_errbound(), Ok(None));
1721 }
1722
1723 #[test]
1724 fn det_errbound_d1_rejects_nonfinite_even_with_zero_bound() {
1725 assert_eq!(
1726 Matrix::<1>::try_from_rows([[f64::INFINITY]]),
1727 Err(LaError::NonFinite {
1728 row: Some(0),
1729 col: 0,
1730 })
1731 );
1732 }
1733
1734 #[test]
1735 fn det_errbound_d5_rejects_nonfinite_before_returning_none() {
1736 let mut m = Matrix::<5>::identity();
1737 assert_eq!(
1738 m.set(4, 1, f64::NAN),
1739 Err(LaError::NonFinite {
1740 row: Some(4),
1741 col: 1,
1742 })
1743 );
1744 }
1745
1746 #[test]
1747 fn det_errbound_rejects_nonfinite_entry_with_coordinates() {
1748 assert_eq!(
1749 Matrix::<2>::try_from_rows([[1.0, f64::INFINITY], [0.0, 1.0]]),
1750 Err(LaError::NonFinite {
1751 row: Some(0),
1752 col: 1,
1753 })
1754 );
1755 }
1756
1757 #[test]
1758 fn det_errbound_rejects_computed_overflow() {
1759 let m = Matrix::<2>::from_rows([[1e300, 0.0], [0.0, 1e300]]);
1760 assert_eq!(
1761 m.det_errbound(),
1762 Err(LaError::NonFinite { row: None, col: 0 })
1763 );
1764 }
1765
1766 // === det_errbound const-evaluability tests (D = 2..=5) ===
1767
1768 macro_rules! gen_det_errbound_const_eval_tests {
1769 ($d:literal) => {
1770 paste! {
1771 /// `Matrix::<D>::det_errbound()` on the identity must const-evaluate
1772 /// to `Ok(Some(bound))` with `bound > 0` for every closed-form dimension
1773 /// `D ∈ {2, 3, 4}`. Each dimension hits a distinct arm of
1774 /// `det_errbound` with a dimension-specific permanent computation.
1775 #[test]
1776 fn [<det_errbound_const_eval_ $d d>]() {
1777 const BOUND: Result<Option<f64>, LaError> = Matrix::<$d>::identity().det_errbound();
1778 assert!(BOUND.unwrap().unwrap() > 0.0);
1779 }
1780 }
1781 };
1782 }
1783
1784 gen_det_errbound_const_eval_tests!(2);
1785 gen_det_errbound_const_eval_tests!(3);
1786 gen_det_errbound_const_eval_tests!(4);
1787
1788 #[test]
1789 fn det_errbound_const_eval_d5_is_none() {
1790 // D ≥ 5 has no fast-filter bound; `det_errbound` returns `Ok(None)`.
1791 const BOUND: Result<Option<f64>, LaError> = Matrix::<5>::identity().det_errbound();
1792 assert_eq!(BOUND, Ok(None));
1793 }
1794
1795 // === inf_norm const-evaluability tests (D = 2..=5) ===
1796
1797 macro_rules! gen_inf_norm_const_eval_tests {
1798 ($d:literal) => {
1799 paste! {
1800 /// `Matrix::<D>::inf_norm()` on the identity must const-evaluate
1801 /// to `1.0` for every `D ≥ 1` — each row has a single `1.0`
1802 /// entry, so the max absolute row sum is exactly `1.0`.
1803 #[test]
1804 fn [<inf_norm_const_eval_ $d d>]() {
1805 const NORM: Result<f64, LaError> = Matrix::<$d>::identity().inf_norm();
1806 assert!((NORM.unwrap() - 1.0).abs() <= 1e-12);
1807 }
1808 }
1809 };
1810 }
1811
1812 gen_inf_norm_const_eval_tests!(2);
1813 gen_inf_norm_const_eval_tests!(3);
1814 gen_inf_norm_const_eval_tests!(4);
1815 gen_inf_norm_const_eval_tests!(5);
1816
1817 // === inf_norm NaN / Inf rejection (regression tests for #85) ===
1818
1819 macro_rules! gen_inf_norm_nonfinite_tests {
1820 ($d:literal) => {
1821 paste! {
1822 #[test]
1823 fn [<inf_norm_all_nan_returns_nonfinite_error_ $d d>]() {
1824 // Before the fix, `NaN > max_row_sum` was always false, so a
1825 // matrix full of NaN silently produced inf_norm == 0.0.
1826 assert_eq!(
1827 Matrix::<$d>::try_from_rows([[f64::NAN; $d]; $d]),
1828 Err(LaError::NonFinite {
1829 row: Some(0),
1830 col: 0,
1831 })
1832 );
1833 }
1834
1835 #[test]
1836 fn [<inf_norm_single_nan_entry_returns_nonfinite_error_ $d d>]() {
1837 // A single NaN entry must surface with its source coordinates.
1838 let mut rows = [[0.0f64; $d]; $d];
1839 rows[0][0] = f64::NAN;
1840 rows[$d - 1][$d - 1] = 1.0;
1841 assert_eq!(
1842 Matrix::<$d>::try_from_rows(rows),
1843 Err(LaError::NonFinite {
1844 row: Some(0),
1845 col: 0,
1846 })
1847 );
1848 }
1849
1850 #[test]
1851 fn [<inf_norm_infinity_entry_returns_nonfinite_error_ $d d>]() {
1852 // Infinity entries should be rejected with their source coordinates.
1853 let mut rows = [[0.0f64; $d]; $d];
1854 rows[0][0] = f64::INFINITY;
1855 assert_eq!(
1856 Matrix::<$d>::try_from_rows(rows),
1857 Err(LaError::NonFinite {
1858 row: Some(0),
1859 col: 0,
1860 })
1861 );
1862 }
1863 }
1864 };
1865 }
1866
1867 gen_inf_norm_nonfinite_tests!(2);
1868 gen_inf_norm_nonfinite_tests!(3);
1869 gen_inf_norm_nonfinite_tests!(4);
1870 gen_inf_norm_nonfinite_tests!(5);
1871
1872 // === is_symmetric / first_asymmetry (public LDLT preconditions helpers) ===
1873
1874 macro_rules! gen_is_symmetric_tests {
1875 ($d:literal) => {
1876 paste! {
1877 #[test]
1878 fn [<is_symmetric_true_for_identity_ $d d>]() {
1879 let m = Matrix::<$d>::identity();
1880 assert!(m.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap());
1881 assert_eq!(m.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(), None);
1882 }
1883
1884 #[test]
1885 fn [<is_symmetric_true_for_zero_ $d d>]() {
1886 let m = Matrix::<$d>::zero();
1887 assert!(m.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap());
1888 assert_eq!(m.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(), None);
1889 }
1890
1891 #[test]
1892 fn [<is_symmetric_true_for_constructed_symmetric_ $d d>]() {
1893 // Construct A = M + Mᵀ so A is provably symmetric.
1894 let mut m = [[0.0f64; $d]; $d];
1895 for r in 0..$d {
1896 for c in 0..$d {
1897 #[allow(clippy::cast_precision_loss)]
1898 {
1899 m[r][c] = (r * $d + c) as f64;
1900 }
1901 }
1902 }
1903 let mut sym = [[0.0f64; $d]; $d];
1904 for r in 0..$d {
1905 for c in 0..$d {
1906 sym[r][c] = m[r][c] + m[c][r];
1907 }
1908 }
1909 let a = Matrix::<$d>::from_rows(sym);
1910 assert!(a.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap());
1911 assert_eq!(a.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(), None);
1912 }
1913
1914 #[test]
1915 fn [<is_symmetric_false_for_asymmetric_offdiagonal_ $d d>]() {
1916 // Perturb a single off-diagonal entry so symmetry fails.
1917 let mut rows = [[0.0f64; $d]; $d];
1918 for i in 0..$d {
1919 rows[i][i] = 1.0;
1920 }
1921 rows[0][$d - 1] = 1.0;
1922 rows[$d - 1][0] = -1.0; // breaks symmetry
1923 let a = Matrix::<$d>::from_rows(rows);
1924 assert!(!a.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap());
1925 assert_eq!(
1926 a.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(),
1927 Some((0, $d - 1))
1928 );
1929 }
1930
1931 #[test]
1932 fn [<is_symmetric_rejects_nan_offdiagonal_ $d d>]() {
1933 // A NaN off-diagonal is a stored non-finite matrix value,
1934 // not merely a symmetry mismatch.
1935 let mut rows = [[0.0f64; $d]; $d];
1936 for i in 0..$d {
1937 rows[i][i] = 1.0;
1938 }
1939 rows[0][1] = f64::NAN;
1940 rows[1][0] = f64::NAN;
1941 assert_eq!(
1942 Matrix::<$d>::try_from_rows(rows),
1943 Err(LaError::NonFinite {
1944 row: Some(0),
1945 col: 1,
1946 })
1947 );
1948 }
1949 }
1950 };
1951 }
1952
1953 gen_is_symmetric_tests!(2);
1954 gen_is_symmetric_tests!(3);
1955 gen_is_symmetric_tests!(4);
1956 gen_is_symmetric_tests!(5);
1957
1958 macro_rules! gen_ldlt_symmetry_proof_tests {
1959 ($d:literal) => {
1960 paste! {
1961 #[test]
1962 fn [<matrix_ldlt_accepts_identity_ $d d>]() {
1963 let ldlt = Matrix::<$d>::identity().ldlt(DEFAULT_PIVOT_TOL).unwrap();
1964 assert_abs_diff_eq!(ldlt.det().unwrap(), 1.0, epsilon = 0.0);
1965 }
1966
1967 #[test]
1968 fn [<symmetric_matrix_try_new_rejects_finite_asymmetric_ $d d>]() {
1969 let mut rows = [[0.0f64; $d]; $d];
1970 for (i, row) in rows.iter_mut().enumerate() {
1971 row[i] = 1.0;
1972 }
1973 rows[0][$d - 1] = 1.0;
1974 rows[$d - 1][0] = -1.0;
1975
1976 assert_eq!(
1977 Matrix::<$d>::try_from_rows(rows)
1978 .and_then(FiniteMatrix::new)
1979 .and_then(SymmetricMatrix::try_new),
1980 Err(LaError::Asymmetric {
1981 row: 0,
1982 col: $d - 1,
1983 dim: $d,
1984 })
1985 );
1986 }
1987 }
1988 };
1989 }
1990
1991 gen_ldlt_symmetry_proof_tests!(2);
1992 gen_ldlt_symmetry_proof_tests!(3);
1993 gen_ldlt_symmetry_proof_tests!(4);
1994 gen_ldlt_symmetry_proof_tests!(5);
1995
1996 #[test]
1997 fn matrix_ldlt_rejects_nonfinite_before_asymmetry() {
1998 assert_eq!(
1999 Matrix::<2>::try_from_rows([[1.0, f64::NAN], [0.0, 1.0]]),
2000 Err(LaError::NonFinite {
2001 row: Some(0),
2002 col: 1,
2003 })
2004 );
2005 }
2006
2007 #[test]
2008 fn symmetric_matrix_into_matrix_roundtrips_storage_internally() {
2009 let a = Matrix::<2>::from_rows([[2.0, 1.0], [1.0, 3.0]]);
2010 let symmetric = SymmetricMatrix::try_new(FiniteMatrix::new(a).unwrap()).unwrap();
2011
2012 assert_eq!(symmetric.into_matrix(), a);
2013 }
2014
2015 #[test]
2016 fn is_symmetric_tolerance_scales_with_inf_norm() {
2017 // Off-diagonal entries differ by 1e-6. With inf_norm ≈ 2e6, the
2018 // relative tolerance 1e-12 yields eps ≈ 2e-6, which accepts the gap;
2019 // a stricter tol of 1e-15 rejects it.
2020 let a = Matrix::<2>::from_rows([[1.0e6, 1.0e6 + 1.0e-6], [1.0e6, 1.0e6]]);
2021 assert!(a.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap());
2022 assert!(!a.is_symmetric(Tolerance::new(1e-15).unwrap()).unwrap());
2023 }
2024
2025 #[test]
2026 fn first_asymmetry_returns_lexicographically_first_pair() {
2027 // Two asymmetric pairs: (0, 2) and (1, 2). We must get (0, 2) first.
2028 let a = Matrix::<3>::from_rows([[1.0, 0.0, 2.0], [0.0, 1.0, 3.0], [-2.0, -3.0, 1.0]]);
2029 assert_eq!(
2030 a.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(),
2031 Some((0, 2))
2032 );
2033 }
2034
2035 #[test]
2036 fn first_asymmetry_rejects_infinite_offdiagonal() {
2037 assert_eq!(
2038 Matrix::<2>::try_from_rows([[1.0, f64::INFINITY], [0.0, 1.0]]),
2039 Err(LaError::NonFinite {
2040 row: Some(0),
2041 col: 1,
2042 })
2043 );
2044 }
2045
2046 #[test]
2047 fn first_asymmetry_rejects_nan_diagonal() {
2048 assert_eq!(
2049 Matrix::<2>::try_from_rows([[f64::NAN, 1.0], [1.0, 1.0]]),
2050 Err(LaError::NonFinite {
2051 row: Some(0),
2052 col: 0,
2053 })
2054 );
2055 }
2056
2057 #[test]
2058 fn first_asymmetry_strict_tol_survives_row_sum_overflow() {
2059 let a = Matrix::<3>::from_rows([
2060 [f64::MAX, 1.0, f64::MAX],
2061 [2.0, f64::MAX, 0.0],
2062 [f64::MAX, 0.0, f64::MAX],
2063 ]);
2064
2065 assert_eq!(a.inf_norm(), Err(LaError::NonFinite { row: None, col: 0 }));
2066 assert_eq!(
2067 a.first_asymmetry(Tolerance::new(0.0).unwrap()).unwrap(),
2068 Some((0, 1))
2069 );
2070 assert!(!a.is_symmetric(Tolerance::new(0.0).unwrap()).unwrap());
2071 }
2072
2073 #[test]
2074 fn first_asymmetry_rejects_scaled_epsilon_overflow() {
2075 let a = Matrix::<2>::from_rows([[2.0, 0.0], [0.0, 1.0]]);
2076 let tol = Tolerance::new(f64::MAX).unwrap();
2077
2078 assert_eq!(
2079 a.first_asymmetry(tol),
2080 Err(LaError::NonFinite { row: None, col: 0 })
2081 );
2082 assert_eq!(
2083 a.is_symmetric(tol),
2084 Err(LaError::NonFinite { row: None, col: 0 })
2085 );
2086 }
2087
2088 #[test]
2089 fn first_asymmetry_flags_overflowed_finite_difference() {
2090 let a = Matrix::<2>::from_rows([[1.0, f64::MAX], [-f64::MAX, 1.0]]);
2091 assert_eq!(
2092 a.first_asymmetry(Tolerance::new(1e-12).unwrap()).unwrap(),
2093 Some((0, 1))
2094 );
2095 assert!(!a.is_symmetric(Tolerance::new(1e-12).unwrap()).unwrap());
2096 }
2097
2098 #[test]
2099 fn is_symmetric_rejects_invalid_tol() {
2100 assert_eq!(
2101 Tolerance::new(-1.0),
2102 Err(LaError::InvalidTolerance { value: -1.0 })
2103 );
2104 assert_matches!(
2105 Tolerance::new(f64::NAN),
2106 Err(LaError::InvalidTolerance { value }) if value.is_nan()
2107 );
2108 assert_eq!(
2109 Tolerance::new(f64::INFINITY),
2110 Err(LaError::InvalidTolerance {
2111 value: f64::INFINITY,
2112 })
2113 );
2114 }
2115
2116 #[test]
2117 fn first_asymmetry_rejects_negative_tol() {
2118 assert_eq!(
2119 Tolerance::new(-1.0),
2120 Err(LaError::InvalidTolerance { value: -1.0 })
2121 );
2122 }
2123
2124 #[test]
2125 fn first_asymmetry_rejects_nonfinite_tol() {
2126 assert_matches!(
2127 Tolerance::new(f64::NAN),
2128 Err(LaError::InvalidTolerance { value }) if value.is_nan()
2129 );
2130 assert_eq!(
2131 Tolerance::new(f64::INFINITY),
2132 Err(LaError::InvalidTolerance {
2133 value: f64::INFINITY,
2134 })
2135 );
2136 }
2137}