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
/workspace/linsolve/src/config.rs
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
ls tests/ src/ 2>/dev/null && echo "---ERROR---" && cat src/error.rs && echo "---LIB---" && cat src/lib.rs
src/:
cholesky.rs
config.rs
error.rs
lib.rs
matrix.rs
tests/:
integration.rs
---ERROR---
//! Error types for the `linsolve` crate.
use thiserror::Error;
/// Errors that can arise while building a matrix, factorizing it, or solving a
/// system.
///
/// These cover the cases where the requested operation cannot be carried out
/// reliably (a malformed matrix, a mismatched right-hand side, or a singular
/// system). They are intentionally distinct from a *successful* factorization,
/// whose result is returned through the [`Factorization`](crate::Factorization)
/// handle.
#[derive(Debug, Error, Clone, PartialEq)]
#[non_exhaustive]
pub enum LinSolveError {
/// A matrix was constructed with zero rows or zero columns. Factorization
/// requires at least a 1×1 matrix.
#[error("matrix must have at least one row and one column")]
EmptyMatrix,
/// The data buffer length did not equal `rows * cols`.
#[error("data length {len} does not match dimensions {rows}x{cols}")]
DataShapeMismatch {
/// Number of rows requested.
rows: usize,
/// Number of columns requested.
cols: usize,
/// Length of the supplied data buffer.
len: usize,
},
/// The factorization requires a square matrix; this one was not.
#[error("matrix must be square, got {rows}x{cols}")]
NonSquare {
/// Number of rows.
rows: usize,
/// Number of columns.
cols: usize,
},
/// The matrix is not symmetric within the configured tolerance, so the
/// symmetric factorization does not apply. Reports the worst offending
/// entry pair `(i, j)`.
#[error("matrix is not symmetric: |A[{i},{j}] - A[{j},{i}]| = {asymmetry} exceeds tolerance")]
NotSymmetric {
/// Row index of the offending entry.
i: usize,
/// Column index of the offending entry.
j: usize,
/// Magnitude of the asymmetry `|A[i,j] - A[j,i]|`.
asymmetry: f64,
},
/// A solve was requested against a singular factorization (one or more
/// negligible pivots, i.e. `rank < dim`), which has no unique solution.
/// Reports the numerical rank that was found.
#[error("matrix is singular: rank {rank} is less than the dimension")]
SingularMatrix {
/// The numerical rank found by the factorization.
rank: usize,
},
/// A right-hand side had a length that did not match the system dimension.
#[error("right-hand side length {got} does not match system dimension {expected}")]
RhsDimensionMismatch {
/// The system dimension `n`.
expected: usize,
/// The supplied right-hand side length.
got: usize,
},
/// A supplied matrix entry was not a finite number (it was `NaN` or an
/// infinity), so the factorization cannot proceed.
#[error("matrix entry at [{i},{j}] is not finite: {value}")]
NonFiniteEntry {
/// Row index of the offending entry.
i: usize,
/// Column index of the offending entry.
j: usize,
/// The non-finite value.
value: f64,
},
}
---LIB---
//! # linsolve
//!
//! Dense linear solver for **real symmetric** systems via a signed
//! diagonally-pivoted root-free factorization
//!
//! ```text
//! P A Pᵀ = L D Lᵀ
//! ```
//!
//! where `P` is a permutation, `L` is unit lower triangular, and `D` is a
//! diagonal matrix with **signed** real entries. This generalizes the classical
//! Cholesky factorization: it uses no square roots and works for **indefinite**
//! and **negative-definite** matrices, not just positive-definite ones. The
//! signs of `D` give the matrix's inertia.
//!
//! The crate exposes a small dense [`Matrix`] type, a [`Config`] for the
//! factorization checks, and the [`Factorization`] handle which factors a
//! matrix once and then solves `A x = b` for any number of right-hand sides.
//! Failure modes (non-square, non-symmetric, singular, mismatched right-hand
//! side, non-finite entry) are reported as a [`LinSolveError`].
//!
//! ```
//! use linsolve::{Config, Factorization, Matrix};
//!
//! // A symmetric *indefinite* system A x = b (eigenvalues of mixed sign).
//! let a = Matrix::from_row_major(3, 3, vec![
//! 1.0, 2.0, 3.0,
//! 2.0, 1.0, 4.0,
//! 3.0, 4.0, 1.0,
//! ]).unwrap();
//! let fac = Factorization::factor(&a, &Config::new()).unwrap();
//! // Mixed inertia: one positive, two negative.
//! assert_eq!(fac.inertia(), (1, 2, 0));
//! let x = fac.solve(&[1.0, 2.0, 3.0]).unwrap();
//! // Residual A x - b is tiny.
//! let r = a.matvec(&x).unwrap();
//! for (ri, bi) in r.iter().zip([1.0, 2.0, 3.0]) {
//! assert!((ri - bi).abs() < 1e-9);
//! }
//! ```
//!
//! Singular (rank-deficient under the pivot floor) matrices are detected rather
//! than silently producing nonsense; solving one is an error:
//!
//! ```
//! use linsolve::{Config, Factorization, LinSolveError, Matrix};
//! // Rank-1 PSD matrix v vᵀ with v = (1, 1): rank 1, one zero pivot.
//! let a = Matrix::from_row_major(2, 2, vec![1.0, 1.0, 1.0, 1.0]).unwrap();
//! let fac = Factorization::factor(&a, &Config::new()).unwrap();
//! assert_eq!(fac.rank(), 1);
//! assert!(!fac.is_nonsingular());
//! let err = fac.solve(&[1.0, 1.0]).unwrap_err();
//! assert!(matches!(err, LinSolveError::SingularMatrix { rank: 1 }));
//! ```
//!
//! See the [`cholesky`] module for the factorization handle and the numerical
//! core.
#![forbid(unsafe_code)]
#![warn(missing_docs)]
mod cholesky;
mod config;
mod error;
mod matrix;
pub use cholesky::{factor_in_place, Factorization};
pub use config::{Config, DEFAULT_PIVOT_TOLERANCE, DEFAULT_SYMMETRY_TOLERANCE};
pub use error::LinSolveError;
pub use matrix::Matrix;
/// Crate-level convenience: factor `a` with the default [`Config`] and solve
/// `A x = b` in one call.
///
/// Equivalent to `Factorization::factor(a, &Config::new())?.solve(b)`.
///
/// ```
/// use linsolve::{solve_symmetric, Matrix};
/// // 2x2 SPD system.
/// let a = Matrix::from_row_major(2, 2, vec![2.0, 0.0, 0.0, 4.0]).unwrap();
/// let x = solve_symmetric(&a, &[6.0, 8.0]).unwrap();
/// assert!((x[0] - 3.0).abs() < 1e-12);
/// assert!((x[1] - 2.0).abs() < 1e-12);
/// ```
pub fn solve_symmetric(a: &Matrix, b: &[f64]) -> Result<Vec<f64>, LinSolveError> {
Factorization::factor(a, &Config::new())?.solve(b)
}
[stdout]
src/:
cholesky.rs
config.rs
error.rs
lib.rs
matrix.rs
tests/:
integration.rs
---ERROR---
//! Error types for the `linsolve` crate.
use thiserror::Error;
/// Errors that can arise while building a matrix, factorizing it, or solving a
/// system.
///
/// These cover the cases where the requested operation cannot be carried out
/// reliably (a malformed matrix, a mismatched right-hand side, or a singular
/// system). They are intentionally distinct from a *successful* factorization,
/// whose result is returned through the [`Factorization`](crate::Factorization)
/// handle.
#[derive(Debug, Error, Clone, PartialEq)]
#[non_exhaustive]
pub enum LinSolveError {
/// A matrix was constructed with zero rows or zero columns. Factorization
/// requires at least a 1×1 matrix.
#[error("matrix must have at least one row and one column")]
EmptyMatrix,
/// The data buffer length did not equal `rows * cols`.
#[error("data length {len} does not match dimensions {rows}x{cols}")]
DataShapeMismatch {
/// Number of rows requested.
rows: usize,
/// Number of columns requested.
cols: usize,
/// Length of the supplied data buffer.
len: usize,
},
/// The factorization requires a square matrix; this one was not.
#[error("matrix must be square, got {rows}x{cols}")]
NonSquare {
/// Number of rows.
rows: usize,
/// Number of columns.
cols: usize,
},
/// The matrix is not symmetric within the configured tolerance, so the
/// symmetric factorization does not apply. Reports the worst offending
/// entry pair `(i, j)`.
#[error("matrix is not symmetric: |A[{i},{j}] - A[{j},{i}]| = {asymmetry} exceeds tolerance")]
NotSymmetric {
/// Row index of the offending entry.
i: usize,
/// Column index of the offending entry.
j: usize,
/// Magnitude of the asymmetry `|A[i,j] - A[j,i]|`.
asymmetry: f64,
},
/// A solve was requested against a singular factorization (one or more
/// negligible pivots, i.e. `rank < dim`), which has no unique solution.
/// Reports the numerical rank that was found.
#[error("matrix is singular: rank {rank} is less than the dimension")]
SingularMatrix {
/// The numerical rank found by the factorization.
rank: usize,
},
/// A right-hand side had a length that did not match the system dimension.
#[error("right-hand side length {got} does not match system dimension {expected}")]
RhsDimensionMismatch {
/// The system dimension `n`.
expected: usize,
/// The supplied right-hand side length.
got: usize,
},
/// A supplied matrix entry was not a finite number (it was `NaN` or an
/// infinity), so the factorization cannot proceed.
#[error("matrix entry at [{i},{j}] is not finite: {value}")]
NonFiniteEntry {
/// Row index of the offending entry.
i: usize,
/// Column index of the offending entry.
j: usize,
/// The non-finite value.
value: f64,
},
}
---LIB---
//! # linsolve
//!
//! Dense linear solver for **real symmetric** systems via a signed
//! diagonally-pivoted root-free factorization
//!
//! ```text
//! P A Pᵀ = L D Lᵀ
//! ```
//!
//! where `P` is a permutation, `L` is unit lower triangular, and `D` is a
//! diagonal matrix with **signed** real entries. This generalizes the classical
//! Cholesky factorization: it uses no square roots and works for **indefinite**
//! and **negative-definite** matrices, not just positive-definite ones. The
//! signs of `D` give the matrix's inertia.
//!
//! The crate exposes a small dense [`Matrix`] type, a [`Config`] for the
//! factorization checks, and the [`Factorization`] handle which factors a
//! matrix once and then solves `A x = b` for any number of right-hand sides.
//! Failure modes (non-square, non-symmetric, singular, mismatched right-hand
//! side, non-finite entry) are reported as a [`LinSolveError`].
//!
//! ```
//! use linsolve::{Config, Factorization, Matrix};
//!
//! // A symmetric *indefinite* system A x = b (eigenvalues of mixed sign).
//! let a = Matrix::from_row_major(3, 3, vec![
//! 1.0, 2.0, 3.0,
//! 2.0, 1.0, 4.0,
//! 3.0, 4.0, 1.0,
//! ]).unwrap();
//! let fac = Factorization::factor(&a, &Config::new()).unwrap();
//! // Mixed inertia: one positive, two negative.
//! assert_eq!(fac.inertia(), (1, 2, 0));
//! let x = fac.solve(&[1.0, 2.0, 3.0]).unwrap();
//! // Residual A x - b is tiny.
//! let r = a.matvec(&x).unwrap();
//! for (ri, bi) in r.iter().zip([1.0, 2.0, 3.0]) {
//! assert!((ri - bi).abs() < 1e-9);
//! }
//! ```
//!
//! Singular (rank-deficient under the pivot floor) matrices are detected rather
//! than silently producing nonsense; solving one is an error:
//!
//! ```
//! use linsolve::{Config, Factorization, LinSolveError, Matrix};
//! // Rank-1 PSD matrix v vᵀ with v = (1, 1): rank 1, one zero pivot.
//! let a = Matrix::from_row_major(2, 2, vec![1.0, 1.0, 1.0, 1.0]).unwrap();
//! let fac = Factorization::factor(&a, &Config::new()).unwrap();
//! assert_eq!(fac.rank(), 1);
//! assert!(!fac.is_nonsingular());
//! let err = fac.solve(&[1.0, 1.0]).unwrap_err();
//! assert!(matches!(err, LinSolveError::SingularMatrix { rank: 1 }));
//! ```
//!
//! See the [`cholesky`] module for the factorization handle and the numerical
//! core.
#![forbid(unsafe_code)]
#![warn(missing_docs)]
mod cholesky;
mod config;
mod error;
mod matrix;
pub use cholesky::{factor_in_place, Factorization};
pub use config::{Config, DEFAULT_PIVOT_TOLERANCE, DEFAULT_SYMMETRY_TOLERANCE};
pub use error::LinSolveError;
pub use matrix::Matrix;
/// Crate-level convenience: factor `a` with the default [`Config`] and solve
/// `A x = b` in one call.
///
/// Equivalent to `Factorization::factor(a, &Config::new())?.solve(b)`.
///
/// ```
/// use linsolve::{solve_symmetric, Matrix};
/// // 2x2 SPD system.
/// let a = Matrix::from_row_major(2, 2, vec![2.0, 0.0, 0.0, 4.0]).unwrap();
/// let x = solve_symmetric(&a, &[6.0, 8.0]).unwrap();
/// assert!((x[0] - 3.0).abs() < 1e-12);
/// assert!((x[1] - 2.0).abs() < 1e-12);
/// ```
pub fn solve_symmetric(a: &Matrix, b: &[f64]) -> Result<Vec<f64>, LinSolveError> {
Factorization::factor(a, &Config::new())?.solve(b)
}/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
/workspace/linsolve/src/cholesky.rs
/// The inertia `(positive, negative, zero)` of the signed diagonal `D`.
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")
} /// The inertia `(positive, negative, zero)` of the signed diagonal `D`.
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 {
let mut acc = 0.0;
for k in 0..n {
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 mut out = Matrix::zeros(n, n);
for i in 0..n {
for j in 0..n {
let mut acc = 0.0;
for k in 0..n {
acc += self.lower.get(i, k) * self.diag[k] * self.lower.get(j, k);
}
// `M[i][j] = A[perm[i]][perm[j]]`; scatter back to original space.
out.set(self.perm[i], self.perm[j], acc);
}
}
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(D) = product of the pivots.
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 });
}
// Solve P A Pᵀ (P x) = P b, i.e. L D Lᵀ y = c with c = P b, x = Pᵀ y.
let mut y: Vec<f64> = (0..n).map(|i| rhs[self.perm[i]]).collect();
// Forward substitution: L z = c (L unit lower triangular).
for i in 0..n {
let mut s = y[i];
for j in 0..i {
s -= self.lower.get(i, j) * y[j];
}
y[i] = s;
}
// 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 s = y[i];
for j in (i + 1)..n {
s -= self.lower.get(j, i) * y[j];
}
y[i] = s;
}
// Scatter back to the original index space: x[perm[i]] = y[i].
for i in 0..n {
rhs[self.perm[i]] = y[i];
}
Ok(())
}/workspace/linsolve/src/cholesky.rs
// Numerical core unimplemented in this build.
let _ = (m, lower, diag, perm, pivot_floor);
todo!("implement the signed diagonally-pivoted LDLᵀ factorization in place")
} let n = diag.len();
let mut rank = 0;
for k in 0..n {
// Symmetric pivoting: pick the remaining diagonal of largest magnitude.
// Ties keep the earliest index (strict `>`), matching the reference.
let mut best = k;
for i in (k + 1)..n {
if m.get(i, i).abs() > m.get(best, best).abs() {
best = i;
}
}
let pivot = m.get(best, best);
// Halt if every remaining candidate is negligible (or non-finite): the
// trailing block is the null space and stays as seeded.
if !pivot.is_finite() || pivot.abs() <= pivot_floor {
break;
}
if best != k {
// Symmetric swap of positions `k` and `best`: whole rows/columns of
// the working matrix, the computed columns of `L`, and `perm`.
perm.swap(k, best);
for c in 0..n {
let a = m.get(k, c);
let b = m.get(best, c);
m.set(k, c, b);
m.set(best, c, a);
}
for r in 0..n {
let a = m.get(r, k);
let b = m.get(r, best);
m.set(r, k, b);
m.set(r, best, a);
}
for c in 0..k {
let a = lower.get(k, c);
let b = lower.get(best, c);
lower.set(k, c, b);
lower.set(best, c, a);
}
}
// Eliminate: column `k` of `L`, then update the Schur complement.
diag[k] = pivot;
for i in (k + 1)..n {
lower.set(i, k, m.get(i, k) / pivot);
}
for i in (k + 1)..n {
let lik = lower.get(i, k);
for j in (k + 1)..n {
let updated = m.get(i, j) - lik * m.get(j, k);
m.set(i, j, updated);
}
}
rank += 1;
}
Ok(rank)
}cargo build --release 2>&1 | tail -20
Compiling linsolve v0.4.2 (/workspace/linsolve)
Finished `release` profile [optimized] target(s) in 1.41s
[stdout]
Compiling linsolve v0.4.2 (/workspace/linsolve)
Finished `release` profile [optimized] target(s) in 1.41scat > /tmp/check.rs <<'EOF'
use linsolve::*;
fn main(){
let cfg=Config::new();
let a=Matrix::from_row_major(3,3,vec![4.,12.,-16.,12.,37.,-43.,-16.,-43.,98.]).unwrap();
let f=Factorization::factor(&a,&cfg).unwrap();
println!("SPD rank={} inertia={:?} det={} perm={:?} diag={:?}",f.rank(),f.inertia(),f.signed_determinant(),f.permutation(),f.diagonal());
let r=f.reconstruct_original();
let mut maxe=0f64; for i in 0..3{for j in 0..3{maxe=maxe.max((r.get(i,j)-a.get(i,j)).abs());}}
println!(" recon err={}", maxe);
let a=Matrix::from_row_major(3,3,vec![2.,0.,0.,0.,4.,0.,0.,0.,5.]).unwrap();
let f=Factorization::factor(&a,&cfg).unwrap();
println!("diag245 rank={} inertia={:?} det={} L_is_I={} perm={:?}",f.rank(),f.inertia(),f.signed_determinant(),(0..3).all(|i|(0..3).all(|j|f.lower().get(i,j)==if i==j{1.}else{0.})),f.permutation());
let a=Matrix::from_row_major(3,3,vec![1.,2.,3.,2.,1.,4.,3.,4.,1.]).unwrap();
let f=Factorization::factor(&a,&cfg).unwrap();
let x=f.solve(&[1.,2.,3.]).unwrap();
println!("indef rank={} inertia={:?} det={} x={:?}",f.rank(),f.inertia(),f.signed_determinant(),x);
let a=Matrix::from_row_major(2,2,vec![-1.,0.,0.,-1.]).unwrap();
let f=Factorization::factor(&a,&cfg).unwrap();
println!("negdef rank={} inertia={:?} x={:?}",f.rank(),f.inertia(),f.solve(&[2.,2.]).unwrap());
let a=Matrix::from_row_major(3,3,vec![0.,3.,1.,3.,5.,2.,1.,2.,9.]).unwrap();
let f=Factorization::factor(&a,&cfg).unwrap();
println!("zerolead rank={} inertia={:?}",f.rank(),f.inertia());
let a=Matrix::from_row_major(2,2,vec![1.,1.,1.,1.]).unwrap();
let f=Factorization::factor(&a,&cfg).unwrap();
println!("rank1 rank={} inertia={:?} solveerr={:?}",f.rank(),f.inertia(),f.solve(&[1.,1.]).unwrap_err());
let a=Matrix::from_row_major(2,2,vec![0.,1.,1.,0.]).unwrap();
let f=Factorization::factor(&a,&cfg).unwrap();
println!("antidiag rank={} inertia={:?}",f.rank(),f.inertia());
let a=Matrix::from_row_major(2,2,vec![1.,0.,0.,1e-12]).unwrap();
let f=Factorization::factor(&a,&cfg).unwrap();
let f2=Factorization::factor(&a,&Config::new().with_pivot_tolerance(1e-9)).unwrap();
println!("tiny rank0={} rank9={} inertia9={:?}",f.rank(),f2.rank(),f2.inertia());
}
EOF
mkdir -p examples && cp /tmp/check.rs examples/check.rs && cargo run --release --example check 2>&1 | tail -20; rm examples/check.rsCompiling linsolve v0.4.2 (/workspace/linsolve)
Finished `release` profile [optimized] target(s) in 0.42s
Running `target/release/examples/check`
SPD rank=3 inertia=(3, 0, 0) det=36.000000000000334 perm=[2, 1, 0] diag=[98.0, 18.13265306122449, 0.02025886325267323]
recon err=0.0000000000000017763568394002505
diag245 rank=3 inertia=(3, 0, 0) det=40 L_is_I=true perm=[2, 1, 0]
indef rank=3 inertia=(1, 2, 0) det=20 x=[1.0, -0.0, 0.0]
negdef rank=2 inertia=(0, 2, 0) x=[-2.0, -2.0]
zerolead rank=3 inertia=(2, 1, 0)
rank1 rank=1 inertia=(1, 0, 1) solveerr=SingularMatrix { rank: 1 }
antidiag rank=0 inertia=(0, 0, 2)
tiny rank0=2 rank9=1 inertia9=(1, 0, 1)
[stdout]
Compiling linsolve v0.4.2 (/workspace/linsolve)
Finished `release` profile [optimized] target(s) in 0.42s
Running `target/release/examples/check`
SPD rank=3 inertia=(3, 0, 0) det=36.000000000000334 perm=[2, 1, 0] diag=[98.0, 18.13265306122449, 0.02025886325267323]
recon err=0.0000000000000017763568394002505
diag245 rank=3 inertia=(3, 0, 0) det=40 L_is_I=true perm=[2, 1, 0]
indef rank=3 inertia=(1, 2, 0) det=20 x=[1.0, -0.0, 0.0]
negdef rank=2 inertia=(0, 2, 0) x=[-2.0, -2.0]
zerolead rank=3 inertia=(2, 1, 0)
rank1 rank=1 inertia=(1, 0, 1) solveerr=SingularMatrix { rank: 1 }
antidiag rank=0 inertia=(0, 0, 2)
tiny rank0=2 rank9=1 inertia9=(1, 0, 1)cargo test --release 2>&1 | tail -25
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.41s
[stdout]
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.41scargo clippy --release 2>&1 | tail -5
cat > examples/check.rs <<'EOF'
use linsolve::*;
fn main(){
// 5x5 Hilbert
let n=5; let mut d=vec![]; for i in 0..n{for j in 0..n{d.push(1.0/((i+j+1) as f64));}}
let a=Matrix::from_row_major(n,n,d).unwrap();
let f=Factorization::factor(&a,&Config::new()).unwrap();
let b:Vec<f64>=(0..n).map(|i|(i+1) as f64).collect();
let x=f.solve(&b).unwrap();
let r=a.matvec(&x).unwrap();
let err:f64=r.iter().zip(&b).map(|(ri,bi)|(ri-bi).abs()).fold(0.,f64::max);
println!("Hilbert residual={}",err);
// reconstruct_permuted == P A Pᵀ and symmetric
let perm=f.permutation(); let rp=f.reconstruct_permuted();
let mut e=0f64; for i in 0..n{for j in 0..n{ e=e.max((rp.get(i,j)-a.get(perm[i],perm[j])).abs()); e=e.max((rp.get(i,j)-rp.get(j,i)).abs()); }}
println!("permuted recon+sym err={}",e);
// linearity
let b1=[1.,2.,3.,4.,5.]; let b2=[5.,4.,3.,2.,1.];
let x1=f.solve(&b1).unwrap(); let x2=f.solve(&b2).unwrap();
let bs:Vec<f64>=b1.iter().zip(&b2).map(|(a,b)|a+b).collect();
let xs=f.solve(&bs).unwrap();
let le:f64=(0..n).map(|i|(xs[i]-x1[i]-x2[i]).abs()).fold(0.,f64::max);
println!("linearity err={}",le);
}
EOF
cargo run --release --example check 2>&1 | tail -4; rm examples/check.rs; rmdir examples 2>/dev/null; trueerror: 'cargo-clippy' is not installed for the toolchain '1.83.0-x86_64-unknown-linux-gnu'.
To install, run `rustup component add clippy`
Running `target/release/examples/check`
Hilbert residual=0.0000000000004547473508864641
permuted recon+sym err=0.000000000000000027755575615628914
linearity err=0.00000000004001776687800884
[stdout]
error: 'cargo-clippy' is not installed for the toolchain '1.83.0-x86_64-unknown-linux-gnu'.
To install, run `rustup component add clippy`
Running `target/release/examples/check`
Hilbert residual=0.0000000000004547473508864641
permuted recon+sym err=0.000000000000000027755575615628914
linearity err=0.00000000004001776687800884No 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.45s
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_negative_definite_identity ... ok
test factor_indefinite_3x3_signed_diagonal ... FAILED
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_4e8355b8b20143f1. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_4e8355b8b20143f1 · verifier authoritative; classifier explanatory.