la_stack/
lib.rs

1#![forbid(unsafe_code)]
2#![warn(missing_docs)]
3#![doc = include_str!("../README.md")]
4
5#[cfg(doc)]
6mod readme_doctests {
7    //! Executable versions of README examples.
8    /// ```rust
9    /// use la_stack::prelude::*;
10    ///
11    /// // This system requires pivoting (a[0][0] = 0), so it's a good LU demo.
12    /// let a = Matrix::<5>::from_rows([
13    ///     [0.0, 1.0, 1.0, 1.0, 1.0],
14    ///     [1.0, 0.0, 1.0, 1.0, 1.0],
15    ///     [1.0, 1.0, 0.0, 1.0, 1.0],
16    ///     [1.0, 1.0, 1.0, 0.0, 1.0],
17    ///     [1.0, 1.0, 1.0, 1.0, 0.0],
18    /// ]);
19    ///
20    /// let b = Vector::<5>::new([14.0, 13.0, 12.0, 11.0, 10.0]);
21    ///
22    /// let lu = a.lu(DEFAULT_PIVOT_TOL).unwrap();
23    /// let x = lu.solve_vec(b).unwrap().into_array();
24    ///
25    /// // Floating-point rounding is expected; compare with a tolerance.
26    /// let expected = [1.0, 2.0, 3.0, 4.0, 5.0];
27    /// for (x_i, e_i) in x.iter().zip(expected.iter()) {
28    ///     assert!((*x_i - *e_i).abs() <= 1e-12);
29    /// }
30    /// ```
31    fn solve_5x5_example() {}
32
33    /// ```rust
34    /// use la_stack::prelude::*;
35    ///
36    /// // This matrix is symmetric positive-definite (A = L*L^T) so LDLT works without pivoting.
37    /// let a = Matrix::<5>::from_rows([
38    ///     [1.0, 1.0, 0.0, 0.0, 0.0],
39    ///     [1.0, 2.0, 1.0, 0.0, 0.0],
40    ///     [0.0, 1.0, 2.0, 1.0, 0.0],
41    ///     [0.0, 0.0, 1.0, 2.0, 1.0],
42    ///     [0.0, 0.0, 0.0, 1.0, 2.0],
43    /// ]);
44    ///
45    /// let det = a.ldlt(DEFAULT_SINGULAR_TOL).unwrap().det();
46    /// assert!((det - 1.0).abs() <= 1e-12);
47    /// ```
48    fn det_5x5_ldlt_example() {}
49}
50
51#[cfg(feature = "exact")]
52mod exact;
53mod ldlt;
54mod lu;
55mod matrix;
56mod vector;
57
58use core::fmt;
59
60/// Default absolute threshold used for singularity/degeneracy detection.
61///
62/// This is intentionally conservative for geometric predicates and small systems.
63///
64/// Conceptually, this is an absolute bound for deciding when a scalar should be treated
65/// as "numerically zero" (e.g. LU pivots, LDLT diagonal entries).
66pub const DEFAULT_SINGULAR_TOL: f64 = 1e-12;
67
68/// Default absolute pivot magnitude threshold used for LU pivot selection / singularity detection.
69///
70/// This name is kept for backwards compatibility; prefer [`DEFAULT_SINGULAR_TOL`] when the
71/// tolerance is not specifically about pivot selection.
72pub const DEFAULT_PIVOT_TOL: f64 = DEFAULT_SINGULAR_TOL;
73
74/// Linear algebra errors.
75#[derive(Clone, Copy, Debug, PartialEq, Eq)]
76pub enum LaError {
77    /// The matrix is (numerically) singular.
78    Singular {
79        /// The factorization column/step where a suitable pivot/diagonal could not be found.
80        pivot_col: usize,
81    },
82    /// A non-finite value (NaN/∞) was encountered.
83    NonFinite {
84        /// The column being processed when non-finite values were detected.
85        pivot_col: usize,
86    },
87}
88
89impl fmt::Display for LaError {
90    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
91        match *self {
92            Self::Singular { pivot_col } => {
93                write!(f, "singular matrix at pivot column {pivot_col}")
94            }
95            Self::NonFinite { pivot_col } => {
96                write!(
97                    f,
98                    "non-finite value encountered at pivot column {pivot_col}"
99                )
100            }
101        }
102    }
103}
104
105impl std::error::Error for LaError {}
106
107pub use ldlt::Ldlt;
108pub use lu::Lu;
109pub use matrix::Matrix;
110pub use vector::Vector;
111
112/// Common imports for ergonomic usage.
113///
114/// This prelude re-exports the primary types and constants: [`Matrix`], [`Vector`], [`Lu`],
115/// [`Ldlt`], [`LaError`], [`DEFAULT_PIVOT_TOL`], and [`DEFAULT_SINGULAR_TOL`].
116pub mod prelude {
117    pub use crate::{DEFAULT_PIVOT_TOL, DEFAULT_SINGULAR_TOL, LaError, Ldlt, Lu, Matrix, Vector};
118}
119
120#[cfg(test)]
121mod tests {
122    use super::*;
123
124    use approx::assert_abs_diff_eq;
125
126    #[test]
127    fn default_singular_tol_is_expected() {
128        assert_abs_diff_eq!(DEFAULT_SINGULAR_TOL, 1e-12, epsilon = 0.0);
129        assert_abs_diff_eq!(DEFAULT_PIVOT_TOL, DEFAULT_SINGULAR_TOL, epsilon = 0.0);
130    }
131
132    #[test]
133    fn laerror_display_formats_singular() {
134        let err = LaError::Singular { pivot_col: 3 };
135        assert_eq!(err.to_string(), "singular matrix at pivot column 3");
136    }
137
138    #[test]
139    fn laerror_display_formats_nonfinite() {
140        let err = LaError::NonFinite { pivot_col: 2 };
141        assert_eq!(
142            err.to_string(),
143            "non-finite value encountered at pivot column 2"
144        );
145    }
146
147    #[test]
148    fn laerror_is_std_error_with_no_source() {
149        let err = LaError::Singular { pivot_col: 0 };
150        let e: &dyn std::error::Error = &err;
151        assert!(e.source().is_none());
152    }
153
154    #[test]
155    fn prelude_reexports_compile_and_work() {
156        use crate::prelude::*;
157
158        // Use the items so we know they are in scope and usable.
159        let m = Matrix::<2>::identity();
160        let v = Vector::<2>::new([1.0, 2.0]);
161        let _ = m.lu(DEFAULT_PIVOT_TOL).unwrap().solve_vec(v).unwrap();
162        let _ = m.ldlt(DEFAULT_SINGULAR_TOL).unwrap().solve_vec(v).unwrap();
163    }
164}
la_stack/lib.rs

la_stack/
lib.rs
