tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it./workspace/linsolve/src/cholesky.rs
1 //! Symmetric **signed diagonally-pivoted root-free factorization** and the
2 //! associated indefinite linear solver.
3 //!
4 //! For a real symmetric matrix `A` this crate computes the factorization
5 //!
6 //! ```text
7 //! P A Pᵀ = L D Lᵀ
8 //! ```
9 //!
10 //! where `P` is a permutation, `L` is **unit lower triangular** (ones on the
11 //! diagonal) and `D` is a **diagonal matrix with signed real entries**. Unlike
12 //! the classical Cholesky factor this takes **no square roots** and handles
13 //! indefinite and negative-definite matrices directly: the signs of `D` are the
14 //! signs of `A`'s eigenvalues.
15 //!
16 //! The exact conventions (pivot rule, permutation meaning, rank / inertia /
17 //! singularity contract, solve and reconstruction semantics, tolerances and
18 //! anchor values) are specified in `instruction.md`. The numerical core of this
19 //! module is **unimplemented** in this build: [`factor_in_place`] and the
20 //! numerical methods of [`Factorization`] call [`todo!`]. Implementing them so
21 //! the behavioral contract holds is the task.
22
23 use crate::config::Config;
24 use crate::error::LinSolveError;
25 use crate::matrix::Matrix;
26
27 /// A computed signed diagonally-pivoted factorization `P A Pᵀ = L D Lᵀ`.
28 ///
29 /// Holds the unit-lower-triangular factor `L`, the signed diagonal `D`, and the
30 /// permutation `P` (as a vector). Reuse one `Factorization` to solve `A x = b`
31 /// for many right-hand sides without re-factorizing.
32 #[derive(Debug, Clone, PartialEq)]
33 pub struct Factorization {
34 /// Unit-lower-triangular factor `L` (diagonal entries are exactly `1.0`,
35 /// strictly-upper entries are `0.0`), stored in the permuted index space.
36 lower: Matrix,
37 /// Signed diagonal entries `D[0..n]` in the permuted index space.
38 diag: Vec<f64>,
39 /// The pivot permutation `P`, in the permuted index space.
40 perm: Vec<usize>,
41 /// Number of accepted (non-negligible) pivots; the rank under the pivot
42 /// floor used at factor time.
43 rank: usize,
44 /// Pivot floor used at factor time (entries with `|d| <= floor` count as
45 /// zero for rank / inertia).
46 floor: f64,
47 n: usize,
48 }
49
50 impl Factorization {
51 /// Factor the symmetric matrix `a` as `P A Pᵀ = L D Lᵀ`.
52 ///
53 /// With the default [`Config`] the input is first checked for squareness
54 /// and symmetry; then it is factored with symmetric pivoting.
55 ///
56 /// # Errors
57 ///
58 /// - [`LinSolveError::NonSquare`] if `a` is not square.
59 /// - [`LinSolveError::NotSymmetric`] if symmetry checking is enabled and
60 /// `a` is not symmetric within tolerance.
61 pub fn factor(a: &Matrix, config: &Config) -> Result<Self, LinSolveError> {
62 let n = a.require_square()?;
63 if config.check_symmetry() {
64 a.check_symmetric(config.symmetry_tolerance())?;
65 }
66
67 // Working symmetric matrix `m` (full n×n) in a permuted index space.
68 let mut m = Matrix::zeros(n, n);
69 for i in 0..n {
70 for j in 0..n {
71 // Mirror the lower triangle into the upper so swaps stay
72 // symmetric regardless of the caller's upper triangle.
73 let v = if i >= j { a.get(i, j) } else { a.get(j, i) };
74 m.set(i, j, v);
75 }
76 }
77
78 let mut lower = Matrix::zeros(n, n);
79 for i in 0..n {
80 lower.set(i, i, 1.0);
81 }
82 let mut diag = vec![0.0_f64; n];
83 let mut perm: Vec<usize> = (0..n).collect();
84 let floor = config.pivot_tolerance();
85
86 let rank = factor_in_place(&mut m, &mut lower, &mut diag, &mut perm, floor)?;
87
88 Ok(Self {
89 lower,
90 diag,
91 perm,
92 rank,
93 floor,
94 n,
95 })
96 }
97
98 /// The system dimension `n`.
99 #[inline]
100 pub fn dim(&self) -> usize {
101 self.n
102 }
103
104 /// Borrow the unit-lower-triangular factor `L` (permuted index space).
105 pub fn lower(&self) -> &Matrix {
106 &self.lower
107 }
108
109 /// Borrow the signed diagonal `D` (permuted index space).
110 pub fn diagonal(&self) -> &[f64] {
111 &self.diag
112 }
113
114 /// Borrow the pivot permutation `P`.
115 pub fn permutation(&self) -> &[usize] {
116 &self.perm
117 }
118
119 /// Numerical rank: number of pivots with magnitude strictly above the pivot
120 /// floor used at factor time.
121 #[inline]
122 pub fn rank(&self) -> usize {
123 self.rank
124 }
125
126 /// Whether the factorization is numerically nonsingular (`rank == dim`).
127 #[inline]
128 pub fn is_nonsingular(&self) -> bool {
129 self.rank == self.n
130 }
131
132 /// The inertia `(positive, negative, zero)` of the signed diagonal `D`.
133 pub fn inertia(&self) -> (usize, usize, usize) {
134 // Numerical core unimplemented in this build.
135 let _ = (&self.diag, self.floor);
136 todo!("compute the inertia (positive, negative, zero) of the signed diagonal")
137 }
138
139 /// Reconstruct `P A Pᵀ = L D Lᵀ` (the *permuted* matrix), as a full
140 /// symmetric matrix.
141 pub fn reconstruct_permuted(&self) -> Matrix {
142 // Numerical core unimplemented in this build.
143 let _ = (&self.lower, &self.diag, self.n);
144 todo!("reconstruct L D Lᵀ in the permuted index space")
145 }
146
147 /// Reconstruct the original matrix `A = Pᵀ (L D Lᵀ) P`, undoing the
148 /// permutation, as a full symmetric matrix.
149 pub fn reconstruct_original(&self) -> Matrix {
150 // Numerical core unimplemented in this build.
151 let _ = (&self.lower, &self.diag, &self.perm, self.n);
152 todo!("reconstruct A by undoing the permutation of L D Lᵀ")
153 }
154
155 /// The signed determinant `det(A)` (sign preserved).
156 pub fn signed_determinant(&self) -> f64 {
157 // Numerical core unimplemented in this build.
158 let _ = &self.diag;
159 todo!("compute the signed determinant from D")
160 }
161
162 /// `ln |det(A)|`, computed in log space. `f64::NEG_INFINITY` if singular.
163 pub fn log_abs_determinant(&self) -> f64 {
164 // Numerical core unimplemented in this build.
165 let _ = &self.diag;
166 todo!("compute ln|det(A)| in log space")
167 }
168
169 /// Solve `A x = b`, returning a fresh vector.
170 ///
171 /// # Errors
172 ///
173 /// - [`LinSolveError::RhsDimensionMismatch`] if `b.len() != dim()`.
174 /// - [`LinSolveError::SingularMatrix`] if the factorization is singular.
175 pub fn solve(&self, b: &[f64]) -> Result<Vec<f64>, LinSolveError> {
176 let mut x = b.to_vec();
177 self.solve_into(&mut x)?;
178 Ok(x)
179 }
180
181 /// Solve `A x = b` in place: `rhs` is overwritten with the solution `x`.
182 ///
183 /// # Errors
184 ///
185 /// - [`LinSolveError::RhsDimensionMismatch`] if `rhs.len() != dim()`.
186 /// - [`LinSolveError::SingularMatrix`] if the factorization is singular.
187 pub fn solve_into(&self, rhs: &mut [f64]) -> Result<(), LinSolveError> {
188 // Numerical core unimplemented in this build.
189 let _ = (&self.lower, &self.diag, &self.perm, self.n, self.rank, &rhs);
190 todo!("solve A x = b via the permuted L D Lᵀ factorization")
191 }
192 }
193
194 /// Numerical core: signed diagonally-pivoted root-free (`LDLᵀ`) factorization.
195 ///
196 /// On entry `m` holds a full symmetric copy of the input matrix `A` (lower and
197 /// upper triangles mirror each other), `lower` is the identity, `diag` is all
198 /// zeros, and `perm` is the identity permutation `[0, 1, …, n-1]`.
199 ///
200 /// On return the routine must have overwritten `lower` with the unit-lower
201 /// factor `L`, `diag` with the signed diagonal `D`, and `perm` with the
202 /// symmetric pivot order, all in the permuted index space, such that
203 /// `P A Pᵀ = L D Lᵀ`. It returns the rank (the number of pivots whose magnitude
204 /// is strictly greater than `pivot_floor`).
205 ///
206 /// The full behavioral contract , pivot rule, permutation meaning, and the
207 /// rank / singularity convention , is specified in `instruction.md`.
208 pub fn factor_in_place(
209 m: &mut Matrix,
210 lower: &mut Matrix,
211 diag: &mut [f64],
212 perm: &mut [usize],
213 pivot_floor: f64,
214 ) -> Result<usize, LinSolveError> {
215 // Numerical core unimplemented in this build.
216 let _ = (m, lower, diag, perm, pivot_floor);
217 todo!("implement the signed diagonally-pivoted LDLᵀ factorization in place")
218 }
219
/workspace/linsolve/src/matrix.rs
1 //! A minimal dense, row-major real matrix used by the solver.
2 //!
3 //! This is deliberately small: just enough structure to express a square
4 //! coefficient matrix, validate its shape, and perform the handful of
5 //! operations the Cholesky routines and their tests need (element access,
6 //! a symmetry check, and a matrix–vector product for residual computation).
7
8 use crate::error::LinSolveError;
9
10 /// A dense matrix of `f64` stored in row-major order.
11 ///
12 /// Construct one with [`Matrix::from_row_major`] (validates shape and
13 /// finiteness) or [`Matrix::zeros`]. Element access is via [`Matrix::get`] /
14 /// [`Matrix::set`]; both panic on out-of-bounds indices, matching the
15 /// convention of the standard library's indexing.
16 #[derive(Debug, Clone, PartialEq)]
17 pub struct Matrix {
18 rows: usize,
19 cols: usize,
20 data: Vec<f64>,
21 }
22
23 impl Matrix {
24 /// Build a matrix from a row-major data buffer.
25 ///
26 /// Returns [`LinSolveError::EmptyMatrix`] if either dimension is zero,
27 /// [`LinSolveError::DataShapeMismatch`] if `data.len() != rows * cols`, and
28 /// [`LinSolveError::NonFiniteEntry`] if any entry is `NaN`/infinite.
29 ///
30 /// ```
31 /// use linsolve::Matrix;
32 /// let m = Matrix::from_row_major(2, 2, vec![4.0, 1.0, 1.0, 3.0]).unwrap();
33 /// assert_eq!(m.get(0, 1), 1.0);
34 /// ```
35 pub fn from_row_major(rows: usize, cols: usize, data: Vec<f64>) -> Result<Self, LinSolveError> {
36 if rows == 0 || cols == 0 {
37 return Err(LinSolveError::EmptyMatrix);
38 }
39 if data.len() != rows * cols {
40 return Err(LinSolveError::DataShapeMismatch {
41 rows,
42 cols,
43 len: data.len(),
44 });
45 }
46 for i in 0..rows {
47 for j in 0..cols {
48 let v = data[i * cols + j];
49 if !v.is_finite() {
50 return Err(LinSolveError::NonFiniteEntry { i, j, value: v });
51 }
52 }
53 }
54 Ok(Self { rows, cols, data })
55 }
56
57 /// Build a `rows × cols` matrix of zeros.
58 pub fn zeros(rows: usize, cols: usize) -> Self {
59 Self {
60 rows,
61 cols,
62 data: vec![0.0; rows * cols],
63 }
64 }
65
66 /// Number of rows.
67 #[inline]
68 pub fn rows(&self) -> usize {
69 self.rows
70 }
71
72 /// Number of columns.
73 #[inline]
74 pub fn cols(&self) -> usize {
75 self.cols
76 }
77
78 /// Whether the matrix is square.
79 #[inline]
80 pub fn is_square(&self) -> bool {
81 self.rows == self.cols
82 }
83
84 /// Get the entry at `(row, col)`. Panics if out of bounds.
85 #[inline]
86 pub fn get(&self, row: usize, col: usize) -> f64 {
87 assert!(row < self.rows && col < self.cols, "index out of bounds");
88 self.data[row * self.cols + col]
89 }
90
91 /// Set the entry at `(row, col)`. Panics if out of bounds.
92 #[inline]
93 pub fn set(&mut self, row: usize, col: usize, value: f64) {
94 assert!(row < self.rows && col < self.cols, "index out of bounds");
95 self.data[row * self.cols + col] = value;
96 }
97
98 /// Check that the matrix is square, returning [`LinSolveError::NonSquare`]
99 /// otherwise.
100 pub fn require_square(&self) -> Result<usize, LinSolveError> {
101 if self.is_square() {
102 Ok(self.rows)
103 } else {
104 Err(LinSolveError::NonSquare {
105 rows: self.rows,
106 cols: self.cols,
107 })
108 }
109 }
110
111 /// Verify the matrix is symmetric to within absolute tolerance `tol`.
112 ///
113 /// Returns the worst asymmetry found (`0.0` for a perfectly symmetric
114 /// matrix), or [`LinSolveError::NotSymmetric`] if any off-diagonal pair
115 /// differs by more than `tol`.
116 pub fn check_symmetric(&self, tol: f64) -> Result<f64, LinSolveError> {
117 let n = self.require_square()?;
118 let mut worst = 0.0_f64;
119 for i in 0..n {
120 for j in (i + 1)..n {
121 let d = (self.get(i, j) - self.get(j, i)).abs();
122 if d > worst {
123 worst = d;
124 }
125 if d > tol {
126 return Err(LinSolveError::NotSymmetric {
127 i,
128 j,
129 asymmetry: d,
130 });
131 }
132 }
133 }
134 Ok(worst)
135 }
136
137 /// Compute the matrix–vector product `A * x`.
138 ///
139 /// Returns [`LinSolveError::RhsDimensionMismatch`] if `x.len() != cols`.
140 /// This is used by the tests to form the residual `A x - b`.
141 ///
142 /// ```
143 /// use linsolve::Matrix;
144 /// let a = Matrix::from_row_major(2, 2, vec![2.0, 0.0, 0.0, 3.0]).unwrap();
145 /// let y = a.matvec(&[1.0, 1.0]).unwrap();
146 /// assert_eq!(y, vec![2.0, 3.0]);
147 /// ```
148 pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, LinSolveError> {
149 if x.len() != self.cols {
150 return Err(LinSolveError::RhsDimensionMismatch {
151 expected: self.cols,
152 got: x.len(),
153 });
154 }
155 let mut out = vec![0.0; self.rows];
156 for i in 0..self.rows {
157 let mut acc = 0.0;
158 for j in 0..self.cols {
159 acc += self.get(i, j) * x[j];
160 }
161 out[i] = acc;
162 }
163 Ok(out)
164 }
165 }
166
1 //! Configuration for a symmetric factorization.
2
3 /// Default absolute tolerance for the symmetry check used by [`Config::new`].
4 pub const DEFAULT_SYMMETRY_TOLERANCE: f64 = 1e-9;
5
6 /// Default pivot floor used by [`Config::new`]: a pivot whose **absolute value**
7 /// does not exceed this is treated as negligible (zero) for rank/inertia, and
8 /// halts elimination. The default `0.0` accepts any nonzero pivot.
9 pub const DEFAULT_PIVOT_TOLERANCE: f64 = 0.0;
10
11 /// Tuning parameters for a symmetric factorization.
12 ///
13 /// A `Config` bundles the checks applied before and during factorization.
14 /// Construct one with [`Config::new`] (sensible defaults) and refine it with
15 /// the chained setters, e.g.
16 ///
17 /// ```
18 /// use linsolve::Config;
19 /// let cfg = Config::new()
20 /// .with_symmetry_tolerance(1e-12)
21 /// .with_pivot_tolerance(1e-14);
22 /// assert_eq!(cfg.symmetry_tolerance(), 1e-12);
23 /// ```
24 #[derive(Debug, Clone, Copy, PartialEq)]
25 pub struct Config {
26 check_symmetry: bool,
27 symmetry_tolerance: f64,
28 pivot_tolerance: f64,
29 }
30
31 impl Config {
32 /// Create a configuration with the crate defaults: symmetry checking is
33 /// **on** with tolerance [`DEFAULT_SYMMETRY_TOLERANCE`], and the pivot
34 /// floor is [`DEFAULT_PIVOT_TOLERANCE`].
35 pub fn new() -> Self {
36 Self {
37 check_symmetry: true,
38 symmetry_tolerance: DEFAULT_SYMMETRY_TOLERANCE,
39 pivot_tolerance: DEFAULT_PIVOT_TOLERANCE,
40 }
41 }
42
43 /// Enable or disable the pre-factorization symmetry check.
44 ///
45 /// When enabled (the default), [`crate::Factorization::factor`] verifies
46 /// that the input is symmetric to within
47 /// [`symmetry_tolerance`](Self::symmetry_tolerance). When disabled, the
48 /// matrix is symmetrized from its lower triangle with no validation.
49 #[must_use]
50 pub fn with_symmetry_check(mut self, enabled: bool) -> Self {
51 self.check_symmetry = enabled;
52 self
53 }
54
55 /// Set the absolute tolerance for the symmetry check.
56 #[must_use]
57 pub fn with_symmetry_tolerance(mut self, tol: f64) -> Self {
58 self.symmetry_tolerance = tol;
59 self
60 }
61
62 /// Set the pivot floor. During factorization, a candidate pivot whose
63 /// **absolute value** does not exceed this (or is non-finite) is treated as
64 /// negligible: it counts as a zero pivot for rank/inertia and halts
65 /// elimination. Raising it above zero rejects marginally-definite or
66 /// near-singular pivots.
67 #[must_use]
68 pub fn with_pivot_tolerance(mut self, tol: f64) -> Self {
69 self.pivot_tolerance = tol;
70 self
71 }
72
73 /// Whether the symmetry check is enabled.
74 pub fn check_symmetry(&self) -> bool {
75 self.check_symmetry
76 }
77
78 /// The configured symmetry tolerance.
79 pub fn symmetry_tolerance(&self) -> f64 {
80 self.symmetry_tolerance
81 }
82
83 /// The configured pivot floor.
84 pub fn pivot_tolerance(&self) -> f64 {
85 self.pivot_tolerance
86 }
87 }
88
89 impl Default for Config {
90 fn default() -> Self {
91 Self::new()
92 }
93 }
94
1 //! # linsolve
2 //!
3 //! Dense linear solver for **real symmetric** systems via a signed
4 //! diagonally-pivoted root-free factorization
5 //!
6 //! ```text
7 //! P A Pᵀ = L D Lᵀ
8 //! ```
9 //!
10 //! where `P` is a permutation, `L` is unit lower triangular, and `D` is a
11 //! diagonal matrix with **signed** real entries. This generalizes the classical
12 //! Cholesky factorization: it uses no square roots and works for **indefinite**
13 //! and **negative-definite** matrices, not just positive-definite ones. The
14 //! signs of `D` give the matrix's inertia.
15 //!
16 //! The crate exposes a small dense [`Matrix`] type, a [`Config`] for the
17 //! factorization checks, and the [`Factorization`] handle which factors a
18 //! matrix once and then solves `A x = b` for any number of right-hand sides.
19 //! Failure modes (non-square, non-symmetric, singular, mismatched right-hand
20 //! side, non-finite entry) are reported as a [`LinSolveError`].
21 //!
22 //! ```
23 //! use linsolve::{Config, Factorization, Matrix};
24 //!
25 //! // A symmetric *indefinite* system A x = b (eigenvalues of mixed sign).
26 //! let a = Matrix::from_row_major(3, 3, vec![
27 //! 1.0, 2.0, 3.0,
28 //! 2.0, 1.0, 4.0,
29 //! 3.0, 4.0, 1.0,
30 //! ]).unwrap();
31 //! let fac = Factorization::factor(&a, &Config::new()).unwrap();
32 //! // Mixed inertia: one positive, two negative.
33 //! assert_eq!(fac.inertia(), (1, 2, 0));
34 //! let x = fac.solve(&[1.0, 2.0, 3.0]).unwrap();
35 //! // Residual A x - b is tiny.
36 //! let r = a.matvec(&x).unwrap();
37 //! for (ri, bi) in r.iter().zip([1.0, 2.0, 3.0]) {
38 //! assert!((ri - bi).abs() < 1e-9);
39 //! }
40 //! ```
41 //!
42 //! Singular (rank-deficient under the pivot floor) matrices are detected rather
43 //! than silently producing nonsense; solving one is an error:
44 //!
45 //! ```
46 //! use linsolve::{Config, Factorization, LinSolveError, Matrix};
47 //! // Rank-1 PSD matrix v vᵀ with v = (1, 1): rank 1, one zero pivot.
48 //! let a = Matrix::from_row_major(2, 2, vec![1.0, 1.0, 1.0, 1.0]).unwrap();
49 //! let fac = Factorization::factor(&a, &Config::new()).unwrap();
50 //! assert_eq!(fac.rank(), 1);
51 //! assert!(!fac.is_nonsingular());
52 //! let err = fac.solve(&[1.0, 1.0]).unwrap_err();
53 //! assert!(matches!(err, LinSolveError::SingularMatrix { rank: 1 }));
54 //! ```
55 //!
56 //! See the [`cholesky`] module for the factorization handle and the numerical
57 //! core.
58
59 #![forbid(unsafe_code)]
60 #![warn(missing_docs)]
61
62 mod cholesky;
63 mod config;
64 mod error;
65 mod matrix;
66
67 pub use cholesky::{factor_in_place, Factorization};
68 pub use config::{Config, DEFAULT_PIVOT_TOLERANCE, DEFAULT_SYMMETRY_TOLERANCE};
69 pub use error::LinSolveError;
70 pub use matrix::Matrix;
71
72 /// Crate-level convenience: factor `a` with the default [`Config`] and solve
73 /// `A x = b` in one call.
74 ///
75 /// Equivalent to `Factorization::factor(a, &Config::new())?.solve(b)`.
76 ///
77 /// ```
78 /// use linsolve::{solve_symmetric, Matrix};
79 /// // 2x2 SPD system.
80 /// let a = Matrix::from_row_major(2, 2, vec![2.0, 0.0, 0.0, 4.0]).unwrap();
81 /// let x = solve_symmetric(&a, &[6.0, 8.0]).unwrap();
82 /// assert!((x[0] - 3.0).abs() < 1e-12);
83 /// assert!((x[1] - 2.0).abs() < 1e-12);
84 /// ```
85 pub fn solve_symmetric(a: &Matrix, b: &[f64]) -> Result<Vec<f64>, LinSolveError> {
86 Factorization::factor(a, &Config::new())?.solve(b)
87 }
88
/workspace/linsolve/src/config.rs
1 //! A minimal dense, row-major real matrix used by the solver.
2 //!
3 //! This is deliberately small: just enough structure to express a square
4 //! coefficient matrix, validate its shape, and perform the handful of
5 //! operations the Cholesky routines and their tests need (element access,
6 //! a symmetry check, and a matrix–vector product for residual computation).
7
8 use crate::error::LinSolveError;
9
10 /// A dense matrix of `f64` stored in row-major order.
11 ///
12 /// Construct one with [`Matrix::from_row_major`] (validates shape and
13 /// finiteness) or [`Matrix::zeros`]. Element access is via [`Matrix::get`] /
14 /// [`Matrix::set`]; both panic on out-of-bounds indices, matching the
15 /// convention of the standard library's indexing.
16 #[derive(Debug, Clone, PartialEq)]
17 pub struct Matrix {
18 rows: usize,
19 cols: usize,
20 data: Vec<f64>,
21 }
22
23 impl Matrix {
24 /// Build a matrix from a row-major data buffer.
25 ///
26 /// Returns [`LinSolveError::EmptyMatrix`] if either dimension is zero,
27 /// [`LinSolveError::DataShapeMismatch`] if `data.len() != rows * cols`, and
28 /// [`LinSolveError::NonFiniteEntry`] if any entry is `NaN`/infinite.
29 ///
30 /// ```
31 /// use linsolve::Matrix;
32 /// let m = Matrix::from_row_major(2, 2, vec![4.0, 1.0, 1.0, 3.0]).unwrap();
33 /// assert_eq!(m.get(0, 1), 1.0);
34 /// ```
35 pub fn from_row_major(rows: usize, cols: usize, data: Vec<f64>) -> Result<Self, LinSolveError> {
36 if rows == 0 || cols == 0 {
37 return Err(LinSolveError::EmptyMatrix);
38 }
39 if data.len() != rows * cols {
40 return Err(LinSolveError::DataShapeMismatch {
41 rows,
42 cols,
43 len: data.len(),
44 });
45 }
46 for i in 0..rows {
47 for j in 0..cols {
48 let v = data[i * cols + j];
49 if !v.is_finite() {
50 return Err(LinSolveError::NonFiniteEntry { i, j, value: v });
51 }
52 }
53 }
54 Ok(Self { rows, cols, data })
55 }
56
57 /// Build a `rows × cols` matrix of zeros.
58 pub fn zeros(rows: usize, cols: usize) -> Self {
59 Self {
60 rows,
61 cols,
62 data: vec![0.0; rows * cols],
63 }
64 }
65
66 /// Number of rows.
67 #[inline]
68 pub fn rows(&self) -> usize {
69 self.rows
70 }
71
72 /// Number of columns.
73 #[inline]
74 pub fn cols(&self) -> usize {
75 self.cols
76 }
77
78 /// Whether the matrix is square.
79 #[inline]
80 pub fn is_square(&self) -> bool {
81 self.rows == self.cols
82 }
83
84 /// Get the entry at `(row, col)`. Panics if out of bounds.
85 #[inline]
86 pub fn get(&self, row: usize, col: usize) -> f64 {
87 assert!(row < self.rows && col < self.cols, "index out of bounds");
88 self.data[row * self.cols + col]
89 }
90
91 /// Set the entry at `(row, col)`. Panics if out of bounds.
92 #[inline]
93 pub fn set(&mut self, row: usize, col: usize, value: f64) {
94 assert!(row < self.rows && col < self.cols, "index out of bounds");
95 self.data[row * self.cols + col] = value;
96 }
97
98 /// Check that the matrix is square, returning [`LinSolveError::NonSquare`]
99 /// otherwise.
100 pub fn require_square(&self) -> Result<usize, LinSolveError> {
101 if self.is_square() {
102 Ok(self.rows)
103 } else {
104 Err(LinSolveError::NonSquare {
105 rows: self.rows,
106 cols: self.cols,
107 })
108 }
109 }
110
111 /// Verify the matrix is symmetric to within absolute tolerance `tol`.
112 ///
113 /// Returns the worst asymmetry found (`0.0` for a perfectly symmetric
114 /// matrix), or [`LinSolveError::NotSymmetric`] if any off-diagonal pair
115 /// differs by more than `tol`.
116 pub fn check_symmetric(&self, tol: f64) -> Result<f64, LinSolveError> {
117 let n = self.require_square()?;
118 let mut worst = 0.0_f64;
119 for i in 0..n {
120 for j in (i + 1)..n {
121 let d = (self.get(i, j) - self.get(j, i)).abs();
122 if d > worst {
123 worst = d;
124 }
125 if d > tol {
126 return Err(LinSolveError::NotSymmetric {
127 i,
128 j,
129 asymmetry: d,
130 });
131 }
132 }
133 }
134 Ok(worst)
135 }
136
137 /// Compute the matrix–vector product `A * x`.
138 ///
139 /// Returns [`LinSolveError::RhsDimensionMismatch`] if `x.len() != cols`.
140 /// This is used by the tests to form the residual `A x - b`.
141 ///
142 /// ```
143 /// use linsolve::Matrix;
144 /// let a = Matrix::from_row_major(2, 2, vec![2.0, 0.0, 0.0, 3.0]).unwrap();
145 /// let y = a.matvec(&[1.0, 1.0]).unwrap();
146 /// assert_eq!(y, vec![2.0, 3.0]);
147 /// ```
148 pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, LinSolveError> {
149 if x.len() != self.cols {
150 return Err(LinSolveError::RhsDimensionMismatch {
151 expected: self.cols,
152 got: x.len(),
153 });
154 }
155 let mut out = vec![0.0; self.rows];
156 for i in 0..self.rows {
157 let mut acc = 0.0;
158 for j in 0..self.cols {
159 acc += self.get(i, j) * x[j];
160 }
161 out[i] = acc;
162 }
163 Ok(out)
164 }
165 }
166
1 //! Configuration for a symmetric factorization.
2
3 /// Default absolute tolerance for the symmetry check used by [`Config::new`].
4 pub const DEFAULT_SYMMETRY_TOLERANCE: f64 = 1e-9;
5
6 /// Default pivot floor used by [`Config::new`]: a pivot whose **absolute value**
7 /// does not exceed this is treated as negligible (zero) for rank/inertia, and
8 /// halts elimination. The default `0.0` accepts any nonzero pivot.
9 pub const DEFAULT_PIVOT_TOLERANCE: f64 = 0.0;
10
11 /// Tuning parameters for a symmetric factorization.
12 ///
13 /// A `Config` bundles the checks applied before and during factorization.
14 /// Construct one with [`Config::new`] (sensible defaults) and refine it with
15 /// the chained setters, e.g.
16 ///
17 /// ```
18 /// use linsolve::Config;
19 /// let cfg = Config::new()
20 /// .with_symmetry_tolerance(1e-12)
21 /// .with_pivot_tolerance(1e-14);
22 /// assert_eq!(cfg.symmetry_tolerance(), 1e-12);
23 /// ```
24 #[derive(Debug, Clone, Copy, PartialEq)]
25 pub struct Config {
26 check_symmetry: bool,
27 symmetry_tolerance: f64,
28 pivot_tolerance: f64,
29 }
30
31 impl Config {
32 /// Create a configuration with the crate defaults: symmetry checking is
33 /// **on** with tolerance [`DEFAULT_SYMMETRY_TOLERANCE`], and the pivot
34 /// floor is [`DEFAULT_PIVOT_TOLERANCE`].
35 pub fn new() -> Self {
36 Self {
37 check_symmetry: true,
38 symmetry_tolerance: DEFAULT_SYMMETRY_TOLERANCE,
39 pivot_tolerance: DEFAULT_PIVOT_TOLERANCE,
40 }
41 }
42
43 /// Enable or disable the pre-factorization symmetry check.
44 ///
45 /// When enabled (the default), [`crate::Factorization::factor`] verifies
46 /// that the input is symmetric to within
47 /// [`symmetry_tolerance`](Self::symmetry_tolerance). When disabled, the
48 /// matrix is symmetrized from its lower triangle with no validation.
49 #[must_use]
50 pub fn with_symmetry_check(mut self, enabled: bool) -> Self {
51 self.check_symmetry = enabled;
52 self
53 }
54
55 /// Set the absolute tolerance for the symmetry check.
56 #[must_use]
57 pub fn with_symmetry_tolerance(mut self, tol: f64) -> Self {
58 self.symmetry_tolerance = tol;
59 self
60 }
61
62 /// Set the pivot floor. During factorization, a candidate pivot whose
63 /// **absolute value** does not exceed this (or is non-finite) is treated as
64 /// negligible: it counts as a zero pivot for rank/inertia and halts
65 /// elimination. Raising it above zero rejects marginally-definite or
66 /// near-singular pivots.
67 #[must_use]
68 pub fn with_pivot_tolerance(mut self, tol: f64) -> Self {
69 self.pivot_tolerance = tol;
70 self
71 }
72
73 /// Whether the symmetry check is enabled.
74 pub fn check_symmetry(&self) -> bool {
75 self.check_symmetry
76 }
77
78 /// The configured symmetry tolerance.
79 pub fn symmetry_tolerance(&self) -> f64 {
80 self.symmetry_tolerance
81 }
82
83 /// The configured pivot floor.
84 pub fn pivot_tolerance(&self) -> f64 {
85 self.pivot_tolerance
86 }
87 }
88
89 impl Default for Config {
90 fn default() -> Self {
91 Self::new()
92 }
93 }
94
1 //! # linsolve
2 //!
3 //! Dense linear solver for **real symmetric** systems via a signed
4 //! diagonally-pivoted root-free factorization
5 //!
6 //! ```text
7 //! P A Pᵀ = L D Lᵀ
8 //! ```
9 //!
10 //! where `P` is a permutation, `L` is unit lower triangular, and `D` is a
11 //! diagonal matrix with **signed** real entries. This generalizes the classical
12 //! Cholesky factorization: it uses no square roots and works for **indefinite**
13 //! and **negative-definite** matrices, not just positive-definite ones. The
14 //! signs of `D` give the matrix's inertia.
15 //!
16 //! The crate exposes a small dense [`Matrix`] type, a [`Config`] for the
17 //! factorization checks, and the [`Factorization`] handle which factors a
18 //! matrix once and then solves `A x = b` for any number of right-hand sides.
19 //! Failure modes (non-square, non-symmetric, singular, mismatched right-hand
20 //! side, non-finite entry) are reported as a [`LinSolveError`].
21 //!
22 //! ```
23 //! use linsolve::{Config, Factorization, Matrix};
24 //!
25 //! // A symmetric *indefinite* system A x = b (eigenvalues of mixed sign).
26 //! let a = Matrix::from_row_major(3, 3, vec![
27 //! 1.0, 2.0, 3.0,
28 //! 2.0, 1.0, 4.0,
29 //! 3.0, 4.0, 1.0,
30 //! ]).unwrap();
31 //! let fac = Factorization::factor(&a, &Config::new()).unwrap();
32 //! // Mixed inertia: one positive, two negative.
33 //! assert_eq!(fac.inertia(), (1, 2, 0));
34 //! let x = fac.solve(&[1.0, 2.0, 3.0]).unwrap();
35 //! // Residual A x - b is tiny.
36 //! let r = a.matvec(&x).unwrap();
37 //! for (ri, bi) in r.iter().zip([1.0, 2.0, 3.0]) {
38 //! assert!((ri - bi).abs() < 1e-9);
39 //! }
40 //! ```
41 //!
42 //! Singular (rank-deficient under the pivot floor) matrices are detected rather
43 //! than silently producing nonsense; solving one is an error:
44 //!
45 //! ```
46 //! use linsolve::{Config, Factorization, LinSolveError, Matrix};
47 //! // Rank-1 PSD matrix v vᵀ with v = (1, 1): rank 1, one zero pivot.
48 //! let a = Matrix::from_row_major(2, 2, vec![1.0, 1.0, 1.0, 1.0]).unwrap();
49 //! let fac = Factorization::factor(&a, &Config::new()).unwrap();
50 //! assert_eq!(fac.rank(), 1);
51 //! assert!(!fac.is_nonsingular());
52 //! let err = fac.solve(&[1.0, 1.0]).unwrap_err();
53 //! assert!(matches!(err, LinSolveError::SingularMatrix { rank: 1 }));
54 //! ```
55 //!
56 //! See the [`cholesky`] module for the factorization handle and the numerical
57 //! core.
58
59 #![forbid(unsafe_code)]
60 #![warn(missing_docs)]
61
62 mod cholesky;
63 mod config;
64 mod error;
65 mod matrix;
66
67 pub use cholesky::{factor_in_place, Factorization};
68 pub use config::{Config, DEFAULT_PIVOT_TOLERANCE, DEFAULT_SYMMETRY_TOLERANCE};
69 pub use error::LinSolveError;
70 pub use matrix::Matrix;
71
72 /// Crate-level convenience: factor `a` with the default [`Config`] and solve
73 /// `A x = b` in one call.
74 ///
75 /// Equivalent to `Factorization::factor(a, &Config::new())?.solve(b)`.
76 ///
77 /// ```
78 /// use linsolve::{solve_symmetric, Matrix};
79 /// // 2x2 SPD system.
80 /// let a = Matrix::from_row_major(2, 2, vec![2.0, 0.0, 0.0, 4.0]).unwrap();
81 /// let x = solve_symmetric(&a, &[6.0, 8.0]).unwrap();
82 /// assert!((x[0] - 3.0).abs() < 1e-12);
83 /// assert!((x[1] - 2.0).abs() < 1e-12);
84 /// ```
85 pub fn solve_symmetric(a: &Matrix, b: &[f64]) -> Result<Vec<f64>, LinSolveError> {
86 Factorization::factor(a, &Config::new())?.solve(b)
87 }
88
/workspace/linsolve/src/lib.rs
1 //! A minimal dense, row-major real matrix used by the solver.
2 //!
3 //! This is deliberately small: just enough structure to express a square
4 //! coefficient matrix, validate its shape, and perform the handful of
5 //! operations the Cholesky routines and their tests need (element access,
6 //! a symmetry check, and a matrix–vector product for residual computation).
7
8 use crate::error::LinSolveError;
9
10 /// A dense matrix of `f64` stored in row-major order.
11 ///
12 /// Construct one with [`Matrix::from_row_major`] (validates shape and
13 /// finiteness) or [`Matrix::zeros`]. Element access is via [`Matrix::get`] /
14 /// [`Matrix::set`]; both panic on out-of-bounds indices, matching the
15 /// convention of the standard library's indexing.
16 #[derive(Debug, Clone, PartialEq)]
17 pub struct Matrix {
18 rows: usize,
19 cols: usize,
20 data: Vec<f64>,
21 }
22
23 impl Matrix {
24 /// Build a matrix from a row-major data buffer.
25 ///
26 /// Returns [`LinSolveError::EmptyMatrix`] if either dimension is zero,
27 /// [`LinSolveError::DataShapeMismatch`] if `data.len() != rows * cols`, and
28 /// [`LinSolveError::NonFiniteEntry`] if any entry is `NaN`/infinite.
29 ///
30 /// ```
31 /// use linsolve::Matrix;
32 /// let m = Matrix::from_row_major(2, 2, vec![4.0, 1.0, 1.0, 3.0]).unwrap();
33 /// assert_eq!(m.get(0, 1), 1.0);
34 /// ```
35 pub fn from_row_major(rows: usize, cols: usize, data: Vec<f64>) -> Result<Self, LinSolveError> {
36 if rows == 0 || cols == 0 {
37 return Err(LinSolveError::EmptyMatrix);
38 }
39 if data.len() != rows * cols {
40 return Err(LinSolveError::DataShapeMismatch {
41 rows,
42 cols,
43 len: data.len(),
44 });
45 }
46 for i in 0..rows {
47 for j in 0..cols {
48 let v = data[i * cols + j];
49 if !v.is_finite() {
50 return Err(LinSolveError::NonFiniteEntry { i, j, value: v });
51 }
52 }
53 }
54 Ok(Self { rows, cols, data })
55 }
56
57 /// Build a `rows × cols` matrix of zeros.
58 pub fn zeros(rows: usize, cols: usize) -> Self {
59 Self {
60 rows,
61 cols,
62 data: vec![0.0; rows * cols],
63 }
64 }
65
66 /// Number of rows.
67 #[inline]
68 pub fn rows(&self) -> usize {
69 self.rows
70 }
71
72 /// Number of columns.
73 #[inline]
74 pub fn cols(&self) -> usize {
75 self.cols
76 }
77
78 /// Whether the matrix is square.
79 #[inline]
80 pub fn is_square(&self) -> bool {
81 self.rows == self.cols
82 }
83
84 /// Get the entry at `(row, col)`. Panics if out of bounds.
85 #[inline]
86 pub fn get(&self, row: usize, col: usize) -> f64 {
87 assert!(row < self.rows && col < self.cols, "index out of bounds");
88 self.data[row * self.cols + col]
89 }
90
91 /// Set the entry at `(row, col)`. Panics if out of bounds.
92 #[inline]
93 pub fn set(&mut self, row: usize, col: usize, value: f64) {
94 assert!(row < self.rows && col < self.cols, "index out of bounds");
95 self.data[row * self.cols + col] = value;
96 }
97
98 /// Check that the matrix is square, returning [`LinSolveError::NonSquare`]
99 /// otherwise.
100 pub fn require_square(&self) -> Result<usize, LinSolveError> {
101 if self.is_square() {
102 Ok(self.rows)
103 } else {
104 Err(LinSolveError::NonSquare {
105 rows: self.rows,
106 cols: self.cols,
107 })
108 }
109 }
110
111 /// Verify the matrix is symmetric to within absolute tolerance `tol`.
112 ///
113 /// Returns the worst asymmetry found (`0.0` for a perfectly symmetric
114 /// matrix), or [`LinSolveError::NotSymmetric`] if any off-diagonal pair
115 /// differs by more than `tol`.
116 pub fn check_symmetric(&self, tol: f64) -> Result<f64, LinSolveError> {
117 let n = self.require_square()?;
118 let mut worst = 0.0_f64;
119 for i in 0..n {
120 for j in (i + 1)..n {
121 let d = (self.get(i, j) - self.get(j, i)).abs();
122 if d > worst {
123 worst = d;
124 }
125 if d > tol {
126 return Err(LinSolveError::NotSymmetric {
127 i,
128 j,
129 asymmetry: d,
130 });
131 }
132 }
133 }
134 Ok(worst)
135 }
136
137 /// Compute the matrix–vector product `A * x`.
138 ///
139 /// Returns [`LinSolveError::RhsDimensionMismatch`] if `x.len() != cols`.
140 /// This is used by the tests to form the residual `A x - b`.
141 ///
142 /// ```
143 /// use linsolve::Matrix;
144 /// let a = Matrix::from_row_major(2, 2, vec![2.0, 0.0, 0.0, 3.0]).unwrap();
145 /// let y = a.matvec(&[1.0, 1.0]).unwrap();
146 /// assert_eq!(y, vec![2.0, 3.0]);
147 /// ```
148 pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, LinSolveError> {
149 if x.len() != self.cols {
150 return Err(LinSolveError::RhsDimensionMismatch {
151 expected: self.cols,
152 got: x.len(),
153 });
154 }
155 let mut out = vec![0.0; self.rows];
156 for i in 0..self.rows {
157 let mut acc = 0.0;
158 for j in 0..self.cols {
159 acc += self.get(i, j) * x[j];
160 }
161 out[i] = acc;
162 }
163 Ok(out)
164 }
165 }
166
1 //! Configuration for a symmetric factorization.
2
3 /// Default absolute tolerance for the symmetry check used by [`Config::new`].
4 pub const DEFAULT_SYMMETRY_TOLERANCE: f64 = 1e-9;
5
6 /// Default pivot floor used by [`Config::new`]: a pivot whose **absolute value**
7 /// does not exceed this is treated as negligible (zero) for rank/inertia, and
8 /// halts elimination. The default `0.0` accepts any nonzero pivot.
9 pub const DEFAULT_PIVOT_TOLERANCE: f64 = 0.0;
10
11 /// Tuning parameters for a symmetric factorization.
12 ///
13 /// A `Config` bundles the checks applied before and during factorization.
14 /// Construct one with [`Config::new`] (sensible defaults) and refine it with
15 /// the chained setters, e.g.
16 ///
17 /// ```
18 /// use linsolve::Config;
19 /// let cfg = Config::new()
20 /// .with_symmetry_tolerance(1e-12)
21 /// .with_pivot_tolerance(1e-14);
22 /// assert_eq!(cfg.symmetry_tolerance(), 1e-12);
23 /// ```
24 #[derive(Debug, Clone, Copy, PartialEq)]
25 pub struct Config {
26 check_symmetry: bool,
27 symmetry_tolerance: f64,
28 pivot_tolerance: f64,
29 }
30
31 impl Config {
32 /// Create a configuration with the crate defaults: symmetry checking is
33 /// **on** with tolerance [`DEFAULT_SYMMETRY_TOLERANCE`], and the pivot
34 /// floor is [`DEFAULT_PIVOT_TOLERANCE`].
35 pub fn new() -> Self {
36 Self {
37 check_symmetry: true,
38 symmetry_tolerance: DEFAULT_SYMMETRY_TOLERANCE,
39 pivot_tolerance: DEFAULT_PIVOT_TOLERANCE,
40 }
41 }
42
43 /// Enable or disable the pre-factorization symmetry check.
44 ///
45 /// When enabled (the default), [`crate::Factorization::factor`] verifies
46 /// that the input is symmetric to within
47 /// [`symmetry_tolerance`](Self::symmetry_tolerance). When disabled, the
48 /// matrix is symmetrized from its lower triangle with no validation.
49 #[must_use]
50 pub fn with_symmetry_check(mut self, enabled: bool) -> Self {
51 self.check_symmetry = enabled;
52 self
53 }
54
55 /// Set the absolute tolerance for the symmetry check.
56 #[must_use]
57 pub fn with_symmetry_tolerance(mut self, tol: f64) -> Self {
58 self.symmetry_tolerance = tol;
59 self
60 }
61
62 /// Set the pivot floor. During factorization, a candidate pivot whose
63 /// **absolute value** does not exceed this (or is non-finite) is treated as
64 /// negligible: it counts as a zero pivot for rank/inertia and halts
65 /// elimination. Raising it above zero rejects marginally-definite or
66 /// near-singular pivots.
67 #[must_use]
68 pub fn with_pivot_tolerance(mut self, tol: f64) -> Self {
69 self.pivot_tolerance = tol;
70 self
71 }
72
73 /// Whether the symmetry check is enabled.
74 pub fn check_symmetry(&self) -> bool {
75 self.check_symmetry
76 }
77
78 /// The configured symmetry tolerance.
79 pub fn symmetry_tolerance(&self) -> f64 {
80 self.symmetry_tolerance
81 }
82
83 /// The configured pivot floor.
84 pub fn pivot_tolerance(&self) -> f64 {
85 self.pivot_tolerance
86 }
87 }
88
89 impl Default for Config {
90 fn default() -> Self {
91 Self::new()
92 }
93 }
94
1 //! # linsolve
2 //!
3 //! Dense linear solver for **real symmetric** systems via a signed
4 //! diagonally-pivoted root-free factorization
5 //!
6 //! ```text
7 //! P A Pᵀ = L D Lᵀ
8 //! ```
9 //!
10 //! where `P` is a permutation, `L` is unit lower triangular, and `D` is a
11 //! diagonal matrix with **signed** real entries. This generalizes the classical
12 //! Cholesky factorization: it uses no square roots and works for **indefinite**
13 //! and **negative-definite** matrices, not just positive-definite ones. The
14 //! signs of `D` give the matrix's inertia.
15 //!
16 //! The crate exposes a small dense [`Matrix`] type, a [`Config`] for the
17 //! factorization checks, and the [`Factorization`] handle which factors a
18 //! matrix once and then solves `A x = b` for any number of right-hand sides.
19 //! Failure modes (non-square, non-symmetric, singular, mismatched right-hand
20 //! side, non-finite entry) are reported as a [`LinSolveError`].
21 //!
22 //! ```
23 //! use linsolve::{Config, Factorization, Matrix};
24 //!
25 //! // A symmetric *indefinite* system A x = b (eigenvalues of mixed sign).
26 //! let a = Matrix::from_row_major(3, 3, vec![
27 //! 1.0, 2.0, 3.0,
28 //! 2.0, 1.0, 4.0,
29 //! 3.0, 4.0, 1.0,
30 //! ]).unwrap();
31 //! let fac = Factorization::factor(&a, &Config::new()).unwrap();
32 //! // Mixed inertia: one positive, two negative.
33 //! assert_eq!(fac.inertia(), (1, 2, 0));
34 //! let x = fac.solve(&[1.0, 2.0, 3.0]).unwrap();
35 //! // Residual A x - b is tiny.
36 //! let r = a.matvec(&x).unwrap();
37 //! for (ri, bi) in r.iter().zip([1.0, 2.0, 3.0]) {
38 //! assert!((ri - bi).abs() < 1e-9);
39 //! }
40 //! ```
41 //!
42 //! Singular (rank-deficient under the pivot floor) matrices are detected rather
43 //! than silently producing nonsense; solving one is an error:
44 //!
45 //! ```
46 //! use linsolve::{Config, Factorization, LinSolveError, Matrix};
47 //! // Rank-1 PSD matrix v vᵀ with v = (1, 1): rank 1, one zero pivot.
48 //! let a = Matrix::from_row_major(2, 2, vec![1.0, 1.0, 1.0, 1.0]).unwrap();
49 //! let fac = Factorization::factor(&a, &Config::new()).unwrap();
50 //! assert_eq!(fac.rank(), 1);
51 //! assert!(!fac.is_nonsingular());
52 //! let err = fac.solve(&[1.0, 1.0]).unwrap_err();
53 //! assert!(matches!(err, LinSolveError::SingularMatrix { rank: 1 }));
54 //! ```
55 //!
56 //! See the [`cholesky`] module for the factorization handle and the numerical
57 //! core.
58
59 #![forbid(unsafe_code)]
60 #![warn(missing_docs)]
61
62 mod cholesky;
63 mod config;
64 mod error;
65 mod matrix;
66
67 pub use cholesky::{factor_in_place, Factorization};
68 pub use config::{Config, DEFAULT_PIVOT_TOLERANCE, DEFAULT_SYMMETRY_TOLERANCE};
69 pub use error::LinSolveError;
70 pub use matrix::Matrix;
71
72 /// Crate-level convenience: factor `a` with the default [`Config`] and solve
73 /// `A x = b` in one call.
74 ///
75 /// Equivalent to `Factorization::factor(a, &Config::new())?.solve(b)`.
76 ///
77 /// ```
78 /// use linsolve::{solve_symmetric, Matrix};
79 /// // 2x2 SPD system.
80 /// let a = Matrix::from_row_major(2, 2, vec![2.0, 0.0, 0.0, 4.0]).unwrap();
81 /// let x = solve_symmetric(&a, &[6.0, 8.0]).unwrap();
82 /// assert!((x[0] - 3.0).abs() < 1e-12);
83 /// assert!((x[1] - 2.0).abs() < 1e-12);
84 /// ```
85 pub fn solve_symmetric(a: &Matrix, b: &[f64]) -> Result<Vec<f64>, LinSolveError> {
86 Factorization::factor(a, &Config::new())?.solve(b)
87 }
88
/workspace/linsolve/src/error.rs
1 //! Error types for the `linsolve` crate.
2
3 use thiserror::Error;
4
5 /// Errors that can arise while building a matrix, factorizing it, or solving a
6 /// system.
7 ///
8 /// These cover the cases where the requested operation cannot be carried out
9 /// reliably (a malformed matrix, a mismatched right-hand side, or a singular
10 /// system). They are intentionally distinct from a *successful* factorization,
11 /// whose result is returned through the [`Factorization`](crate::Factorization)
12 /// handle.
13 #[derive(Debug, Error, Clone, PartialEq)]
14 #[non_exhaustive]
15 pub enum LinSolveError {
16 /// A matrix was constructed with zero rows or zero columns. Factorization
17 /// requires at least a 1×1 matrix.
18 #[error("matrix must have at least one row and one column")]
19 EmptyMatrix,
20
21 /// The data buffer length did not equal `rows * cols`.
22 #[error("data length {len} does not match dimensions {rows}x{cols}")]
23 DataShapeMismatch {
24 /// Number of rows requested.
25 rows: usize,
26 /// Number of columns requested.
27 cols: usize,
28 /// Length of the supplied data buffer.
29 len: usize,
30 },
31
32 /// The factorization requires a square matrix; this one was not.
33 #[error("matrix must be square, got {rows}x{cols}")]
34 NonSquare {
35 /// Number of rows.
36 rows: usize,
37 /// Number of columns.
38 cols: usize,
39 },
40
41 /// The matrix is not symmetric within the configured tolerance, so the
42 /// symmetric factorization does not apply. Reports the worst offending
43 /// entry pair `(i, j)`.
44 #[error("matrix is not symmetric: |A[{i},{j}] - A[{j},{i}]| = {asymmetry} exceeds tolerance")]
45 NotSymmetric {
46 /// Row index of the offending entry.
47 i: usize,
48 /// Column index of the offending entry.
49 j: usize,
50 /// Magnitude of the asymmetry `|A[i,j] - A[j,i]|`.
51 asymmetry: f64,
52 },
53
54 /// A solve was requested against a singular factorization (one or more
55 /// negligible pivots, i.e. `rank < dim`), which has no unique solution.
56 /// Reports the numerical rank that was found.
57 #[error("matrix is singular: rank {rank} is less than the dimension")]
58 SingularMatrix {
59 /// The numerical rank found by the factorization.
60 rank: usize,
61 },
62
63 /// A right-hand side had a length that did not match the system dimension.
64 #[error("right-hand side length {got} does not match system dimension {expected}")]
65 RhsDimensionMismatch {
66 /// The system dimension `n`.
67 expected: usize,
68 /// The supplied right-hand side length.
69 got: usize,
70 },
71
72 /// A supplied matrix entry was not a finite number (it was `NaN` or an
73 /// infinity), so the factorization cannot proceed.
74 #[error("matrix entry at [{i},{j}] is not finite: {value}")]
75 NonFiniteEntry {
76 /// Row index of the offending entry.
77 i: usize,
78 /// Column index of the offending entry.
79 j: usize,
80 /// The non-finite value.
81 value: f64,
82 },
83 }
84
/workspace/linsolve:
Cargo.lock
Cargo.toml
src
tests
/workspace/linsolve/src:
cholesky.rs
config.rs
error.rs
lib.rs
matrix.rs
/workspace/linsolve/tests:
integration.rs
[stdout]
/workspace/linsolve:
Cargo.lock
Cargo.toml
src
tests
/workspace/linsolve/src:
cholesky.rs
config.rs
error.rs
lib.rs
matrix.rs
/workspace/linsolve/tests:
integration.rsls -R /workspace/linsolve --ignore=target
1 //! Error types for the `linsolve` crate.
2
3 use thiserror::Error;
4
5 /// Errors that can arise while building a matrix, factorizing it, or solving a
6 /// system.
7 ///
8 /// These cover the cases where the requested operation cannot be carried out
9 /// reliably (a malformed matrix, a mismatched right-hand side, or a singular
10 /// system). They are intentionally distinct from a *successful* factorization,
11 /// whose result is returned through the [`Factorization`](crate::Factorization)
12 /// handle.
13 #[derive(Debug, Error, Clone, PartialEq)]
14 #[non_exhaustive]
15 pub enum LinSolveError {
16 /// A matrix was constructed with zero rows or zero columns. Factorization
17 /// requires at least a 1×1 matrix.
18 #[error("matrix must have at least one row and one column")]
19 EmptyMatrix,
20
21 /// The data buffer length did not equal `rows * cols`.
22 #[error("data length {len} does not match dimensions {rows}x{cols}")]
23 DataShapeMismatch {
24 /// Number of rows requested.
25 rows: usize,
26 /// Number of columns requested.
27 cols: usize,
28 /// Length of the supplied data buffer.
29 len: usize,
30 },
31
32 /// The factorization requires a square matrix; this one was not.
33 #[error("matrix must be square, got {rows}x{cols}")]
34 NonSquare {
35 /// Number of rows.
36 rows: usize,
37 /// Number of columns.
38 cols: usize,
39 },
40
41 /// The matrix is not symmetric within the configured tolerance, so the
42 /// symmetric factorization does not apply. Reports the worst offending
43 /// entry pair `(i, j)`.
44 #[error("matrix is not symmetric: |A[{i},{j}] - A[{j},{i}]| = {asymmetry} exceeds tolerance")]
45 NotSymmetric {
46 /// Row index of the offending entry.
47 i: usize,
48 /// Column index of the offending entry.
49 j: usize,
50 /// Magnitude of the asymmetry `|A[i,j] - A[j,i]|`.
51 asymmetry: f64,
52 },
53
54 /// A solve was requested against a singular factorization (one or more
55 /// negligible pivots, i.e. `rank < dim`), which has no unique solution.
56 /// Reports the numerical rank that was found.
57 #[error("matrix is singular: rank {rank} is less than the dimension")]
58 SingularMatrix {
59 /// The numerical rank found by the factorization.
60 rank: usize,
61 },
62
63 /// A right-hand side had a length that did not match the system dimension.
64 #[error("right-hand side length {got} does not match system dimension {expected}")]
65 RhsDimensionMismatch {
66 /// The system dimension `n`.
67 expected: usize,
68 /// The supplied right-hand side length.
69 got: usize,
70 },
71
72 /// A supplied matrix entry was not a finite number (it was `NaN` or an
73 /// infinity), so the factorization cannot proceed.
74 #[error("matrix entry at [{i},{j}] is not finite: {value}")]
75 NonFiniteEntry {
76 /// Row index of the offending entry.
77 i: usize,
78 /// Column index of the offending entry.
79 j: usize,
80 /// The non-finite value.
81 value: f64,
82 },
83 }
84
/workspace/linsolve:
Cargo.lock
Cargo.toml
src
tests
/workspace/linsolve/src:
cholesky.rs
config.rs
error.rs
lib.rs
matrix.rs
/workspace/linsolve/tests:
integration.rs
[stdout]
/workspace/linsolve:
Cargo.lock
Cargo.toml
src
tests
/workspace/linsolve/src:
cholesky.rs
config.rs
error.rs
lib.rs
matrix.rs
/workspace/linsolve/tests:
integration.rs/workspace/linsolve/tests/integration.rs
1 //! Smoke test for the `linsolve` crate.
2 //!
3 //! This is a minimal placeholder that only checks the data plumbing that does
4 //! not depend on the (currently unimplemented) numerical core: matrix
5 //! construction, shape validation, and the matrix–vector product. The full
6 //! behavioral test suite is applied separately by the grader.
7
8 use linsolve::{Config, LinSolveError, Matrix};
9
10 #[test]
11 fn matrix_construction_and_access() {
12 let m = Matrix::from_row_major(2, 2, vec![4.0, 1.0, 1.0, 3.0]).expect("valid");
13 assert_eq!(m.rows(), 2);
14 assert_eq!(m.cols(), 2);
15 assert!(m.is_square());
16 assert_eq!(m.get(0, 1), 1.0);
17 }
18
19 #[test]
20 fn matvec_smoke() {
21 let a = Matrix::from_row_major(2, 2, vec![2.0, 0.0, 0.0, 3.0]).unwrap();
22 let y = a.matvec(&[1.0, 1.0]).unwrap();
23 assert_eq!(y, vec![2.0, 3.0]);
24 }
25
26 #[test]
27 fn construction_errors() {
28 assert_eq!(
29 Matrix::from_row_major(0, 0, vec![]).unwrap_err(),
30 LinSolveError::EmptyMatrix
31 );
32 assert!(matches!(
33 Matrix::from_row_major(2, 2, vec![1.0, 2.0, 3.0]).unwrap_err(),
34 LinSolveError::DataShapeMismatch { len: 3, .. }
35 ));
36 }
37
38 #[test]
39 fn config_defaults() {
40 let cfg = Config::new();
41 assert!(cfg.check_symmetry());
42 assert_eq!(cfg.pivot_tolerance(), 0.0);
43 }
44
cargo build --release 2>&1 | tail -5
Finished `release` profile [optimized] target(s) in 0.21s
[stdout]
Finished `release` profile [optimized] target(s) in 0.21s/workspace/linsolve/src/cholesky.rs
pub fn factor_in_place(
m: &mut Matrix,
lower: &mut Matrix,
diag: &mut [f64],
perm: &mut [usize],
pivot_floor: f64,
) -> Result<usize, LinSolveError> {
// Numerical core unimplemented in this build.
let _ = (m, lower, diag, perm, pivot_floor);
todo!("implement the signed diagonally-pivoted LDLᵀ factorization in place")
}pub fn factor_in_place(
m: &mut Matrix,
lower: &mut Matrix,
diag: &mut [f64],
perm: &mut [usize],
pivot_floor: f64,
) -> Result<usize, LinSolveError> {
let n = perm.len();
// Default: full rank. Overwritten if elimination halts early.
let mut rank = n;
for k in 0..n {
// Symmetric (diagonal) pivot selection: pick the remaining position
// whose diagonal entry of the running Schur complement has the largest
// absolute value. Exact ties resolve to the smallest index (strict `>`).
let mut p = k;
let mut best = m.get(k, k).abs();
for i in (k + 1)..n {
let cand = m.get(i, i).abs();
if cand > best {
best = cand;
p = i;
}
}
// Halt if every remaining candidate pivot is negligible. The trailing
// block stays the seeded null space (zero diagonal, identity columns).
let pivot = m.get(p, p);
if !pivot.is_finite() || pivot.abs() <= pivot_floor {
rank = k;
break;
}
// Bring the chosen pivot to position k with a symmetric row/column swap.
if p != k {
for c in 0..n {
let a = m.get(k, c);
let b = m.get(p, c);
m.set(k, c, b);
m.set(p, c, a);
}
for r in 0..n {
let a = m.get(r, k);
let b = m.get(r, p);
m.set(r, k, b);
m.set(r, p, a);
}
perm.swap(k, p);
// Keep the already-computed columns of L aligned with the swap.
for c in 0..k {
let a = lower.get(k, c);
let b = lower.get(p, c);
lower.set(k, c, b);
lower.set(p, c, a);
}
}
let d = m.get(k, k);
diag[k] = d;
// Column k of the unit-lower factor.
for i in (k + 1)..n {
lower.set(i, k, m.get(i, k) / d);
}
// Rank-1 symmetric update of the trailing Schur complement.
for i in (k + 1)..n {
let mik = m.get(i, k);
for j in (k + 1)..n {
let updated = m.get(i, j) - mik * m.get(k, j) / d;
m.set(i, j, updated);
}
}
}
Ok(rank)
}/workspace/linsolve/src/cholesky.rs
pub fn inertia(&self) -> (usize, usize, usize) {
// Numerical core unimplemented in this build.
let _ = (&self.diag, self.floor);
todo!("compute the inertia (positive, negative, zero) of the signed diagonal")
} pub fn inertia(&self) -> (usize, usize, usize) {
let mut positive = 0;
let mut negative = 0;
let mut zero = 0;
for &d in &self.diag {
if d > self.floor {
positive += 1;
} else if d < -self.floor {
negative += 1;
} else {
zero += 1;
}
}
(positive, negative, zero)
}/workspace/linsolve/src/cholesky.rs
pub fn reconstruct_permuted(&self) -> Matrix {
// Numerical core unimplemented in this build.
let _ = (&self.lower, &self.diag, self.n);
todo!("reconstruct L D Lᵀ in the permuted index space")
} pub fn reconstruct_permuted(&self) -> Matrix {
let n = self.n;
let mut out = Matrix::zeros(n, n);
for i in 0..n {
for j in 0..n {
// L is unit lower triangular, so L[i][k] = 0 for k > i.
let kmax = i.min(j);
let mut acc = 0.0;
for k in 0..=kmax {
acc += self.lower.get(i, k) * self.diag[k] * self.lower.get(j, k);
}
out.set(i, j, acc);
}
}
out
}/workspace/linsolve/src/cholesky.rs
pub fn reconstruct_original(&self) -> Matrix {
// Numerical core unimplemented in this build.
let _ = (&self.lower, &self.diag, &self.perm, self.n);
todo!("reconstruct A by undoing the permutation of L D Lᵀ")
} pub fn reconstruct_original(&self) -> Matrix {
let n = self.n;
let permuted = self.reconstruct_permuted();
// `perm[i]` is the original index occupying factored position `i`, so
// scatter each permuted entry back to its original location.
let mut out = Matrix::zeros(n, n);
for i in 0..n {
for j in 0..n {
out.set(self.perm[i], self.perm[j], permuted.get(i, j));
}
}
out
}/workspace/linsolve/src/cholesky.rs
pub fn signed_determinant(&self) -> f64 {
// Numerical core unimplemented in this build.
let _ = &self.diag;
todo!("compute the signed determinant from D")
}
/// `ln |det(A)|`, computed in log space. `f64::NEG_INFINITY` if singular.
pub fn log_abs_determinant(&self) -> f64 {
// Numerical core unimplemented in this build.
let _ = &self.diag;
todo!("compute ln|det(A)| in log space")
} pub fn signed_determinant(&self) -> f64 {
// det(A) = det(P A Pᵀ) = det(L D Lᵀ) = ∏ D[k] (L is unit triangular).
self.diag.iter().product()
}
/// `ln |det(A)|`, computed in log space. `f64::NEG_INFINITY` if singular.
pub fn log_abs_determinant(&self) -> f64 {
if self.rank < self.n {
return f64::NEG_INFINITY;
}
self.diag.iter().map(|d| d.abs().ln()).sum()
}/workspace/linsolve/src/cholesky.rs
pub fn solve_into(&self, rhs: &mut [f64]) -> Result<(), LinSolveError> {
// Numerical core unimplemented in this build.
let _ = (&self.lower, &self.diag, &self.perm, self.n, self.rank, &rhs);
todo!("solve A x = b via the permuted L D Lᵀ factorization")
} pub fn solve_into(&self, rhs: &mut [f64]) -> Result<(), LinSolveError> {
let n = self.n;
if rhs.len() != n {
return Err(LinSolveError::RhsDimensionMismatch {
expected: n,
got: rhs.len(),
});
}
if self.rank < n {
return Err(LinSolveError::SingularMatrix { rank: self.rank });
}
// A = Pᵀ (L D Lᵀ) P, so with y = P x and c = P b we solve L D Lᵀ y = c.
// Permute the right-hand side into the factored index space.
let mut y: Vec<f64> = (0..n).map(|k| rhs[self.perm[k]]).collect();
// Forward substitution: L z = c (L is unit lower triangular).
for i in 0..n {
let mut acc = y[i];
for j in 0..i {
acc -= self.lower.get(i, j) * y[j];
}
y[i] = acc;
}
// Diagonal solve: D w = z.
for i in 0..n {
y[i] /= self.diag[i];
}
// Back substitution: Lᵀ y = w.
for i in (0..n).rev() {
let mut acc = y[i];
for j in (i + 1)..n {
acc -= self.lower.get(j, i) * y[j];
}
y[i] = acc;
}
// Scatter back: x = Pᵀ y, i.e. x[perm[k]] = y[k].
for k in 0..n {
rhs[self.perm[k]] = y[k];
}
Ok(())
}cargo build --release 2>&1 | tail -15 && echo "---TEST---" && cargo test --release 2>&1 | tail -25
Compiling linsolve v0.4.2 (/workspace/linsolve)
Finished `release` profile [optimized] target(s) in 1.23s
---TEST---
test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Running tests/integration.rs (target/release/deps/integration-126f5e260c8b96c6)
running 4 tests
test config_defaults ... ok
test construction_errors ... ok
test matrix_construction_and_access ... ok
test matvec_smoke ... ok
test result: ok. 4 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Doc-tests linsolve
running 6 tests
test src/config.rs - config::Config (line 17) ... ok
test src/lib.rs - (line 22) ... ok
test src/lib.rs - (line 45) ... ok
test src/lib.rs - solve_symmetric (line 77) ... ok
test src/matrix.rs - matrix::Matrix::from_row_major (line 30) ... ok
test src/matrix.rs - matrix::Matrix::matvec (line 142) ... ok
test result: ok. 6 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.40s
[stdout]
Compiling linsolve v0.4.2 (/workspace/linsolve)
Finished `release` profile [optimized] target(s) in 1.23s
---TEST---
test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Running tests/integration.rs (target/release/deps/integration-126f5e260c8b96c6)
running 4 tests
test config_defaults ... ok
test construction_errors ... ok
test matrix_construction_and_access ... ok
test matvec_smoke ... ok
test result: ok. 4 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Doc-tests linsolve
running 6 tests
test src/config.rs - config::Config (line 17) ... ok
test src/lib.rs - (line 22) ... ok
test src/lib.rs - (line 45) ... ok
test src/lib.rs - solve_symmetric (line 77) ... ok
test src/matrix.rs - matrix::Matrix::from_row_major (line 30) ... ok
test src/matrix.rs - matrix::Matrix::matvec (line 142) ... ok
test result: ok. 6 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.40scat > /tmp/anchor_test.rs << 'EOF'
use linsolve::{Config, Factorization, Matrix, LinSolveError};
fn approx(a: f64, b: f64) -> bool { (a-b).abs() < 1e-7 }
fn main() {
// SPD
let a = Matrix::from_row_major(3,3, vec![4.,12.,-16.,12.,37.,-43.,-16.,-43.,98.]).unwrap();
let f = Factorization::factor(&a, &Config::new()).unwrap();
assert_eq!(f.rank(),3); assert_eq!(f.inertia(),(3,0,0));
assert!(approx(f.signed_determinant(),36.), "spd det {}", f.signed_determinant());
// reconstruct_original == A
let r = f.reconstruct_original();
for i in 0..3 { for j in 0..3 { assert!(approx(r.get(i,j), a.get(i,j))); }}
// diag
let a = Matrix::from_row_major(3,3, vec![2.,0.,0.,0.,4.,0.,0.,0.,5.]).unwrap();
let f = Factorization::factor(&a, &Config::new()).unwrap();
assert_eq!(f.rank(),3); assert_eq!(f.inertia(),(3,0,0));
assert!(approx(f.signed_determinant(),40.));
// L = I
for i in 0..3 { for j in 0..3 { let e = if i==j {1.0} else {0.0}; assert!(approx(f.lower().get(i,j), e)); }}
// indefinite
let a = Matrix::from_row_major(3,3, vec![1.,2.,3.,2.,1.,4.,3.,4.,1.]).unwrap();
let f = Factorization::factor(&a, &Config::new()).unwrap();
assert_eq!(f.rank(),3); assert_eq!(f.inertia(),(1,2,0));
assert!(approx(f.signed_determinant(),20.));
let x = f.solve(&[1.,2.,3.]).unwrap();
assert!(approx(x[0],1.)&&approx(x[1],0.)&&approx(x[2],0.), "x={:?}",x);
// neg def -I2
let a = Matrix::from_row_major(2,2, vec![-1.,0.,0.,-1.]).unwrap();
let f = Factorization::factor(&a, &Config::new()).unwrap();
assert_eq!(f.inertia(),(0,2,0));
let x = f.solve(&[2.,2.]).unwrap();
assert!(approx(x[0],-2.)&&approx(x[1],-2.));
// zero leading diag
let a = Matrix::from_row_major(3,3, vec![0.,3.,1.,3.,5.,2.,1.,2.,9.]).unwrap();
let f = Factorization::factor(&a, &Config::new()).unwrap();
assert_eq!(f.rank(),3); assert_eq!(f.inertia(),(2,1,0));
let r = f.reconstruct_original();
for i in 0..3 { for j in 0..3 { assert!(approx(r.get(i,j), a.get(i,j)),"recon {} {}",i,j); }}
// rank1 psd
let a = Matrix::from_row_major(2,2, vec![1.,1.,1.,1.]).unwrap();
let f = Factorization::factor(&a, &Config::new()).unwrap();
assert_eq!(f.rank(),1); assert_eq!(f.inertia(),(1,0,1));
assert!(matches!(f.solve(&[1.,1.]).unwrap_err(), LinSolveError::SingularMatrix{rank:1}));
// [[0,1],[1,0]]
let a = Matrix::from_row_major(2,2, vec![0.,1.,1.,0.]).unwrap();
let f = Factorization::factor(&a, &Config::new()).unwrap();
assert_eq!(f.rank(),0); assert_eq!(f.inertia(),(0,0,2));
// diag(1,1e-12)
let a = Matrix::from_row_major(2,2, vec![1.,0.,0.,1e-12]).unwrap();
let f = Factorization::factor(&a, &Config::new()).unwrap();
assert_eq!(f.rank(),2);
let f = Factorization::factor(&a, &Config::new().with_pivot_tolerance(1e-9)).unwrap();
assert_eq!(f.rank(),1); assert_eq!(f.inertia(),(1,0,1));
// Hilbert 5x5 solve accuracy + non-involution perm check via reconstruct
let n=5; let mut data=vec![0.0;n*n];
for i in 0..n { for j in 0..n { data[i*n+j]=1.0/((i+j+1) as f64); }}
let a = Matrix::from_row_major(n,n,data).unwrap();
let f = Factorization::factor(&a,&Config::new()).unwrap();
let xtrue=vec![1.,2.,3.,4.,5.];
let b=a.matvec(&xtrue).unwrap();
let x=f.solve(&b).unwrap();
for i in 0..n { assert!((x[i]-xtrue[i]).abs()<1e-6,"hilbert x[{}]={}",i,x[i]); }
let r=f.reconstruct_original();
for i in 0..n { for j in 0..n { assert!(approx(r.get(i,j),a.get(i,j)),"hilb recon"); }}
// reconstruct_permuted symmetric
let p=f.reconstruct_permuted();
for i in 0..n { for j in 0..n { assert_eq!(p.get(i,j),p.get(j,i)); }}
println!("perm = {:?}", f.permutation());
println!("ALL ANCHORS PASS");
}
EOF
mkdir -p examples && cp /tmp/anchor_test.rs examples/anchor_test.rs && cargo run --release --example anchor_test 2>&1 | tail -5; rm -f examples/anchor_test.rs; rmdir examples 2>/dev/null; truethread 'main' panicked at examples/anchor_test.rs:76:37: assertion `left == right` failed left: 0.125 right: 0.12499999999999999 note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace [stdout] thread 'main' panicked at examples/anchor_test.rs:76:37: assertion `left == right` failed left: 0.125 right: 0.12499999999999999 note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
/workspace/linsolve/src/cholesky.rs
let mut out = Matrix::zeros(n, n);
for i in 0..n {
for j in 0..n {
// L is unit lower triangular, so L[i][k] = 0 for k > i.
let kmax = i.min(j);
let mut acc = 0.0;
for k in 0..=kmax {
acc += self.lower.get(i, k) * self.diag[k] * self.lower.get(j, k);
}
out.set(i, j, acc);
}
}
out let mut out = Matrix::zeros(n, n);
// Compute the lower triangle and mirror it, so the result is exactly
// symmetric (immune to floating-point association order).
for i in 0..n {
for j in 0..=i {
// L is unit lower triangular, so L[i][k] = 0 for k > i.
let mut acc = 0.0;
for k in 0..=j {
acc += self.lower.get(i, k) * self.diag[k] * self.lower.get(j, k);
}
out.set(i, j, acc);
out.set(j, i, acc);
}
}
outmkdir -p examples && cp /tmp/anchor_test.rs examples/anchor_test.rs && cargo run --release --example anchor_test 2>&1 | tail -5; rm -f examples/anchor_test.rs; rmdir examples 2>/dev/null; true
Compiling linsolve v0.4.2 (/workspace/linsolve)
Finished `release` profile [optimized] target(s) in 0.65s
Running `target/release/examples/anchor_test`
perm = [0, 2, 4, 1, 3]
ALL ANCHORS PASS
[stdout]
Compiling linsolve v0.4.2 (/workspace/linsolve)
Finished `release` profile [optimized] target(s) in 0.65s
Running `target/release/examples/anchor_test`
perm = [0, 2, 4, 1, 3]
ALL ANCHORS PASScargo clean -p linsolve 2>/dev/null; cargo build --release 2>&1 | grep -iE "warning|error" | head; echo "build done (empty above = no warnings)"; cargo test --release 2>&1 | grep -E "test result|error"
build done (empty above = no warnings) test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s test construction_errors ... ok test result: ok. 4 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s test result: ok. 6 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.40s [stdout] build done (empty above = no warnings) test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s test construction_errors ... ok test result: ok. 4 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s test result: ok. 6 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.40s
No file changes, the agent left the workspace untouched.
[test] injected hidden suite into /workspace/linsolve/tests/integration.rs
[test] running cargo test --release (offline)
Compiling linsolve v0.4.2 (/workspace/linsolve)
Finished `release` profile [optimized] target(s) in 1.40s
Running unittests src/lib.rs (target/release/deps/linsolve-235cd0f53c08b5b3)
running 0 tests
test result: ok. 0 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Running tests/integration.rs (target/release/deps/integration-126f5e260c8b96c6)
running 50 tests
test discriminator_permutation_is_nontrivial ... ok
test discriminator_factors_indefinite_without_error ... ok
test discriminator_reconstruct_differs_from_naive_llt ... ok
test discriminator_unit_diagonal_not_sqrt ... ok
test factor_4x4_spd_anchored ... ok
test factor_diagonal_orders_by_descending_magnitude ... ok
test factor_indefinite_3x3_signed_diagonal ... FAILED
test factor_negative_definite_identity ... ok
test factor_pivots_around_zero_diagonal ... ok
test factor_spd_3x3_pivots_to_largest_diagonal ... ok
test factor_tie_break_orientation ... FAILED
test factor_unit_lower_2x2 ... ok
test factors_4x4_indefinite_with_full_diagnostics ... ok
test indefinite_with_zero_usable_diagonal_is_singular ... ok
test inertia_spd_all_positive ... ok
test log_abs_determinant_is_computed_in_log_space ... ok
test multiway_pivot_ties_follow_reference_orientation ... FAILED
test one_factorization_serves_multiple_rhs ... ok
test pivot_floor_boundary_classifies_inertia ... ok
test pivot_floor_rejects_marginally_definite ... ok
test reconstruct_original_recovers_a ... ok
test reconstruct_original_recovers_indefinite_a ... ok
test reconstruct_permuted_equals_p_a_pt ... ok
test rejects_asymmetric ... ok
test rejects_empty_matrix_at_construction ... ok
test rejects_non_finite_entry_at_construction ... ok
test rejects_non_square ... ok
test rejects_shape_mismatch_at_construction ... ok
test residual_tiny_across_many_rhs ... ok
test round_trip_constructed_rhs_spd ... ok
test signed_determinant_is_negative_for_indefinite ... ok
test signed_determinant_negative_definite ... ok
test signed_determinant_spd ... ok
test singular_rank1_psd_detected ... ok
test solution_is_linear_in_rhs_indefinite ... ok
test solve_into_overwrites_rhs ... ok
test solve_rhs_dimension_mismatch ... ok
test solve_singular_is_error ... ok
test solves_diagonal_system ... ok
test solves_hilbert_5_within_tolerance ... ok
test solves_indefinite_system ... ok
test solves_known_2x2 ... ok
test solves_negative_definite_system ... ok
test solves_one_by_one ... ok
test solves_pivot_around_zero_diagonal_system ... ok
test symmetry_check_can_be_disabled_uses_lower_triangle ... ok
test tie_orientation_all_equal_3x3 ... FAILED
test tie_orientation_all_equal_4x4 ... FAILED
test tie_orientation_interior_pair ... FAILED
test tie_orientation_paired_magnitudes ... FAILED
failures:
---- factor_indefinite_3x3_signed_diagonal stdout ----
thread 'factor_indefinite_3x3_signed_diagonal' panicked at tests/integration.rs:163:5:
assertion `left == right` failed
left: [0, 2, 1]
right: [2, 1, 0]
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
---- factor_tie_break_orientation stdout ----
thread 'factor_tie_break_orientation' panicked at tests/integration.rs:104:5:
assertion `left == right` failed
left: [0, 1, 2]
right: [2, 0, 1]
---- multiway_pivot_ties_follow_reference_orientation stdout ----
thread 'multiway_pivot_ties_follow_reference_orientation' panicked at tests/integration.rs:700:5:
assertion `left == right` failed
left: [0, 1, 2]
right: [1, 0, 2]
---- tie_orientation_all_equal_3x3 stdout ----
thread 'tie_orientation_all_equal_3x3' panicked at tests/integration.rs:752:5:
assertion `left == right` failed
left: [0, 1, 2]
right: [2, 0, 1]
---- tie_orientation_all_equal_4x4 stdout ----
thread 'tie_orientation_all_equal_4x4' panicked at tests/integration.rs:723:5:
assertion `left == right` failed
left: [0, 1, 2, 3]
right: [3, 0, 1, 2]
---- tie_orientation_interior_pair stdout ----
thread 'tie_orientation_interior_pair' panicked at tests/integration.rs:744:5:
assertion `left == right` failed
left: [1, 2, 3, 0]
right: [2, 1, 3, 0]
---- tie_orientation_paired_magnitudes stdout ----
thread 'tie_orientation_paired_magnitudes' panicked at tests/integration.rs:733:5:
assertion `left == right` failed
left: [2, 3, 0, 1]
right: [3, 2, 0, 1]
failures:
factor_indefinite_3x3_signed_diagonal
factor_tie_break_orientation
multiway_pivot_ties_follow_reference_orientation
tie_orientation_all_equal_3x3
tie_orientation_all_equal_4x4
tie_orientation_interior_pair
tie_orientation_paired_magnitudes
test result: FAILED. 43 passed; 7 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
error: test failed, to rerun pass `--test integration`
[test] tests FAILED
[test] reward = 0Reproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_433c38a8db16412f. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_433c38a8db16412f · verifier authoritative; classifier explanatory.