uncertain_numerics/
linear_system.rs

1use nalgebra::DMatrix;
2
3use crate::LinearSolverError;
4
5const SYMMETRY_TOLERANCE: f64 = 64.0 * f64::EPSILON;
6
7/// Dense symmetric positive-definite linear system `A x = b`.
8#[derive(Debug, Clone, PartialEq)]
9pub struct SpdLinearSystem {
10    dimension: usize,
11    matrix: Vec<f64>,
12    rhs: Vec<f64>,
13}
14
15impl SpdLinearSystem {
16    /// Construct and validate an SPD linear system from a row-major matrix.
17    ///
18    /// # Errors
19    ///
20    /// Returns [`LinearSolverError`] for invalid dimensions, non-finite values,
21    /// non-symmetry, or failure of Cholesky positive-definiteness validation.
22    pub fn new(matrix: &[f64], rhs: &[f64], dimension: usize) -> Result<Self, LinearSolverError> {
23        if dimension == 0 {
24            return Err(LinearSolverError::ZeroDimension);
25        }
26        let expected_len = dimension
27            .checked_mul(dimension)
28            .ok_or(LinearSolverError::MatrixDimensionMismatch)?;
29        if matrix.len() != expected_len {
30            return Err(LinearSolverError::MatrixDimensionMismatch);
31        }
32        if rhs.len() != dimension {
33            return Err(LinearSolverError::VectorDimensionMismatch);
34        }
35        if matrix.iter().any(|value| !value.is_finite()) {
36            return Err(LinearSolverError::NonFiniteMatrixEntry);
37        }
38        if rhs.iter().any(|value| !value.is_finite()) {
39            return Err(LinearSolverError::NonFiniteVectorEntry);
40        }
41
42        let matrix_view = DMatrix::from_row_slice(dimension, dimension, matrix);
43        let scale = matrix
44            .iter()
45            .fold(1.0_f64, |acc, value| acc.max(value.abs()));
46        let tolerance = SYMMETRY_TOLERANCE * scale;
47        for row in 0..dimension {
48            for column in (row + 1)..dimension {
49                if (matrix_view[(row, column)] - matrix_view[(column, row)]).abs() > tolerance {
50                    return Err(LinearSolverError::NonSymmetricSystemMatrix);
51                }
52            }
53        }
54        if matrix_view.cholesky().is_none() {
55            return Err(LinearSolverError::SystemMatrixNotPositiveDefinite);
56        }
57
58        Ok(Self {
59            dimension,
60            matrix: matrix.to_vec(),
61            rhs: rhs.to_vec(),
62        })
63    }
64
65    /// Return the system dimension.
66    #[must_use]
67    pub const fn dimension(&self) -> usize {
68        self.dimension
69    }
70
71    /// Return the row-major system matrix.
72    #[must_use]
73    pub fn matrix(&self) -> &[f64] {
74        &self.matrix
75    }
76
77    /// Return the right-hand side.
78    #[must_use]
79    pub fn rhs(&self) -> &[f64] {
80        &self.rhs
81    }
82}
83
84#[cfg(test)]
85mod tests {
86    use super::SpdLinearSystem;
87    use crate::LinearSolverError;
88
89    #[test]
90    fn accepts_spd_system() {
91        let system =
92            SpdLinearSystem::new(&[4.0, 1.0, 1.0, 3.0], &[1.0, 2.0], 2).expect("matrix is SPD");
93        assert_eq!(system.dimension(), 2);
94        assert_eq!(system.rhs(), &[1.0, 2.0]);
95    }
96
97    #[test]
98    fn rejects_non_symmetric_matrix() {
99        assert_eq!(
100            SpdLinearSystem::new(&[2.0, 1.0, 0.0, 2.0], &[1.0, 1.0], 2),
101            Err(LinearSolverError::NonSymmetricSystemMatrix)
102        );
103    }
104
105    #[test]
106    fn rejects_non_positive_definite_matrix() {
107        assert_eq!(
108            SpdLinearSystem::new(&[1.0, 2.0, 2.0, 1.0], &[1.0, 1.0], 2),
109            Err(LinearSolverError::SystemMatrixNotPositiveDefinite)
110        );
111    }
112}
uncertain_numerics/linear_system.rs

uncertain_numerics/
linear_system.rs
