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