tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it.List source and test files
find /workspace/cgsolve/src -type f | sort && echo "---" && find /workspace/cgsolve/tests -type f | sort
/workspace/cgsolve/src/cg.rs /workspace/cgsolve/src/config.rs /workspace/cgsolve/src/error.rs /workspace/cgsolve/src/lib.rs /workspace/cgsolve/src/matrix.rs /workspace/cgsolve/src/outcome.rs --- /workspace/cgsolve/tests/integration.rs [stdout] /workspace/cgsolve/src/cg.rs /workspace/cgsolve/src/config.rs /workspace/cgsolve/src/error.rs /workspace/cgsolve/src/lib.rs /workspace/cgsolve/src/matrix.rs /workspace/cgsolve/src/outcome.rs --- /workspace/cgsolve/tests/integration.rs
/workspace/cgsolve/src/matrix.rs
1 //! A minimal sparse, square real matrix in compressed-sparse-row (CSR) form.
2 //!
3 //! This is deliberately small: just enough structure to express a sparse
4 //! symmetric coefficient matrix, validate its shape, and perform the single
5 //! operation the conjugate-gradient routine needs , a matrix–vector product.
6 //! Only the nonzero entries are stored, so large sparse systems stay cheap.
7
8 use crate::error::CgError;
9
10 /// A square sparse matrix of `f64` stored in compressed-sparse-row (CSR) order.
11 ///
12 /// Build one with [`SparseMatrix::from_triplets`] (validates the dimension and
13 /// finiteness, and sums duplicate `(row, col)` entries). The matrix is not
14 /// required to be symmetric at construction , symmetry is the caller's
15 /// responsibility for a meaningful conjugate-gradient solve , but a
16 /// [`SparseMatrix::is_symmetric`] check is provided.
17 #[derive(Debug, Clone, PartialEq)]
18 pub struct SparseMatrix {
19 n: usize,
20 /// `row_ptr[i] .. row_ptr[i + 1]` indexes the entries of row `i`.
21 row_ptr: Vec<usize>,
22 /// Column index of each stored entry.
23 col_idx: Vec<usize>,
24 /// Value of each stored entry.
25 values: Vec<f64>,
26 }
27
28 impl SparseMatrix {
29 /// Build an `n × n` sparse matrix from `(row, col, value)` triplets.
30 ///
31 /// Duplicate coordinates are **summed**. Returns [`CgError::EmptyMatrix`] if
32 /// `n == 0`, [`CgError::IndexOutOfBounds`] if any coordinate is `>= n`, and
33 /// [`CgError::NonFiniteEntry`] if any value is `NaN`/infinite.
34 ///
35 /// ```
36 /// use cgsolve::SparseMatrix;
37 /// // The 2x2 identity.
38 /// let a = SparseMatrix::from_triplets(2, &[(0, 0, 1.0), (1, 1, 1.0)]).unwrap();
39 /// assert_eq!(a.dim(), 2);
40 /// assert_eq!(a.nnz(), 2);
41 /// ```
42 pub fn from_triplets(n: usize, triplets: &[(usize, usize, f64)]) -> Result<Self, CgError> {
43 let _ = (n, triplets);
44 todo!("implement from_triplets (sci-4519)")
45 }
46
47 /// Build an `n × n` sparse matrix from a dense row-major buffer, keeping
48 /// only the structurally nonzero entries.
49 ///
50 /// Returns [`CgError::EmptyMatrix`] if `n == 0`,
51 /// [`CgError::DataShapeMismatch`] if `data.len() != n * n`, and
52 /// [`CgError::NonFiniteEntry`] for any non-finite value.
53 pub fn from_dense(n: usize, data: &[f64]) -> Result<Self, CgError> {
54 if n == 0 {
55 return Err(CgError::EmptyMatrix);
56 }
57 if data.len() != n * n {
58 return Err(CgError::DataShapeMismatch { dim: n, len: data.len() });
59 }
60 let mut triplets = Vec::new();
61 for i in 0..n {
62 for j in 0..n {
63 let v = data[i * n + j];
64 if !v.is_finite() {
65 return Err(CgError::NonFiniteEntry { row: i, col: j, value: v });
66 }
67 if v != 0.0 {
68 triplets.push((i, j, v));
69 }
70 }
71 }
72 Self::from_triplets(n, &triplets)
73 }
74
75 /// The system dimension `n`.
76 #[inline]
77 pub fn dim(&self) -> usize {
78 self.n
79 }
80
81 /// The number of stored (structurally nonzero) entries.
82 #[inline]
83 pub fn nnz(&self) -> usize {
84 self.values.len()
85 }
86
87 /// The entry at `(row, col)`, or `0.0` if not stored. Panics if out of
88 /// bounds. `O(nnz in row)`; intended for tests and symmetry checks, not the
89 /// hot loop.
90 pub fn get(&self, row: usize, col: usize) -> f64 {
91 assert!(row < self.n && col < self.n, "index out of bounds");
92 let start = self.row_ptr[row];
93 let end = self.row_ptr[row + 1];
94 for k in start..end {
95 if self.col_idx[k] == col {
96 return self.values[k];
97 }
98 }
99 0.0
100 }
101
102 /// Whether the matrix is symmetric to within absolute tolerance `tol`,
103 /// i.e. `|A[i][j] - A[j][i]| <= tol` for all `i, j`.
104 pub fn is_symmetric(&self, tol: f64) -> bool {
105 for i in 0..self.n {
106 let start = self.row_ptr[i];
107 let end = self.row_ptr[i + 1];
108 for k in start..end {
109 let j = self.col_idx[k];
110 if (self.values[k] - self.get(j, i)).abs() > tol {
111 return false;
112 }
113 }
114 }
115 true
116 }
117
118 /// Compute the matrix–vector product `A * x` into a fresh vector.
119 ///
120 /// Returns [`CgError::DimensionMismatch`] if `x.len() != dim()`.
121 ///
122 /// ```
123 /// use cgsolve::SparseMatrix;
124 /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
125 /// let y = a.matvec(&[1.0, 1.0]).unwrap();
126 /// assert_eq!(y, vec![2.0, 3.0]);
127 /// ```
128 pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, CgError> {
129 if x.len() != self.n {
130 return Err(CgError::DimensionMismatch {
131 expected: self.n,
132 got: x.len(),
133 });
134 }
135 let mut out = vec![0.0; self.n];
136 self.matvec_into(x, &mut out)
137 .expect("output buffer sized to n");
138 Ok(out)
139 }
140
141 /// Compute `A * x`, writing the result into the preallocated `out` buffer.
142 ///
143 /// This is the allocation-free product used inside the conjugate-gradient
144 /// iteration. Returns [`CgError::DimensionMismatch`] if either `x` or `out`
145 /// has the wrong length.
146 pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
147 let _ = (x, out);
148 todo!("implement matvec_into (sci-4519)")
149 }
150 }
151
1 //! The conjugate-gradient method for sparse symmetric-positive-definite systems.
2 //!
3 //! Given a symmetric-positive-definite (SPD) matrix `A` and a right-hand side
4 //! `b`, the conjugate-gradient (CG) method finds `x` solving `A x = b` by a
5 //! sequence of line minimizations of the quadratic `½ xᵀ A x - bᵀ x` along
6 //! mutually `A`-conjugate search directions. Each iteration costs one
7 //! matrix–vector product, so it is the method of choice for large *sparse* SPD
8 //! systems where a direct factorization would fill in.
9 //!
10 //! In exact arithmetic CG converges in at most `n` steps; in floating point
11 //! it is run to a residual tolerance. Each step needs one matrix–vector
12 //! product, and a non-positive curvature `pᵀ A p` signals a non-SPD matrix.
13 //!
14 //! The public entry point is [`ConjugateGradient::solve`]. The numerical core,
15 //! [`cg_iterate`], runs the iteration in place and is invoked once per solve.
16
17 use crate::config::Config;
18 use crate::error::CgError;
19 use crate::matrix::SparseMatrix;
20 use crate::outcome::CgOutcome;
21
22 /// A conjugate-gradient solver bound to a sparse SPD matrix.
23 ///
24 /// Holds a reference to the coefficient matrix `A` and the stopping criteria.
25 /// Reuse a single solver to solve `A x = b` for many right-hand sides.
26 ///
27 /// ```
28 /// use cgsolve::{Config, ConjugateGradient, SparseMatrix};
29 /// // A = [[4, 1], [1, 3]] is SPD.
30 /// let a = SparseMatrix::from_dense(2, &[4.0, 1.0, 1.0, 3.0]).unwrap();
31 /// let cg = ConjugateGradient::new(&a, Config::new());
32 /// // Solve A x = [1, 2].
33 /// let out = cg.solve(&[1.0, 2.0]).unwrap();
34 /// let r = a.matvec(out.solution()).unwrap();
35 /// assert!((r[0] - 1.0).abs() < 1e-9 && (r[1] - 2.0).abs() < 1e-9);
36 /// ```
37 #[derive(Debug, Clone)]
38 pub struct ConjugateGradient<'a> {
39 a: &'a SparseMatrix,
40 config: Config,
41 }
42
43 impl<'a> ConjugateGradient<'a> {
44 /// Create a solver for the matrix `a` with the given [`Config`].
45 pub fn new(a: &'a SparseMatrix, config: Config) -> Self {
46 Self { a, config }
47 }
48
49 /// The system dimension `n`.
50 #[inline]
51 pub fn dim(&self) -> usize {
52 self.a.dim()
53 }
54
55 /// Solve `A x = b`, returning the solution and diagnostics.
56 ///
57 /// The iteration starts from the zero vector. On success the returned
58 /// [`CgOutcome`] carries the solution, the number of iterations, and the
59 /// final residual norm.
60 ///
61 /// # Errors
62 ///
63 /// - [`CgError::InvalidTolerance`] if the configured tolerance is not a
64 /// positive finite number.
65 /// - [`CgError::DimensionMismatch`] if `b.len() != dim()`.
66 /// - [`CgError::NotPositiveDefinite`] if the curvature `pᵀ A p` becomes
67 /// non-positive (the matrix is not SPD).
68 /// - [`CgError::NotConverged`] if the residual tolerance is not reached
69 /// within the iteration budget.
70 pub fn solve(&self, b: &[f64]) -> Result<CgOutcome, CgError> {
71 self.config.validate()?;
72 let n = self.a.dim();
73 if b.len() != n {
74 return Err(CgError::DimensionMismatch {
75 expected: n,
76 got: b.len(),
77 });
78 }
79
80 let max_iterations = self.config.effective_max_iterations(n);
81 let mut x = vec![0.0; n];
82 let (iterations, residual_norm) =
83 cg_iterate(self.a, b, &mut x, self.config.tolerance(), max_iterations)?;
84
85 Ok(CgOutcome {
86 solution: x,
87 iterations,
88 residual_norm,
89 converged: true,
90 })
91 }
92 }
93
94 /// Euclidean dot product of two equal-length slices.
95 #[inline]
96 pub(crate) fn dot(u: &[f64], v: &[f64]) -> f64 {
97 u.iter().zip(v).map(|(a, b)| a * b).sum()
98 }
99
100 /// Numerical core of the conjugate-gradient method.
101 ///
102 /// Solves `A x = b` in place: `x` holds the initial guess on entry (the public
103 /// wrapper passes the zero vector) and the solution on successful return. The
104 /// matrix `A` is assumed symmetric-positive-definite.
105 ///
106 /// Returns `(iterations, residual_norm)` , the number of iterations performed
107 /// and the Euclidean norm of the final residual `b - A x`.
108 ///
109 /// The exact stopping rule, iteration accounting, the non-SPD curvature
110 /// guard, and the error returns are specified at the crate level and pinned
111 /// by `tests/integration.rs`. The module-private `dot` helper is available;
112 /// all work is `O(nnz)` per iteration.
113 pub fn cg_iterate(
114 a: &SparseMatrix,
115 b: &[f64],
116 x: &mut [f64],
117 tolerance: f64,
118 max_iterations: usize,
119 ) -> Result<(usize, f64), CgError> {
120 // TODO(sci-4519): implement the conjugate-gradient iteration described in
121 // the doc comment above. Form the initial residual r = b - A x, iterate the
122 // alpha/x/r/beta/p updates using one matvec per step, stop on the relative
123 // residual threshold, return NotPositiveDefinite on non-positive curvature
124 // and NotConverged if the budget is exhausted. See `tests/integration.rs`
125 // for the contract under test. The `dot` helper below computes the
126 // Euclidean dot products (rᵀr, pᵀAp, ...) you will need.
127 let _ = (a, b, x, tolerance, max_iterations, dot);
128 todo!("implement cg_iterate (sci-4519)")
129 }
130
1 //! Error types for the `cgsolve` crate.
2
3 use thiserror::Error;
4
5 /// Errors that can arise while building a sparse matrix or running a
6 /// conjugate-gradient solve.
7 ///
8 /// These cover the cases where the requested operation cannot be carried out
9 /// reliably (a malformed matrix, a mismatched right-hand side, a breakdown in
10 /// the iteration that indicates the operator is not positive-definite, or
11 /// failure to converge within the iteration budget). A *successful* solve is
12 /// reported through [`CgOutcome`](crate::CgOutcome) instead.
13 #[derive(Debug, Error, Clone, PartialEq)]
14 #[non_exhaustive]
15 pub enum CgError {
16 /// A matrix was constructed with dimension zero. Solving requires at least
17 /// a 1×1 system.
18 #[error("matrix must have dimension at least 1")]
19 EmptyMatrix,
20
21 /// A dense buffer length did not equal `dim * dim`.
22 #[error("dense data length {len} does not match dimension {dim}x{dim}")]
23 DataShapeMismatch {
24 /// The square dimension `n`.
25 dim: usize,
26 /// Length of the supplied data buffer.
27 len: usize,
28 },
29
30 /// A triplet coordinate referenced a row or column outside `0..dim`.
31 #[error("triplet ({row},{col}) is out of bounds for dimension {dim}")]
32 IndexOutOfBounds {
33 /// Offending row index.
34 row: usize,
35 /// Offending column index.
36 col: usize,
37 /// The matrix dimension.
38 dim: usize,
39 },
40
41 /// A supplied matrix entry was not a finite number (it was `NaN` or an
42 /// infinity).
43 #[error("matrix entry at ({row},{col}) is not finite: {value}")]
44 NonFiniteEntry {
45 /// Row index of the offending entry.
46 row: usize,
47 /// Column index of the offending entry.
48 col: usize,
49 /// The non-finite value.
50 value: f64,
51 },
52
53 /// The right-hand side (or an operand) had a length that did not match the
54 /// system dimension.
55 #[error("dimension mismatch: expected length {expected}, got {got}")]
56 DimensionMismatch {
57 /// The dimension required.
58 expected: usize,
59 /// The length actually supplied.
60 got: usize,
61 },
62
63 /// The configured tolerance was not a strictly positive, finite number.
64 #[error("tolerance must be finite and strictly positive, got {0}")]
65 InvalidTolerance(f64),
66
67 /// The conjugate-gradient iteration broke down because a curvature term
68 /// `pᵀ A p` was non-positive (or non-finite). For a genuinely
69 /// symmetric-positive-definite operator this cannot happen; it signals that
70 /// the matrix is indefinite or not positive-definite. Reports the iteration
71 /// at which the breakdown occurred and the offending curvature value.
72 #[error("conjugate-gradient breakdown at iteration {iteration}: pᵀ A p = {curvature} is not positive (matrix not SPD?)")]
73 NotPositiveDefinite {
74 /// The iteration index at which the breakdown occurred.
75 iteration: usize,
76 /// The offending curvature value `pᵀ A p`.
77 curvature: f64,
78 },
79
80 /// The iteration did not reach the requested residual tolerance within the
81 /// allotted number of iterations. Reports the best residual norm achieved.
82 #[error("failed to converge within {max_iterations} iterations (residual norm {residual_norm})")]
83 NotConverged {
84 /// The iteration budget that was exhausted.
85 max_iterations: usize,
86 /// The residual norm at the final iterate.
87 residual_norm: f64,
88 },
89 }
90
/workspace/cgsolve/src/cg.rs
1 //! A minimal sparse, square real matrix in compressed-sparse-row (CSR) form.
2 //!
3 //! This is deliberately small: just enough structure to express a sparse
4 //! symmetric coefficient matrix, validate its shape, and perform the single
5 //! operation the conjugate-gradient routine needs , a matrix–vector product.
6 //! Only the nonzero entries are stored, so large sparse systems stay cheap.
7
8 use crate::error::CgError;
9
10 /// A square sparse matrix of `f64` stored in compressed-sparse-row (CSR) order.
11 ///
12 /// Build one with [`SparseMatrix::from_triplets`] (validates the dimension and
13 /// finiteness, and sums duplicate `(row, col)` entries). The matrix is not
14 /// required to be symmetric at construction , symmetry is the caller's
15 /// responsibility for a meaningful conjugate-gradient solve , but a
16 /// [`SparseMatrix::is_symmetric`] check is provided.
17 #[derive(Debug, Clone, PartialEq)]
18 pub struct SparseMatrix {
19 n: usize,
20 /// `row_ptr[i] .. row_ptr[i + 1]` indexes the entries of row `i`.
21 row_ptr: Vec<usize>,
22 /// Column index of each stored entry.
23 col_idx: Vec<usize>,
24 /// Value of each stored entry.
25 values: Vec<f64>,
26 }
27
28 impl SparseMatrix {
29 /// Build an `n × n` sparse matrix from `(row, col, value)` triplets.
30 ///
31 /// Duplicate coordinates are **summed**. Returns [`CgError::EmptyMatrix`] if
32 /// `n == 0`, [`CgError::IndexOutOfBounds`] if any coordinate is `>= n`, and
33 /// [`CgError::NonFiniteEntry`] if any value is `NaN`/infinite.
34 ///
35 /// ```
36 /// use cgsolve::SparseMatrix;
37 /// // The 2x2 identity.
38 /// let a = SparseMatrix::from_triplets(2, &[(0, 0, 1.0), (1, 1, 1.0)]).unwrap();
39 /// assert_eq!(a.dim(), 2);
40 /// assert_eq!(a.nnz(), 2);
41 /// ```
42 pub fn from_triplets(n: usize, triplets: &[(usize, usize, f64)]) -> Result<Self, CgError> {
43 let _ = (n, triplets);
44 todo!("implement from_triplets (sci-4519)")
45 }
46
47 /// Build an `n × n` sparse matrix from a dense row-major buffer, keeping
48 /// only the structurally nonzero entries.
49 ///
50 /// Returns [`CgError::EmptyMatrix`] if `n == 0`,
51 /// [`CgError::DataShapeMismatch`] if `data.len() != n * n`, and
52 /// [`CgError::NonFiniteEntry`] for any non-finite value.
53 pub fn from_dense(n: usize, data: &[f64]) -> Result<Self, CgError> {
54 if n == 0 {
55 return Err(CgError::EmptyMatrix);
56 }
57 if data.len() != n * n {
58 return Err(CgError::DataShapeMismatch { dim: n, len: data.len() });
59 }
60 let mut triplets = Vec::new();
61 for i in 0..n {
62 for j in 0..n {
63 let v = data[i * n + j];
64 if !v.is_finite() {
65 return Err(CgError::NonFiniteEntry { row: i, col: j, value: v });
66 }
67 if v != 0.0 {
68 triplets.push((i, j, v));
69 }
70 }
71 }
72 Self::from_triplets(n, &triplets)
73 }
74
75 /// The system dimension `n`.
76 #[inline]
77 pub fn dim(&self) -> usize {
78 self.n
79 }
80
81 /// The number of stored (structurally nonzero) entries.
82 #[inline]
83 pub fn nnz(&self) -> usize {
84 self.values.len()
85 }
86
87 /// The entry at `(row, col)`, or `0.0` if not stored. Panics if out of
88 /// bounds. `O(nnz in row)`; intended for tests and symmetry checks, not the
89 /// hot loop.
90 pub fn get(&self, row: usize, col: usize) -> f64 {
91 assert!(row < self.n && col < self.n, "index out of bounds");
92 let start = self.row_ptr[row];
93 let end = self.row_ptr[row + 1];
94 for k in start..end {
95 if self.col_idx[k] == col {
96 return self.values[k];
97 }
98 }
99 0.0
100 }
101
102 /// Whether the matrix is symmetric to within absolute tolerance `tol`,
103 /// i.e. `|A[i][j] - A[j][i]| <= tol` for all `i, j`.
104 pub fn is_symmetric(&self, tol: f64) -> bool {
105 for i in 0..self.n {
106 let start = self.row_ptr[i];
107 let end = self.row_ptr[i + 1];
108 for k in start..end {
109 let j = self.col_idx[k];
110 if (self.values[k] - self.get(j, i)).abs() > tol {
111 return false;
112 }
113 }
114 }
115 true
116 }
117
118 /// Compute the matrix–vector product `A * x` into a fresh vector.
119 ///
120 /// Returns [`CgError::DimensionMismatch`] if `x.len() != dim()`.
121 ///
122 /// ```
123 /// use cgsolve::SparseMatrix;
124 /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
125 /// let y = a.matvec(&[1.0, 1.0]).unwrap();
126 /// assert_eq!(y, vec![2.0, 3.0]);
127 /// ```
128 pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, CgError> {
129 if x.len() != self.n {
130 return Err(CgError::DimensionMismatch {
131 expected: self.n,
132 got: x.len(),
133 });
134 }
135 let mut out = vec![0.0; self.n];
136 self.matvec_into(x, &mut out)
137 .expect("output buffer sized to n");
138 Ok(out)
139 }
140
141 /// Compute `A * x`, writing the result into the preallocated `out` buffer.
142 ///
143 /// This is the allocation-free product used inside the conjugate-gradient
144 /// iteration. Returns [`CgError::DimensionMismatch`] if either `x` or `out`
145 /// has the wrong length.
146 pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
147 let _ = (x, out);
148 todo!("implement matvec_into (sci-4519)")
149 }
150 }
151
1 //! The conjugate-gradient method for sparse symmetric-positive-definite systems.
2 //!
3 //! Given a symmetric-positive-definite (SPD) matrix `A` and a right-hand side
4 //! `b`, the conjugate-gradient (CG) method finds `x` solving `A x = b` by a
5 //! sequence of line minimizations of the quadratic `½ xᵀ A x - bᵀ x` along
6 //! mutually `A`-conjugate search directions. Each iteration costs one
7 //! matrix–vector product, so it is the method of choice for large *sparse* SPD
8 //! systems where a direct factorization would fill in.
9 //!
10 //! In exact arithmetic CG converges in at most `n` steps; in floating point
11 //! it is run to a residual tolerance. Each step needs one matrix–vector
12 //! product, and a non-positive curvature `pᵀ A p` signals a non-SPD matrix.
13 //!
14 //! The public entry point is [`ConjugateGradient::solve`]. The numerical core,
15 //! [`cg_iterate`], runs the iteration in place and is invoked once per solve.
16
17 use crate::config::Config;
18 use crate::error::CgError;
19 use crate::matrix::SparseMatrix;
20 use crate::outcome::CgOutcome;
21
22 /// A conjugate-gradient solver bound to a sparse SPD matrix.
23 ///
24 /// Holds a reference to the coefficient matrix `A` and the stopping criteria.
25 /// Reuse a single solver to solve `A x = b` for many right-hand sides.
26 ///
27 /// ```
28 /// use cgsolve::{Config, ConjugateGradient, SparseMatrix};
29 /// // A = [[4, 1], [1, 3]] is SPD.
30 /// let a = SparseMatrix::from_dense(2, &[4.0, 1.0, 1.0, 3.0]).unwrap();
31 /// let cg = ConjugateGradient::new(&a, Config::new());
32 /// // Solve A x = [1, 2].
33 /// let out = cg.solve(&[1.0, 2.0]).unwrap();
34 /// let r = a.matvec(out.solution()).unwrap();
35 /// assert!((r[0] - 1.0).abs() < 1e-9 && (r[1] - 2.0).abs() < 1e-9);
36 /// ```
37 #[derive(Debug, Clone)]
38 pub struct ConjugateGradient<'a> {
39 a: &'a SparseMatrix,
40 config: Config,
41 }
42
43 impl<'a> ConjugateGradient<'a> {
44 /// Create a solver for the matrix `a` with the given [`Config`].
45 pub fn new(a: &'a SparseMatrix, config: Config) -> Self {
46 Self { a, config }
47 }
48
49 /// The system dimension `n`.
50 #[inline]
51 pub fn dim(&self) -> usize {
52 self.a.dim()
53 }
54
55 /// Solve `A x = b`, returning the solution and diagnostics.
56 ///
57 /// The iteration starts from the zero vector. On success the returned
58 /// [`CgOutcome`] carries the solution, the number of iterations, and the
59 /// final residual norm.
60 ///
61 /// # Errors
62 ///
63 /// - [`CgError::InvalidTolerance`] if the configured tolerance is not a
64 /// positive finite number.
65 /// - [`CgError::DimensionMismatch`] if `b.len() != dim()`.
66 /// - [`CgError::NotPositiveDefinite`] if the curvature `pᵀ A p` becomes
67 /// non-positive (the matrix is not SPD).
68 /// - [`CgError::NotConverged`] if the residual tolerance is not reached
69 /// within the iteration budget.
70 pub fn solve(&self, b: &[f64]) -> Result<CgOutcome, CgError> {
71 self.config.validate()?;
72 let n = self.a.dim();
73 if b.len() != n {
74 return Err(CgError::DimensionMismatch {
75 expected: n,
76 got: b.len(),
77 });
78 }
79
80 let max_iterations = self.config.effective_max_iterations(n);
81 let mut x = vec![0.0; n];
82 let (iterations, residual_norm) =
83 cg_iterate(self.a, b, &mut x, self.config.tolerance(), max_iterations)?;
84
85 Ok(CgOutcome {
86 solution: x,
87 iterations,
88 residual_norm,
89 converged: true,
90 })
91 }
92 }
93
94 /// Euclidean dot product of two equal-length slices.
95 #[inline]
96 pub(crate) fn dot(u: &[f64], v: &[f64]) -> f64 {
97 u.iter().zip(v).map(|(a, b)| a * b).sum()
98 }
99
100 /// Numerical core of the conjugate-gradient method.
101 ///
102 /// Solves `A x = b` in place: `x` holds the initial guess on entry (the public
103 /// wrapper passes the zero vector) and the solution on successful return. The
104 /// matrix `A` is assumed symmetric-positive-definite.
105 ///
106 /// Returns `(iterations, residual_norm)` , the number of iterations performed
107 /// and the Euclidean norm of the final residual `b - A x`.
108 ///
109 /// The exact stopping rule, iteration accounting, the non-SPD curvature
110 /// guard, and the error returns are specified at the crate level and pinned
111 /// by `tests/integration.rs`. The module-private `dot` helper is available;
112 /// all work is `O(nnz)` per iteration.
113 pub fn cg_iterate(
114 a: &SparseMatrix,
115 b: &[f64],
116 x: &mut [f64],
117 tolerance: f64,
118 max_iterations: usize,
119 ) -> Result<(usize, f64), CgError> {
120 // TODO(sci-4519): implement the conjugate-gradient iteration described in
121 // the doc comment above. Form the initial residual r = b - A x, iterate the
122 // alpha/x/r/beta/p updates using one matvec per step, stop on the relative
123 // residual threshold, return NotPositiveDefinite on non-positive curvature
124 // and NotConverged if the budget is exhausted. See `tests/integration.rs`
125 // for the contract under test. The `dot` helper below computes the
126 // Euclidean dot products (rᵀr, pᵀAp, ...) you will need.
127 let _ = (a, b, x, tolerance, max_iterations, dot);
128 todo!("implement cg_iterate (sci-4519)")
129 }
130
1 //! Error types for the `cgsolve` crate.
2
3 use thiserror::Error;
4
5 /// Errors that can arise while building a sparse matrix or running a
6 /// conjugate-gradient solve.
7 ///
8 /// These cover the cases where the requested operation cannot be carried out
9 /// reliably (a malformed matrix, a mismatched right-hand side, a breakdown in
10 /// the iteration that indicates the operator is not positive-definite, or
11 /// failure to converge within the iteration budget). A *successful* solve is
12 /// reported through [`CgOutcome`](crate::CgOutcome) instead.
13 #[derive(Debug, Error, Clone, PartialEq)]
14 #[non_exhaustive]
15 pub enum CgError {
16 /// A matrix was constructed with dimension zero. Solving requires at least
17 /// a 1×1 system.
18 #[error("matrix must have dimension at least 1")]
19 EmptyMatrix,
20
21 /// A dense buffer length did not equal `dim * dim`.
22 #[error("dense data length {len} does not match dimension {dim}x{dim}")]
23 DataShapeMismatch {
24 /// The square dimension `n`.
25 dim: usize,
26 /// Length of the supplied data buffer.
27 len: usize,
28 },
29
30 /// A triplet coordinate referenced a row or column outside `0..dim`.
31 #[error("triplet ({row},{col}) is out of bounds for dimension {dim}")]
32 IndexOutOfBounds {
33 /// Offending row index.
34 row: usize,
35 /// Offending column index.
36 col: usize,
37 /// The matrix dimension.
38 dim: usize,
39 },
40
41 /// A supplied matrix entry was not a finite number (it was `NaN` or an
42 /// infinity).
43 #[error("matrix entry at ({row},{col}) is not finite: {value}")]
44 NonFiniteEntry {
45 /// Row index of the offending entry.
46 row: usize,
47 /// Column index of the offending entry.
48 col: usize,
49 /// The non-finite value.
50 value: f64,
51 },
52
53 /// The right-hand side (or an operand) had a length that did not match the
54 /// system dimension.
55 #[error("dimension mismatch: expected length {expected}, got {got}")]
56 DimensionMismatch {
57 /// The dimension required.
58 expected: usize,
59 /// The length actually supplied.
60 got: usize,
61 },
62
63 /// The configured tolerance was not a strictly positive, finite number.
64 #[error("tolerance must be finite and strictly positive, got {0}")]
65 InvalidTolerance(f64),
66
67 /// The conjugate-gradient iteration broke down because a curvature term
68 /// `pᵀ A p` was non-positive (or non-finite). For a genuinely
69 /// symmetric-positive-definite operator this cannot happen; it signals that
70 /// the matrix is indefinite or not positive-definite. Reports the iteration
71 /// at which the breakdown occurred and the offending curvature value.
72 #[error("conjugate-gradient breakdown at iteration {iteration}: pᵀ A p = {curvature} is not positive (matrix not SPD?)")]
73 NotPositiveDefinite {
74 /// The iteration index at which the breakdown occurred.
75 iteration: usize,
76 /// The offending curvature value `pᵀ A p`.
77 curvature: f64,
78 },
79
80 /// The iteration did not reach the requested residual tolerance within the
81 /// allotted number of iterations. Reports the best residual norm achieved.
82 #[error("failed to converge within {max_iterations} iterations (residual norm {residual_norm})")]
83 NotConverged {
84 /// The iteration budget that was exhausted.
85 max_iterations: usize,
86 /// The residual norm at the final iterate.
87 residual_norm: f64,
88 },
89 }
90
/workspace/cgsolve/src/error.rs
1 //! A minimal sparse, square real matrix in compressed-sparse-row (CSR) form.
2 //!
3 //! This is deliberately small: just enough structure to express a sparse
4 //! symmetric coefficient matrix, validate its shape, and perform the single
5 //! operation the conjugate-gradient routine needs , a matrix–vector product.
6 //! Only the nonzero entries are stored, so large sparse systems stay cheap.
7
8 use crate::error::CgError;
9
10 /// A square sparse matrix of `f64` stored in compressed-sparse-row (CSR) order.
11 ///
12 /// Build one with [`SparseMatrix::from_triplets`] (validates the dimension and
13 /// finiteness, and sums duplicate `(row, col)` entries). The matrix is not
14 /// required to be symmetric at construction , symmetry is the caller's
15 /// responsibility for a meaningful conjugate-gradient solve , but a
16 /// [`SparseMatrix::is_symmetric`] check is provided.
17 #[derive(Debug, Clone, PartialEq)]
18 pub struct SparseMatrix {
19 n: usize,
20 /// `row_ptr[i] .. row_ptr[i + 1]` indexes the entries of row `i`.
21 row_ptr: Vec<usize>,
22 /// Column index of each stored entry.
23 col_idx: Vec<usize>,
24 /// Value of each stored entry.
25 values: Vec<f64>,
26 }
27
28 impl SparseMatrix {
29 /// Build an `n × n` sparse matrix from `(row, col, value)` triplets.
30 ///
31 /// Duplicate coordinates are **summed**. Returns [`CgError::EmptyMatrix`] if
32 /// `n == 0`, [`CgError::IndexOutOfBounds`] if any coordinate is `>= n`, and
33 /// [`CgError::NonFiniteEntry`] if any value is `NaN`/infinite.
34 ///
35 /// ```
36 /// use cgsolve::SparseMatrix;
37 /// // The 2x2 identity.
38 /// let a = SparseMatrix::from_triplets(2, &[(0, 0, 1.0), (1, 1, 1.0)]).unwrap();
39 /// assert_eq!(a.dim(), 2);
40 /// assert_eq!(a.nnz(), 2);
41 /// ```
42 pub fn from_triplets(n: usize, triplets: &[(usize, usize, f64)]) -> Result<Self, CgError> {
43 let _ = (n, triplets);
44 todo!("implement from_triplets (sci-4519)")
45 }
46
47 /// Build an `n × n` sparse matrix from a dense row-major buffer, keeping
48 /// only the structurally nonzero entries.
49 ///
50 /// Returns [`CgError::EmptyMatrix`] if `n == 0`,
51 /// [`CgError::DataShapeMismatch`] if `data.len() != n * n`, and
52 /// [`CgError::NonFiniteEntry`] for any non-finite value.
53 pub fn from_dense(n: usize, data: &[f64]) -> Result<Self, CgError> {
54 if n == 0 {
55 return Err(CgError::EmptyMatrix);
56 }
57 if data.len() != n * n {
58 return Err(CgError::DataShapeMismatch { dim: n, len: data.len() });
59 }
60 let mut triplets = Vec::new();
61 for i in 0..n {
62 for j in 0..n {
63 let v = data[i * n + j];
64 if !v.is_finite() {
65 return Err(CgError::NonFiniteEntry { row: i, col: j, value: v });
66 }
67 if v != 0.0 {
68 triplets.push((i, j, v));
69 }
70 }
71 }
72 Self::from_triplets(n, &triplets)
73 }
74
75 /// The system dimension `n`.
76 #[inline]
77 pub fn dim(&self) -> usize {
78 self.n
79 }
80
81 /// The number of stored (structurally nonzero) entries.
82 #[inline]
83 pub fn nnz(&self) -> usize {
84 self.values.len()
85 }
86
87 /// The entry at `(row, col)`, or `0.0` if not stored. Panics if out of
88 /// bounds. `O(nnz in row)`; intended for tests and symmetry checks, not the
89 /// hot loop.
90 pub fn get(&self, row: usize, col: usize) -> f64 {
91 assert!(row < self.n && col < self.n, "index out of bounds");
92 let start = self.row_ptr[row];
93 let end = self.row_ptr[row + 1];
94 for k in start..end {
95 if self.col_idx[k] == col {
96 return self.values[k];
97 }
98 }
99 0.0
100 }
101
102 /// Whether the matrix is symmetric to within absolute tolerance `tol`,
103 /// i.e. `|A[i][j] - A[j][i]| <= tol` for all `i, j`.
104 pub fn is_symmetric(&self, tol: f64) -> bool {
105 for i in 0..self.n {
106 let start = self.row_ptr[i];
107 let end = self.row_ptr[i + 1];
108 for k in start..end {
109 let j = self.col_idx[k];
110 if (self.values[k] - self.get(j, i)).abs() > tol {
111 return false;
112 }
113 }
114 }
115 true
116 }
117
118 /// Compute the matrix–vector product `A * x` into a fresh vector.
119 ///
120 /// Returns [`CgError::DimensionMismatch`] if `x.len() != dim()`.
121 ///
122 /// ```
123 /// use cgsolve::SparseMatrix;
124 /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
125 /// let y = a.matvec(&[1.0, 1.0]).unwrap();
126 /// assert_eq!(y, vec![2.0, 3.0]);
127 /// ```
128 pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, CgError> {
129 if x.len() != self.n {
130 return Err(CgError::DimensionMismatch {
131 expected: self.n,
132 got: x.len(),
133 });
134 }
135 let mut out = vec![0.0; self.n];
136 self.matvec_into(x, &mut out)
137 .expect("output buffer sized to n");
138 Ok(out)
139 }
140
141 /// Compute `A * x`, writing the result into the preallocated `out` buffer.
142 ///
143 /// This is the allocation-free product used inside the conjugate-gradient
144 /// iteration. Returns [`CgError::DimensionMismatch`] if either `x` or `out`
145 /// has the wrong length.
146 pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
147 let _ = (x, out);
148 todo!("implement matvec_into (sci-4519)")
149 }
150 }
151
1 //! The conjugate-gradient method for sparse symmetric-positive-definite systems.
2 //!
3 //! Given a symmetric-positive-definite (SPD) matrix `A` and a right-hand side
4 //! `b`, the conjugate-gradient (CG) method finds `x` solving `A x = b` by a
5 //! sequence of line minimizations of the quadratic `½ xᵀ A x - bᵀ x` along
6 //! mutually `A`-conjugate search directions. Each iteration costs one
7 //! matrix–vector product, so it is the method of choice for large *sparse* SPD
8 //! systems where a direct factorization would fill in.
9 //!
10 //! In exact arithmetic CG converges in at most `n` steps; in floating point
11 //! it is run to a residual tolerance. Each step needs one matrix–vector
12 //! product, and a non-positive curvature `pᵀ A p` signals a non-SPD matrix.
13 //!
14 //! The public entry point is [`ConjugateGradient::solve`]. The numerical core,
15 //! [`cg_iterate`], runs the iteration in place and is invoked once per solve.
16
17 use crate::config::Config;
18 use crate::error::CgError;
19 use crate::matrix::SparseMatrix;
20 use crate::outcome::CgOutcome;
21
22 /// A conjugate-gradient solver bound to a sparse SPD matrix.
23 ///
24 /// Holds a reference to the coefficient matrix `A` and the stopping criteria.
25 /// Reuse a single solver to solve `A x = b` for many right-hand sides.
26 ///
27 /// ```
28 /// use cgsolve::{Config, ConjugateGradient, SparseMatrix};
29 /// // A = [[4, 1], [1, 3]] is SPD.
30 /// let a = SparseMatrix::from_dense(2, &[4.0, 1.0, 1.0, 3.0]).unwrap();
31 /// let cg = ConjugateGradient::new(&a, Config::new());
32 /// // Solve A x = [1, 2].
33 /// let out = cg.solve(&[1.0, 2.0]).unwrap();
34 /// let r = a.matvec(out.solution()).unwrap();
35 /// assert!((r[0] - 1.0).abs() < 1e-9 && (r[1] - 2.0).abs() < 1e-9);
36 /// ```
37 #[derive(Debug, Clone)]
38 pub struct ConjugateGradient<'a> {
39 a: &'a SparseMatrix,
40 config: Config,
41 }
42
43 impl<'a> ConjugateGradient<'a> {
44 /// Create a solver for the matrix `a` with the given [`Config`].
45 pub fn new(a: &'a SparseMatrix, config: Config) -> Self {
46 Self { a, config }
47 }
48
49 /// The system dimension `n`.
50 #[inline]
51 pub fn dim(&self) -> usize {
52 self.a.dim()
53 }
54
55 /// Solve `A x = b`, returning the solution and diagnostics.
56 ///
57 /// The iteration starts from the zero vector. On success the returned
58 /// [`CgOutcome`] carries the solution, the number of iterations, and the
59 /// final residual norm.
60 ///
61 /// # Errors
62 ///
63 /// - [`CgError::InvalidTolerance`] if the configured tolerance is not a
64 /// positive finite number.
65 /// - [`CgError::DimensionMismatch`] if `b.len() != dim()`.
66 /// - [`CgError::NotPositiveDefinite`] if the curvature `pᵀ A p` becomes
67 /// non-positive (the matrix is not SPD).
68 /// - [`CgError::NotConverged`] if the residual tolerance is not reached
69 /// within the iteration budget.
70 pub fn solve(&self, b: &[f64]) -> Result<CgOutcome, CgError> {
71 self.config.validate()?;
72 let n = self.a.dim();
73 if b.len() != n {
74 return Err(CgError::DimensionMismatch {
75 expected: n,
76 got: b.len(),
77 });
78 }
79
80 let max_iterations = self.config.effective_max_iterations(n);
81 let mut x = vec![0.0; n];
82 let (iterations, residual_norm) =
83 cg_iterate(self.a, b, &mut x, self.config.tolerance(), max_iterations)?;
84
85 Ok(CgOutcome {
86 solution: x,
87 iterations,
88 residual_norm,
89 converged: true,
90 })
91 }
92 }
93
94 /// Euclidean dot product of two equal-length slices.
95 #[inline]
96 pub(crate) fn dot(u: &[f64], v: &[f64]) -> f64 {
97 u.iter().zip(v).map(|(a, b)| a * b).sum()
98 }
99
100 /// Numerical core of the conjugate-gradient method.
101 ///
102 /// Solves `A x = b` in place: `x` holds the initial guess on entry (the public
103 /// wrapper passes the zero vector) and the solution on successful return. The
104 /// matrix `A` is assumed symmetric-positive-definite.
105 ///
106 /// Returns `(iterations, residual_norm)` , the number of iterations performed
107 /// and the Euclidean norm of the final residual `b - A x`.
108 ///
109 /// The exact stopping rule, iteration accounting, the non-SPD curvature
110 /// guard, and the error returns are specified at the crate level and pinned
111 /// by `tests/integration.rs`. The module-private `dot` helper is available;
112 /// all work is `O(nnz)` per iteration.
113 pub fn cg_iterate(
114 a: &SparseMatrix,
115 b: &[f64],
116 x: &mut [f64],
117 tolerance: f64,
118 max_iterations: usize,
119 ) -> Result<(usize, f64), CgError> {
120 // TODO(sci-4519): implement the conjugate-gradient iteration described in
121 // the doc comment above. Form the initial residual r = b - A x, iterate the
122 // alpha/x/r/beta/p updates using one matvec per step, stop on the relative
123 // residual threshold, return NotPositiveDefinite on non-positive curvature
124 // and NotConverged if the budget is exhausted. See `tests/integration.rs`
125 // for the contract under test. The `dot` helper below computes the
126 // Euclidean dot products (rᵀr, pᵀAp, ...) you will need.
127 let _ = (a, b, x, tolerance, max_iterations, dot);
128 todo!("implement cg_iterate (sci-4519)")
129 }
130
1 //! Error types for the `cgsolve` crate.
2
3 use thiserror::Error;
4
5 /// Errors that can arise while building a sparse matrix or running a
6 /// conjugate-gradient solve.
7 ///
8 /// These cover the cases where the requested operation cannot be carried out
9 /// reliably (a malformed matrix, a mismatched right-hand side, a breakdown in
10 /// the iteration that indicates the operator is not positive-definite, or
11 /// failure to converge within the iteration budget). A *successful* solve is
12 /// reported through [`CgOutcome`](crate::CgOutcome) instead.
13 #[derive(Debug, Error, Clone, PartialEq)]
14 #[non_exhaustive]
15 pub enum CgError {
16 /// A matrix was constructed with dimension zero. Solving requires at least
17 /// a 1×1 system.
18 #[error("matrix must have dimension at least 1")]
19 EmptyMatrix,
20
21 /// A dense buffer length did not equal `dim * dim`.
22 #[error("dense data length {len} does not match dimension {dim}x{dim}")]
23 DataShapeMismatch {
24 /// The square dimension `n`.
25 dim: usize,
26 /// Length of the supplied data buffer.
27 len: usize,
28 },
29
30 /// A triplet coordinate referenced a row or column outside `0..dim`.
31 #[error("triplet ({row},{col}) is out of bounds for dimension {dim}")]
32 IndexOutOfBounds {
33 /// Offending row index.
34 row: usize,
35 /// Offending column index.
36 col: usize,
37 /// The matrix dimension.
38 dim: usize,
39 },
40
41 /// A supplied matrix entry was not a finite number (it was `NaN` or an
42 /// infinity).
43 #[error("matrix entry at ({row},{col}) is not finite: {value}")]
44 NonFiniteEntry {
45 /// Row index of the offending entry.
46 row: usize,
47 /// Column index of the offending entry.
48 col: usize,
49 /// The non-finite value.
50 value: f64,
51 },
52
53 /// The right-hand side (or an operand) had a length that did not match the
54 /// system dimension.
55 #[error("dimension mismatch: expected length {expected}, got {got}")]
56 DimensionMismatch {
57 /// The dimension required.
58 expected: usize,
59 /// The length actually supplied.
60 got: usize,
61 },
62
63 /// The configured tolerance was not a strictly positive, finite number.
64 #[error("tolerance must be finite and strictly positive, got {0}")]
65 InvalidTolerance(f64),
66
67 /// The conjugate-gradient iteration broke down because a curvature term
68 /// `pᵀ A p` was non-positive (or non-finite). For a genuinely
69 /// symmetric-positive-definite operator this cannot happen; it signals that
70 /// the matrix is indefinite or not positive-definite. Reports the iteration
71 /// at which the breakdown occurred and the offending curvature value.
72 #[error("conjugate-gradient breakdown at iteration {iteration}: pᵀ A p = {curvature} is not positive (matrix not SPD?)")]
73 NotPositiveDefinite {
74 /// The iteration index at which the breakdown occurred.
75 iteration: usize,
76 /// The offending curvature value `pᵀ A p`.
77 curvature: f64,
78 },
79
80 /// The iteration did not reach the requested residual tolerance within the
81 /// allotted number of iterations. Reports the best residual norm achieved.
82 #[error("failed to converge within {max_iterations} iterations (residual norm {residual_norm})")]
83 NotConverged {
84 /// The iteration budget that was exhausted.
85 max_iterations: usize,
86 /// The residual norm at the final iterate.
87 residual_norm: f64,
88 },
89 }
90
/workspace/cgsolve/tests/integration.rs
1 //! Integration tests for the `cgsolve` conjugate-gradient solver.
2 //!
3 //! These exercise the public API against systems with known solutions, a large
4 //! sparse system, edge cases, every error path, and three structural
5 //! properties , the residual `b - A x ≈ 0`, the finite-termination convergence
6 //! guarantee (CG converges within `n` iterations), and linearity in the
7 //! right-hand side. A correct CG implementation passes all of them; a hardcoded
8 //! or trivially constant implementation does not.
9
10 use approx::assert_relative_eq;
11 use cgsolve::{solve_spd, CgError, Config, ConjugateGradient, SparseMatrix};
12
13 /// Build a dense-specified SPD matrix, panicking on malformed input.
14 fn dense(n: usize, data: Vec<f64>) -> SparseMatrix {
15 SparseMatrix::from_dense(n, &data).expect("valid matrix literal")
16 }
17
18 /// Euclidean norm of a vector.
19 fn norm(v: &[f64]) -> f64 {
20 v.iter().map(|x| x * x).sum::<f64>().sqrt()
21 }
22
23 /// Residual norm ‖A x - b‖ for a candidate solution `x`.
24 fn residual_norm(a: &SparseMatrix, x: &[f64], b: &[f64]) -> f64 {
25 let ax = a.matvec(x).expect("dimensions match");
26 let r: Vec<f64> = ax.iter().zip(b).map(|(p, q)| p - q).collect();
27 norm(&r)
28 }
29
30 /// Build the 1-D discrete Laplacian `tridiag(-1, 2, -1)` of dimension `n`.
31 /// This is a classic large sparse SPD operator (the matrix of the second
32 /// difference), positive-definite with eigenvalues `2 - 2 cos(kπ/(n+1))`.
33 fn laplacian(n: usize) -> SparseMatrix {
34 let mut triplets = Vec::new();
35 for i in 0..n {
36 triplets.push((i, i, 2.0));
37 if i + 1 < n {
38 triplets.push((i, i + 1, -1.0));
39 triplets.push((i + 1, i, -1.0));
40 }
41 }
42 SparseMatrix::from_triplets(n, &triplets).expect("valid Laplacian")
43 }
44
45 // ---------------------------------------------------------------------------
46 // Known systems with hand-checked solutions
47 // ---------------------------------------------------------------------------
48
49 #[test]
50 fn solves_diagonal_system() {
51 // diag(2,4,5) x = (6,8,10) => x = (3,2,2).
52 let a = dense(3, vec![2.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 5.0]);
53 let out = solve_spd(&a, &[6.0, 8.0, 10.0]).unwrap();
54 assert_relative_eq!(out.solution[0], 3.0, epsilon = 1e-9);
55 assert_relative_eq!(out.solution[1], 2.0, epsilon = 1e-9);
56 assert_relative_eq!(out.solution[2], 2.0, epsilon = 1e-9);
57 }
58
59 #[test]
60 fn solves_known_2x2_system() {
61 // A = [[4, 1], [1, 3]], b = [1, 2]. Closed form: det = 11,
62 // x0 = (1*3 - 1*2)/11 = 1/11, x1 = (4*2 - 1*1)/11 = 7/11.
63 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
64 let out = solve_spd(&a, &[1.0, 2.0]).unwrap();
65 assert_relative_eq!(out.solution[0], 1.0 / 11.0, epsilon = 1e-9);
66 assert_relative_eq!(out.solution[1], 7.0 / 11.0, epsilon = 1e-9);
67 }
68
69 #[test]
70 fn solves_one_by_one_system() {
71 let a = dense(1, vec![9.0]);
72 let out = solve_spd(&a, &[18.0]).unwrap();
73 assert_relative_eq!(out.solution[0], 2.0, epsilon = 1e-9);
74 }
75
76 #[test]
77 fn solves_via_constructed_rhs() {
78 // Pick A SPD and a known x; form b = A x; recover x.
79 let a = dense(3, vec![6.0, 2.0, 1.0, 2.0, 5.0, 2.0, 1.0, 2.0, 4.0]);
80 let x_true = [1.0, -2.0, 0.5];
81 let b = a.matvec(&x_true).unwrap();
82 let out = solve_spd(&a, &b).unwrap();
83 for (xi, ti) in out.solution.iter().zip(x_true) {
84 assert_relative_eq!(*xi, ti, epsilon = 1e-9);
85 }
86 }
87
88 #[test]
89 fn one_solver_serves_multiple_rhs() {
90 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
91 let cg = ConjugateGradient::new(&a, Config::new());
92 // First column of A as RHS recovers e1; second column recovers e2.
93 let x1 = cg.solve(&[4.0, 1.0]).unwrap().solution;
94 let x2 = cg.solve(&[1.0, 3.0]).unwrap().solution;
95 assert_relative_eq!(x1[0], 1.0, epsilon = 1e-9);
96 assert_relative_eq!(x1[1], 0.0, epsilon = 1e-9);
97 assert_relative_eq!(x2[0], 0.0, epsilon = 1e-9);
98 assert_relative_eq!(x2[1], 1.0, epsilon = 1e-9);
99 }
100
101 // ---------------------------------------------------------------------------
102 // Larger sparse system: the discrete Laplacian
103 // ---------------------------------------------------------------------------
104
105 #[test]
106 fn solves_large_laplacian_via_constructed_rhs() {
107 // A 50x50 tridiagonal Laplacian with a known solution. Sparse, SPD, and
108 // large enough that finite-termination and sparsity matter.
109 let n = 50;
110 let a = laplacian(n);
111 assert_eq!(a.nnz(), 3 * n - 2); // tridiagonal nnz count
112 let x_true: Vec<f64> = (0..n).map(|i| ((i as f64) * 0.1).sin()).collect();
113 let b = a.matvec(&x_true).unwrap();
114 let out = solve_spd(&a, &b).unwrap();
115 for (xi, ti) in out.solution.iter().zip(&x_true) {
116 assert!((xi - ti).abs() < 1e-7, "x={xi} expected {ti}");
117 }
118 // Residual must be tiny.
119 assert!(residual_norm(&a, &out.solution, &b) < 1e-8);
120 }
121
122 // ---------------------------------------------------------------------------
123 // Structural invariant 1: residual ‖A x - b‖ is tiny
124 // ---------------------------------------------------------------------------
125
126 #[test]
127 fn residual_is_tiny_for_spd_system() {
128 let a = dense(3, vec![25.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, 11.0]);
129 let b = [3.0, -7.0, 2.0];
130 let out = solve_spd(&a, &b).unwrap();
131 let r = residual_norm(&a, &out.solution, &b);
132 assert!(r < 1e-9, "residual norm {r} should be ~0");
133 // The reported residual norm agrees with the recomputed one.
134 assert_relative_eq!(out.residual_norm, r, epsilon = 1e-9);
135 }
136
137 #[test]
138 fn residual_tiny_across_many_rhs() {
139 let a = dense(
140 4,
141 vec![
142 10.0, 2.0, 3.0, 1.0, //
143 2.0, 9.0, 1.0, 2.0, //
144 3.0, 1.0, 12.0, 4.0, //
145 1.0, 2.0, 4.0, 8.0,
146 ],
147 );
148 let cg = ConjugateGradient::new(&a, Config::new());
149 for k in 0..6 {
150 let b: Vec<f64> = (0..4).map(|i| ((i * 7 + k * 3) as f64).sin()).collect();
151 let out = cg.solve(&b).unwrap();
152 let r = residual_norm(&a, &out.solution, &b);
153 assert!(r < 1e-9, "k={k} residual norm {r} should be ~0");
154 }
155 }
156
157 // ---------------------------------------------------------------------------
158 // Structural invariant 2: finite-termination convergence guarantee
159 // ---------------------------------------------------------------------------
160
161 #[test]
162 fn converges_within_n_iterations() {
163 // In exact arithmetic CG converges in at most n steps. Allowing a generous
164 // cap, a correct implementation reaches tolerance in <= n iterations for a
165 // well-conditioned small system.
166 let a = dense(3, vec![6.0, 2.0, 1.0, 2.0, 5.0, 2.0, 1.0, 2.0, 4.0]);
167 let out = solve_spd(&a, &[1.0, 2.0, 3.0]).unwrap();
168 assert!(
169 out.iterations <= 3,
170 "CG should converge within n=3 iterations, took {}",
171 out.iterations
172 );
173 assert!(out.converged);
174 }
175
176 #[test]
177 fn distinct_eigenvalues_converge_in_few_steps() {
178 // A diagonal matrix with two distinct eigenvalues: CG converges in at most
179 // 2 iterations (the number of distinct eigenvalues).
180 let a = dense(4, vec![
181 3.0, 0.0, 0.0, 0.0, //
182 0.0, 3.0, 0.0, 0.0, //
183 0.0, 0.0, 7.0, 0.0, //
184 0.0, 0.0, 0.0, 7.0,
185 ]);
186 let out = solve_spd(&a, &[1.0, 1.0, 1.0, 1.0]).unwrap();
187 assert!(
188 out.iterations <= 2,
189 "two distinct eigenvalues => <= 2 iterations, took {}",
190 out.iterations
191 );
192 }
193
194 #[test]
195 fn already_solved_rhs_terminates_immediately() {
196 // If the zero start is already the solution (b = 0 => x = 0), the iteration
197 // stops at 0 steps with zero residual.
198 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
199 let out = solve_spd(&a, &[0.0, 0.0]).unwrap();
200 assert_eq!(out.iterations, 0);
201 assert_relative_eq!(out.residual_norm, 0.0, epsilon = 1e-15);
202 assert_eq!(out.solution, vec![0.0, 0.0]);
203 }
204
205 // ---------------------------------------------------------------------------
206 // Structural invariant 3: linearity in the right-hand side
207 // ---------------------------------------------------------------------------
208
209 #[test]
210 fn solution_is_linear_in_rhs() {
211 // x(b1 + b2) = x(b1) + x(b2): a property any genuine linear solver has and
212 // a constant implementation violates.
213 let a = dense(3, vec![6.0, 2.0, 1.0, 2.0, 5.0, 2.0, 1.0, 2.0, 4.0]);
214 let cg = ConjugateGradient::new(&a, Config::new());
215 let b1 = [1.0, 0.0, 2.0];
216 let b2 = [-3.0, 4.0, 1.0];
217 let x1 = cg.solve(&b1).unwrap().solution;
218 let x2 = cg.solve(&b2).unwrap().solution;
219 let bsum: Vec<f64> = b1.iter().zip(b2).map(|(p, q)| p + q).collect();
220 let xsum = cg.solve(&bsum).unwrap().solution;
221 for i in 0..3 {
222 assert_relative_eq!(xsum[i], x1[i] + x2[i], epsilon = 1e-9);
223 }
224 }
225
226 #[test]
227 fn scaling_rhs_scales_solution() {
228 // x(c b) = c x(b).
229 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
230 let cg = ConjugateGradient::new(&a, Config::new());
231 let b = [1.0, 2.0];
232 let x = cg.solve(&b).unwrap().solution;
233 let scaled: Vec<f64> = b.iter().map(|v| 2.5 * v).collect();
234 let xs = cg.solve(&scaled).unwrap().solution;
235 for i in 0..2 {
236 assert_relative_eq!(xs[i], 2.5 * x[i], epsilon = 1e-9);
237 }
238 }
239
240 // ---------------------------------------------------------------------------
241 // Non-SPD detection (the core failure mode)
242 // ---------------------------------------------------------------------------
243
244 #[test]
245 fn rejects_indefinite_matrix() {
246 // Symmetric but indefinite (eigenvalues ±1, eigenvectors [1,1] and [1,-1]).
247 // A right-hand side along the negative-eigenvalue direction [1,-1] drives
248 // the search direction onto negative curvature pᵀAp < 0, so the iteration
249 // must report a breakdown.
250 let a = dense(2, vec![0.0, 1.0, 1.0, 0.0]);
251 let cg = ConjugateGradient::new(&a, Config::new());
252 let err = cg.solve(&[1.0, -1.0]).unwrap_err();
253 assert!(matches!(err, CgError::NotPositiveDefinite { .. }));
254 }
255
256 #[test]
257 fn rejects_negative_definite_matrix() {
258 // -I is negative-definite; the first curvature pᵀAp = -‖p‖² < 0.
259 let a = dense(2, vec![-1.0, 0.0, 0.0, -1.0]);
260 let cg = ConjugateGradient::new(&a, Config::new());
261 let err = cg.solve(&[1.0, 2.0]).unwrap_err();
262 match err {
263 CgError::NotPositiveDefinite { iteration, curvature } => {
264 assert_eq!(iteration, 0);
265 assert!(curvature <= 0.0, "curvature {curvature} should be non-positive");
266 }
267 other => panic!("expected NotPositiveDefinite, got {other:?}"),
268 }
269 }
270
271 // ---------------------------------------------------------------------------
272 // Non-convergence within the iteration budget
273 // ---------------------------------------------------------------------------
274
275 #[test]
276 fn reports_non_convergence_with_tight_iteration_cap() {
277 // A demanding tolerance with a one-iteration cap on a system that needs
278 // more steps cannot converge; the solve reports NotConverged.
279 let a = laplacian(20);
280 let x_true: Vec<f64> = (0..20).map(|i| (i as f64) + 1.0).collect();
281 let b = a.matvec(&x_true).unwrap();
282 let cfg = Config::new().with_tolerance(1e-14).with_max_iterations(1);
283 let cg = ConjugateGradient::new(&a, cfg);
284 let err = cg.solve(&b).unwrap_err();
285 match err {
286 CgError::NotConverged { max_iterations, residual_norm } => {
287 assert_eq!(max_iterations, 1);
288 assert!(residual_norm.is_finite());
289 }
290 other => panic!("expected NotConverged, got {other:?}"),
291 }
292 }
293
294 // ---------------------------------------------------------------------------
295 // Config / input validation
296 // ---------------------------------------------------------------------------
297
298 #[test]
299 fn rejects_non_positive_tolerance() {
300 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
301 let cfg = Config::new().with_tolerance(0.0);
302 let err = ConjugateGradient::new(&a, cfg).solve(&[1.0, 2.0]).unwrap_err();
303 assert!(matches!(err, CgError::InvalidTolerance(_)));
304 }
305
306 #[test]
307 fn rejects_rhs_dimension_mismatch() {
308 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
309 let cg = ConjugateGradient::new(&a, Config::new());
310 let err = cg.solve(&[1.0, 2.0, 3.0]).unwrap_err();
311 assert!(matches!(
312 err,
313 CgError::DimensionMismatch {
314 expected: 2,
315 got: 3
316 }
317 ));
318 }
319
320 #[test]
321 fn rejects_empty_matrix_at_construction() {
322 let err = SparseMatrix::from_triplets(0, &[]).unwrap_err();
323 assert_eq!(err, CgError::EmptyMatrix);
324 }
325
326 #[test]
327 fn rejects_dense_shape_mismatch_at_construction() {
328 let err = SparseMatrix::from_dense(2, &[1.0, 2.0, 3.0]).unwrap_err();
329 assert!(matches!(err, CgError::DataShapeMismatch { len: 3, .. }));
330 }
331
332 #[test]
333 fn rejects_out_of_bounds_triplet() {
334 let err = SparseMatrix::from_triplets(2, &[(0, 5, 1.0)]).unwrap_err();
335 assert!(matches!(
336 err,
337 CgError::IndexOutOfBounds { row: 0, col: 5, dim: 2 }
338 ));
339 }
340
341 #[test]
342 fn rejects_non_finite_entry_at_construction() {
343 let err = SparseMatrix::from_triplets(2, &[(1, 1, f64::NAN)]).unwrap_err();
344 assert!(matches!(
345 err,
346 CgError::NonFiniteEntry { row: 1, col: 1, .. }
347 ));
348 }
349
350 // ---------------------------------------------------------------------------
351 // Sparse matrix behaviour
352 // ---------------------------------------------------------------------------
353
354 #[test]
355 fn duplicate_triplets_are_summed() {
356 // Two contributions to (0,0) sum to 5; off-diagonal stays separate.
357 let a = SparseMatrix::from_triplets(2, &[(0, 0, 2.0), (0, 0, 3.0), (1, 1, 1.0)]).unwrap();
358 assert_relative_eq!(a.get(0, 0), 5.0, epsilon = 1e-15);
359 assert_relative_eq!(a.get(1, 1), 1.0, epsilon = 1e-15);
360 }
361
362 #[test]
363 fn symmetry_check_distinguishes() {
364 let sym = dense(2, vec![2.0, 1.0, 1.0, 3.0]);
365 assert!(sym.is_symmetric(1e-12));
366 let asym = dense(2, vec![2.0, 1.0, 9.0, 3.0]);
367 assert!(!asym.is_symmetric(1e-12));
368 }
369
370 #[test]
371 fn matvec_matches_dense_definition() {
372 let a = dense(3, vec![1.0, 2.0, 0.0, 0.0, 3.0, 4.0, 5.0, 0.0, 6.0]);
373 let y = a.matvec(&[1.0, 1.0, 1.0]).unwrap();
374 assert_eq!(y, vec![3.0, 7.0, 11.0]);
375 }
376
377 // ---------------------------------------------------------------------------
378 // CSR assembly + product: out-of-order columns and duplicate coordinates.
379 // ---------------------------------------------------------------------------
380
381 #[test]
382 fn from_triplets_sums_dups_orders_columns_and_multiplies() {
383 // A = [[3,0,1],[0,4,0],[1,0,5]] given with a duplicated (0,0) and columns
384 // out of order. from_triplets must sum the duplicate and order columns; the
385 // accessors and matvec must then be correct.
386 let a = SparseMatrix::from_triplets(3, &[
387 (0, 2, 1.0),
388 (2, 2, 5.0),
389 (0, 0, 2.0),
390 (1, 1, 4.0),
391 (0, 0, 1.0),
392 (2, 0, 1.0),
393 ])
394 .unwrap();
395 assert_eq!(a.get(0, 0), 3.0);
396 assert_eq!(a.get(0, 2), 1.0);
397 assert_eq!(a.get(2, 0), 1.0);
398 assert!(a.is_symmetric(1e-12));
399 let y = a.matvec(&[1.0, 1.0, 1.0]).unwrap();
400 assert_relative_eq!(y[0], 4.0, epsilon = 1e-12);
401 assert_relative_eq!(y[1], 4.0, epsilon = 1e-12);
402 assert_relative_eq!(y[2], 6.0, epsilon = 1e-12);
403 // And the solver works on the assembled operator.
404 let out = solve_spd(&a, &[4.0, 4.0, 6.0]).unwrap();
405 assert_relative_eq!(out.solution[0], 1.0, epsilon = 1e-9);
406 assert_relative_eq!(out.solution[1], 1.0, epsilon = 1e-9);
407 assert_relative_eq!(out.solution[2], 1.0, epsilon = 1e-9);
408 }
409
410 #[test]
411 fn three_distinct_eigenvalues_converge_in_at_most_three() {
412 // diag with eigenvalues {2,5,9} (each repeated): a Krylov method terminates
413 // in at most #distinct = 3 iterations. A non-Krylov iterative scheme
414 // (steepest descent, Jacobi) would need many more, so this pins CG behaviour.
415 let a = dense(5, vec![
416 2.0, 0.0, 0.0, 0.0, 0.0,
417 0.0, 2.0, 0.0, 0.0, 0.0,
418 0.0, 0.0, 5.0, 0.0, 0.0,
419 0.0, 0.0, 0.0, 9.0, 0.0,
420 0.0, 0.0, 0.0, 0.0, 9.0,
421 ]);
422 let out = solve_spd(&a, &[1.0, 1.0, 1.0, 1.0, 1.0]).unwrap();
423 assert!(out.iterations <= 3, "three distinct eigenvalues => <= 3 iterations, took {}", out.iterations);
424 assert_relative_eq!(out.solution[2], 0.2, epsilon = 1e-9);
425 assert_relative_eq!(out.solution[3], 1.0 / 9.0, epsilon = 1e-9);
426 }
427
1 //! Configuration for a conjugate-gradient solve.
2
3 use crate::error::CgError;
4
5 /// Default relative residual tolerance used by [`Config::new`].
6 ///
7 /// The iteration stops when `‖b - A x‖ <= tolerance · ‖b‖`.
8 pub const DEFAULT_TOLERANCE: f64 = 1e-10;
9
10 /// Default cap, as a multiple of the system dimension `n`, on the number of
11 /// iterations used by [`Config::new`]. In exact arithmetic conjugate gradients
12 /// converge in at most `n` steps; the extra slack absorbs rounding.
13 pub const DEFAULT_MAX_ITER_FACTOR: usize = 2;
14
15 /// Tuning parameters for a conjugate-gradient solve.
16 ///
17 /// A `Config` bundles the stopping criteria. Construct one with [`Config::new`]
18 /// (sensible defaults derived from the system dimension) and refine it with the
19 /// chained setters, e.g.
20 ///
21 /// ```
22 /// use cgsolve::Config;
23 /// let cfg = Config::new()
24 /// .with_tolerance(1e-8)
25 /// .with_max_iterations(100);
26 /// assert_eq!(cfg.tolerance(), 1e-8);
27 /// assert_eq!(cfg.max_iterations(), Some(100));
28 /// ```
29 #[derive(Debug, Clone, Copy, PartialEq)]
30 pub struct Config {
31 tolerance: f64,
32 max_iterations: Option<usize>,
33 }
34
35 impl Config {
36 /// Create a configuration with the crate default tolerance
37 /// ([`DEFAULT_TOLERANCE`]) and an automatic iteration cap (derived from the
38 /// system dimension when the solve runs).
39 pub fn new() -> Self {
40 Self {
41 tolerance: DEFAULT_TOLERANCE,
42 max_iterations: None,
43 }
44 }
45
46 /// Set the relative residual tolerance.
47 ///
48 /// Smaller values demand a more accurate result at the cost of more
49 /// iterations.
50 #[must_use]
51 pub fn with_tolerance(mut self, tolerance: f64) -> Self {
52 self.tolerance = tolerance;
53 self
54 }
55
56 /// Set an explicit maximum iteration count, overriding the automatic cap.
57 #[must_use]
58 pub fn with_max_iterations(mut self, max_iterations: usize) -> Self {
59 self.max_iterations = Some(max_iterations);
60 self
61 }
62
63 /// The configured relative residual tolerance.
64 pub fn tolerance(&self) -> f64 {
65 self.tolerance
66 }
67
68 /// The explicit iteration cap, if one was set.
69 pub fn max_iterations(&self) -> Option<usize> {
70 self.max_iterations
71 }
72
73 /// Resolve the effective iteration cap for a system of dimension `n`:
74 /// the explicit cap if set, otherwise `DEFAULT_MAX_ITER_FACTOR · n + 1`.
75 pub fn effective_max_iterations(&self, n: usize) -> usize {
76 self.max_iterations
77 .unwrap_or(DEFAULT_MAX_ITER_FACTOR * n + 1)
78 }
79
80 /// Validate the configuration, returning an error if the tolerance is out
81 /// of range. Called internally before a solve begins.
82 pub(crate) fn validate(&self) -> Result<(), CgError> {
83 if !self.tolerance.is_finite() || self.tolerance <= 0.0 {
84 return Err(CgError::InvalidTolerance(self.tolerance));
85 }
86 Ok(())
87 }
88 }
89
90 impl Default for Config {
91 fn default() -> Self {
92 Self::new()
93 }
94 }
95
/workspace/cgsolve/src/config.rs
1 //! Integration tests for the `cgsolve` conjugate-gradient solver.
2 //!
3 //! These exercise the public API against systems with known solutions, a large
4 //! sparse system, edge cases, every error path, and three structural
5 //! properties , the residual `b - A x ≈ 0`, the finite-termination convergence
6 //! guarantee (CG converges within `n` iterations), and linearity in the
7 //! right-hand side. A correct CG implementation passes all of them; a hardcoded
8 //! or trivially constant implementation does not.
9
10 use approx::assert_relative_eq;
11 use cgsolve::{solve_spd, CgError, Config, ConjugateGradient, SparseMatrix};
12
13 /// Build a dense-specified SPD matrix, panicking on malformed input.
14 fn dense(n: usize, data: Vec<f64>) -> SparseMatrix {
15 SparseMatrix::from_dense(n, &data).expect("valid matrix literal")
16 }
17
18 /// Euclidean norm of a vector.
19 fn norm(v: &[f64]) -> f64 {
20 v.iter().map(|x| x * x).sum::<f64>().sqrt()
21 }
22
23 /// Residual norm ‖A x - b‖ for a candidate solution `x`.
24 fn residual_norm(a: &SparseMatrix, x: &[f64], b: &[f64]) -> f64 {
25 let ax = a.matvec(x).expect("dimensions match");
26 let r: Vec<f64> = ax.iter().zip(b).map(|(p, q)| p - q).collect();
27 norm(&r)
28 }
29
30 /// Build the 1-D discrete Laplacian `tridiag(-1, 2, -1)` of dimension `n`.
31 /// This is a classic large sparse SPD operator (the matrix of the second
32 /// difference), positive-definite with eigenvalues `2 - 2 cos(kπ/(n+1))`.
33 fn laplacian(n: usize) -> SparseMatrix {
34 let mut triplets = Vec::new();
35 for i in 0..n {
36 triplets.push((i, i, 2.0));
37 if i + 1 < n {
38 triplets.push((i, i + 1, -1.0));
39 triplets.push((i + 1, i, -1.0));
40 }
41 }
42 SparseMatrix::from_triplets(n, &triplets).expect("valid Laplacian")
43 }
44
45 // ---------------------------------------------------------------------------
46 // Known systems with hand-checked solutions
47 // ---------------------------------------------------------------------------
48
49 #[test]
50 fn solves_diagonal_system() {
51 // diag(2,4,5) x = (6,8,10) => x = (3,2,2).
52 let a = dense(3, vec![2.0, 0.0, 0.0, 0.0, 4.0, 0.0, 0.0, 0.0, 5.0]);
53 let out = solve_spd(&a, &[6.0, 8.0, 10.0]).unwrap();
54 assert_relative_eq!(out.solution[0], 3.0, epsilon = 1e-9);
55 assert_relative_eq!(out.solution[1], 2.0, epsilon = 1e-9);
56 assert_relative_eq!(out.solution[2], 2.0, epsilon = 1e-9);
57 }
58
59 #[test]
60 fn solves_known_2x2_system() {
61 // A = [[4, 1], [1, 3]], b = [1, 2]. Closed form: det = 11,
62 // x0 = (1*3 - 1*2)/11 = 1/11, x1 = (4*2 - 1*1)/11 = 7/11.
63 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
64 let out = solve_spd(&a, &[1.0, 2.0]).unwrap();
65 assert_relative_eq!(out.solution[0], 1.0 / 11.0, epsilon = 1e-9);
66 assert_relative_eq!(out.solution[1], 7.0 / 11.0, epsilon = 1e-9);
67 }
68
69 #[test]
70 fn solves_one_by_one_system() {
71 let a = dense(1, vec![9.0]);
72 let out = solve_spd(&a, &[18.0]).unwrap();
73 assert_relative_eq!(out.solution[0], 2.0, epsilon = 1e-9);
74 }
75
76 #[test]
77 fn solves_via_constructed_rhs() {
78 // Pick A SPD and a known x; form b = A x; recover x.
79 let a = dense(3, vec![6.0, 2.0, 1.0, 2.0, 5.0, 2.0, 1.0, 2.0, 4.0]);
80 let x_true = [1.0, -2.0, 0.5];
81 let b = a.matvec(&x_true).unwrap();
82 let out = solve_spd(&a, &b).unwrap();
83 for (xi, ti) in out.solution.iter().zip(x_true) {
84 assert_relative_eq!(*xi, ti, epsilon = 1e-9);
85 }
86 }
87
88 #[test]
89 fn one_solver_serves_multiple_rhs() {
90 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
91 let cg = ConjugateGradient::new(&a, Config::new());
92 // First column of A as RHS recovers e1; second column recovers e2.
93 let x1 = cg.solve(&[4.0, 1.0]).unwrap().solution;
94 let x2 = cg.solve(&[1.0, 3.0]).unwrap().solution;
95 assert_relative_eq!(x1[0], 1.0, epsilon = 1e-9);
96 assert_relative_eq!(x1[1], 0.0, epsilon = 1e-9);
97 assert_relative_eq!(x2[0], 0.0, epsilon = 1e-9);
98 assert_relative_eq!(x2[1], 1.0, epsilon = 1e-9);
99 }
100
101 // ---------------------------------------------------------------------------
102 // Larger sparse system: the discrete Laplacian
103 // ---------------------------------------------------------------------------
104
105 #[test]
106 fn solves_large_laplacian_via_constructed_rhs() {
107 // A 50x50 tridiagonal Laplacian with a known solution. Sparse, SPD, and
108 // large enough that finite-termination and sparsity matter.
109 let n = 50;
110 let a = laplacian(n);
111 assert_eq!(a.nnz(), 3 * n - 2); // tridiagonal nnz count
112 let x_true: Vec<f64> = (0..n).map(|i| ((i as f64) * 0.1).sin()).collect();
113 let b = a.matvec(&x_true).unwrap();
114 let out = solve_spd(&a, &b).unwrap();
115 for (xi, ti) in out.solution.iter().zip(&x_true) {
116 assert!((xi - ti).abs() < 1e-7, "x={xi} expected {ti}");
117 }
118 // Residual must be tiny.
119 assert!(residual_norm(&a, &out.solution, &b) < 1e-8);
120 }
121
122 // ---------------------------------------------------------------------------
123 // Structural invariant 1: residual ‖A x - b‖ is tiny
124 // ---------------------------------------------------------------------------
125
126 #[test]
127 fn residual_is_tiny_for_spd_system() {
128 let a = dense(3, vec![25.0, 15.0, -5.0, 15.0, 18.0, 0.0, -5.0, 0.0, 11.0]);
129 let b = [3.0, -7.0, 2.0];
130 let out = solve_spd(&a, &b).unwrap();
131 let r = residual_norm(&a, &out.solution, &b);
132 assert!(r < 1e-9, "residual norm {r} should be ~0");
133 // The reported residual norm agrees with the recomputed one.
134 assert_relative_eq!(out.residual_norm, r, epsilon = 1e-9);
135 }
136
137 #[test]
138 fn residual_tiny_across_many_rhs() {
139 let a = dense(
140 4,
141 vec![
142 10.0, 2.0, 3.0, 1.0, //
143 2.0, 9.0, 1.0, 2.0, //
144 3.0, 1.0, 12.0, 4.0, //
145 1.0, 2.0, 4.0, 8.0,
146 ],
147 );
148 let cg = ConjugateGradient::new(&a, Config::new());
149 for k in 0..6 {
150 let b: Vec<f64> = (0..4).map(|i| ((i * 7 + k * 3) as f64).sin()).collect();
151 let out = cg.solve(&b).unwrap();
152 let r = residual_norm(&a, &out.solution, &b);
153 assert!(r < 1e-9, "k={k} residual norm {r} should be ~0");
154 }
155 }
156
157 // ---------------------------------------------------------------------------
158 // Structural invariant 2: finite-termination convergence guarantee
159 // ---------------------------------------------------------------------------
160
161 #[test]
162 fn converges_within_n_iterations() {
163 // In exact arithmetic CG converges in at most n steps. Allowing a generous
164 // cap, a correct implementation reaches tolerance in <= n iterations for a
165 // well-conditioned small system.
166 let a = dense(3, vec![6.0, 2.0, 1.0, 2.0, 5.0, 2.0, 1.0, 2.0, 4.0]);
167 let out = solve_spd(&a, &[1.0, 2.0, 3.0]).unwrap();
168 assert!(
169 out.iterations <= 3,
170 "CG should converge within n=3 iterations, took {}",
171 out.iterations
172 );
173 assert!(out.converged);
174 }
175
176 #[test]
177 fn distinct_eigenvalues_converge_in_few_steps() {
178 // A diagonal matrix with two distinct eigenvalues: CG converges in at most
179 // 2 iterations (the number of distinct eigenvalues).
180 let a = dense(4, vec![
181 3.0, 0.0, 0.0, 0.0, //
182 0.0, 3.0, 0.0, 0.0, //
183 0.0, 0.0, 7.0, 0.0, //
184 0.0, 0.0, 0.0, 7.0,
185 ]);
186 let out = solve_spd(&a, &[1.0, 1.0, 1.0, 1.0]).unwrap();
187 assert!(
188 out.iterations <= 2,
189 "two distinct eigenvalues => <= 2 iterations, took {}",
190 out.iterations
191 );
192 }
193
194 #[test]
195 fn already_solved_rhs_terminates_immediately() {
196 // If the zero start is already the solution (b = 0 => x = 0), the iteration
197 // stops at 0 steps with zero residual.
198 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
199 let out = solve_spd(&a, &[0.0, 0.0]).unwrap();
200 assert_eq!(out.iterations, 0);
201 assert_relative_eq!(out.residual_norm, 0.0, epsilon = 1e-15);
202 assert_eq!(out.solution, vec![0.0, 0.0]);
203 }
204
205 // ---------------------------------------------------------------------------
206 // Structural invariant 3: linearity in the right-hand side
207 // ---------------------------------------------------------------------------
208
209 #[test]
210 fn solution_is_linear_in_rhs() {
211 // x(b1 + b2) = x(b1) + x(b2): a property any genuine linear solver has and
212 // a constant implementation violates.
213 let a = dense(3, vec![6.0, 2.0, 1.0, 2.0, 5.0, 2.0, 1.0, 2.0, 4.0]);
214 let cg = ConjugateGradient::new(&a, Config::new());
215 let b1 = [1.0, 0.0, 2.0];
216 let b2 = [-3.0, 4.0, 1.0];
217 let x1 = cg.solve(&b1).unwrap().solution;
218 let x2 = cg.solve(&b2).unwrap().solution;
219 let bsum: Vec<f64> = b1.iter().zip(b2).map(|(p, q)| p + q).collect();
220 let xsum = cg.solve(&bsum).unwrap().solution;
221 for i in 0..3 {
222 assert_relative_eq!(xsum[i], x1[i] + x2[i], epsilon = 1e-9);
223 }
224 }
225
226 #[test]
227 fn scaling_rhs_scales_solution() {
228 // x(c b) = c x(b).
229 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
230 let cg = ConjugateGradient::new(&a, Config::new());
231 let b = [1.0, 2.0];
232 let x = cg.solve(&b).unwrap().solution;
233 let scaled: Vec<f64> = b.iter().map(|v| 2.5 * v).collect();
234 let xs = cg.solve(&scaled).unwrap().solution;
235 for i in 0..2 {
236 assert_relative_eq!(xs[i], 2.5 * x[i], epsilon = 1e-9);
237 }
238 }
239
240 // ---------------------------------------------------------------------------
241 // Non-SPD detection (the core failure mode)
242 // ---------------------------------------------------------------------------
243
244 #[test]
245 fn rejects_indefinite_matrix() {
246 // Symmetric but indefinite (eigenvalues ±1, eigenvectors [1,1] and [1,-1]).
247 // A right-hand side along the negative-eigenvalue direction [1,-1] drives
248 // the search direction onto negative curvature pᵀAp < 0, so the iteration
249 // must report a breakdown.
250 let a = dense(2, vec![0.0, 1.0, 1.0, 0.0]);
251 let cg = ConjugateGradient::new(&a, Config::new());
252 let err = cg.solve(&[1.0, -1.0]).unwrap_err();
253 assert!(matches!(err, CgError::NotPositiveDefinite { .. }));
254 }
255
256 #[test]
257 fn rejects_negative_definite_matrix() {
258 // -I is negative-definite; the first curvature pᵀAp = -‖p‖² < 0.
259 let a = dense(2, vec![-1.0, 0.0, 0.0, -1.0]);
260 let cg = ConjugateGradient::new(&a, Config::new());
261 let err = cg.solve(&[1.0, 2.0]).unwrap_err();
262 match err {
263 CgError::NotPositiveDefinite { iteration, curvature } => {
264 assert_eq!(iteration, 0);
265 assert!(curvature <= 0.0, "curvature {curvature} should be non-positive");
266 }
267 other => panic!("expected NotPositiveDefinite, got {other:?}"),
268 }
269 }
270
271 // ---------------------------------------------------------------------------
272 // Non-convergence within the iteration budget
273 // ---------------------------------------------------------------------------
274
275 #[test]
276 fn reports_non_convergence_with_tight_iteration_cap() {
277 // A demanding tolerance with a one-iteration cap on a system that needs
278 // more steps cannot converge; the solve reports NotConverged.
279 let a = laplacian(20);
280 let x_true: Vec<f64> = (0..20).map(|i| (i as f64) + 1.0).collect();
281 let b = a.matvec(&x_true).unwrap();
282 let cfg = Config::new().with_tolerance(1e-14).with_max_iterations(1);
283 let cg = ConjugateGradient::new(&a, cfg);
284 let err = cg.solve(&b).unwrap_err();
285 match err {
286 CgError::NotConverged { max_iterations, residual_norm } => {
287 assert_eq!(max_iterations, 1);
288 assert!(residual_norm.is_finite());
289 }
290 other => panic!("expected NotConverged, got {other:?}"),
291 }
292 }
293
294 // ---------------------------------------------------------------------------
295 // Config / input validation
296 // ---------------------------------------------------------------------------
297
298 #[test]
299 fn rejects_non_positive_tolerance() {
300 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
301 let cfg = Config::new().with_tolerance(0.0);
302 let err = ConjugateGradient::new(&a, cfg).solve(&[1.0, 2.0]).unwrap_err();
303 assert!(matches!(err, CgError::InvalidTolerance(_)));
304 }
305
306 #[test]
307 fn rejects_rhs_dimension_mismatch() {
308 let a = dense(2, vec![4.0, 1.0, 1.0, 3.0]);
309 let cg = ConjugateGradient::new(&a, Config::new());
310 let err = cg.solve(&[1.0, 2.0, 3.0]).unwrap_err();
311 assert!(matches!(
312 err,
313 CgError::DimensionMismatch {
314 expected: 2,
315 got: 3
316 }
317 ));
318 }
319
320 #[test]
321 fn rejects_empty_matrix_at_construction() {
322 let err = SparseMatrix::from_triplets(0, &[]).unwrap_err();
323 assert_eq!(err, CgError::EmptyMatrix);
324 }
325
326 #[test]
327 fn rejects_dense_shape_mismatch_at_construction() {
328 let err = SparseMatrix::from_dense(2, &[1.0, 2.0, 3.0]).unwrap_err();
329 assert!(matches!(err, CgError::DataShapeMismatch { len: 3, .. }));
330 }
331
332 #[test]
333 fn rejects_out_of_bounds_triplet() {
334 let err = SparseMatrix::from_triplets(2, &[(0, 5, 1.0)]).unwrap_err();
335 assert!(matches!(
336 err,
337 CgError::IndexOutOfBounds { row: 0, col: 5, dim: 2 }
338 ));
339 }
340
341 #[test]
342 fn rejects_non_finite_entry_at_construction() {
343 let err = SparseMatrix::from_triplets(2, &[(1, 1, f64::NAN)]).unwrap_err();
344 assert!(matches!(
345 err,
346 CgError::NonFiniteEntry { row: 1, col: 1, .. }
347 ));
348 }
349
350 // ---------------------------------------------------------------------------
351 // Sparse matrix behaviour
352 // ---------------------------------------------------------------------------
353
354 #[test]
355 fn duplicate_triplets_are_summed() {
356 // Two contributions to (0,0) sum to 5; off-diagonal stays separate.
357 let a = SparseMatrix::from_triplets(2, &[(0, 0, 2.0), (0, 0, 3.0), (1, 1, 1.0)]).unwrap();
358 assert_relative_eq!(a.get(0, 0), 5.0, epsilon = 1e-15);
359 assert_relative_eq!(a.get(1, 1), 1.0, epsilon = 1e-15);
360 }
361
362 #[test]
363 fn symmetry_check_distinguishes() {
364 let sym = dense(2, vec![2.0, 1.0, 1.0, 3.0]);
365 assert!(sym.is_symmetric(1e-12));
366 let asym = dense(2, vec![2.0, 1.0, 9.0, 3.0]);
367 assert!(!asym.is_symmetric(1e-12));
368 }
369
370 #[test]
371 fn matvec_matches_dense_definition() {
372 let a = dense(3, vec![1.0, 2.0, 0.0, 0.0, 3.0, 4.0, 5.0, 0.0, 6.0]);
373 let y = a.matvec(&[1.0, 1.0, 1.0]).unwrap();
374 assert_eq!(y, vec![3.0, 7.0, 11.0]);
375 }
376
377 // ---------------------------------------------------------------------------
378 // CSR assembly + product: out-of-order columns and duplicate coordinates.
379 // ---------------------------------------------------------------------------
380
381 #[test]
382 fn from_triplets_sums_dups_orders_columns_and_multiplies() {
383 // A = [[3,0,1],[0,4,0],[1,0,5]] given with a duplicated (0,0) and columns
384 // out of order. from_triplets must sum the duplicate and order columns; the
385 // accessors and matvec must then be correct.
386 let a = SparseMatrix::from_triplets(3, &[
387 (0, 2, 1.0),
388 (2, 2, 5.0),
389 (0, 0, 2.0),
390 (1, 1, 4.0),
391 (0, 0, 1.0),
392 (2, 0, 1.0),
393 ])
394 .unwrap();
395 assert_eq!(a.get(0, 0), 3.0);
396 assert_eq!(a.get(0, 2), 1.0);
397 assert_eq!(a.get(2, 0), 1.0);
398 assert!(a.is_symmetric(1e-12));
399 let y = a.matvec(&[1.0, 1.0, 1.0]).unwrap();
400 assert_relative_eq!(y[0], 4.0, epsilon = 1e-12);
401 assert_relative_eq!(y[1], 4.0, epsilon = 1e-12);
402 assert_relative_eq!(y[2], 6.0, epsilon = 1e-12);
403 // And the solver works on the assembled operator.
404 let out = solve_spd(&a, &[4.0, 4.0, 6.0]).unwrap();
405 assert_relative_eq!(out.solution[0], 1.0, epsilon = 1e-9);
406 assert_relative_eq!(out.solution[1], 1.0, epsilon = 1e-9);
407 assert_relative_eq!(out.solution[2], 1.0, epsilon = 1e-9);
408 }
409
410 #[test]
411 fn three_distinct_eigenvalues_converge_in_at_most_three() {
412 // diag with eigenvalues {2,5,9} (each repeated): a Krylov method terminates
413 // in at most #distinct = 3 iterations. A non-Krylov iterative scheme
414 // (steepest descent, Jacobi) would need many more, so this pins CG behaviour.
415 let a = dense(5, vec![
416 2.0, 0.0, 0.0, 0.0, 0.0,
417 0.0, 2.0, 0.0, 0.0, 0.0,
418 0.0, 0.0, 5.0, 0.0, 0.0,
419 0.0, 0.0, 0.0, 9.0, 0.0,
420 0.0, 0.0, 0.0, 0.0, 9.0,
421 ]);
422 let out = solve_spd(&a, &[1.0, 1.0, 1.0, 1.0, 1.0]).unwrap();
423 assert!(out.iterations <= 3, "three distinct eigenvalues => <= 3 iterations, took {}", out.iterations);
424 assert_relative_eq!(out.solution[2], 0.2, epsilon = 1e-9);
425 assert_relative_eq!(out.solution[3], 1.0 / 9.0, epsilon = 1e-9);
426 }
427
1 //! Configuration for a conjugate-gradient solve.
2
3 use crate::error::CgError;
4
5 /// Default relative residual tolerance used by [`Config::new`].
6 ///
7 /// The iteration stops when `‖b - A x‖ <= tolerance · ‖b‖`.
8 pub const DEFAULT_TOLERANCE: f64 = 1e-10;
9
10 /// Default cap, as a multiple of the system dimension `n`, on the number of
11 /// iterations used by [`Config::new`]. In exact arithmetic conjugate gradients
12 /// converge in at most `n` steps; the extra slack absorbs rounding.
13 pub const DEFAULT_MAX_ITER_FACTOR: usize = 2;
14
15 /// Tuning parameters for a conjugate-gradient solve.
16 ///
17 /// A `Config` bundles the stopping criteria. Construct one with [`Config::new`]
18 /// (sensible defaults derived from the system dimension) and refine it with the
19 /// chained setters, e.g.
20 ///
21 /// ```
22 /// use cgsolve::Config;
23 /// let cfg = Config::new()
24 /// .with_tolerance(1e-8)
25 /// .with_max_iterations(100);
26 /// assert_eq!(cfg.tolerance(), 1e-8);
27 /// assert_eq!(cfg.max_iterations(), Some(100));
28 /// ```
29 #[derive(Debug, Clone, Copy, PartialEq)]
30 pub struct Config {
31 tolerance: f64,
32 max_iterations: Option<usize>,
33 }
34
35 impl Config {
36 /// Create a configuration with the crate default tolerance
37 /// ([`DEFAULT_TOLERANCE`]) and an automatic iteration cap (derived from the
38 /// system dimension when the solve runs).
39 pub fn new() -> Self {
40 Self {
41 tolerance: DEFAULT_TOLERANCE,
42 max_iterations: None,
43 }
44 }
45
46 /// Set the relative residual tolerance.
47 ///
48 /// Smaller values demand a more accurate result at the cost of more
49 /// iterations.
50 #[must_use]
51 pub fn with_tolerance(mut self, tolerance: f64) -> Self {
52 self.tolerance = tolerance;
53 self
54 }
55
56 /// Set an explicit maximum iteration count, overriding the automatic cap.
57 #[must_use]
58 pub fn with_max_iterations(mut self, max_iterations: usize) -> Self {
59 self.max_iterations = Some(max_iterations);
60 self
61 }
62
63 /// The configured relative residual tolerance.
64 pub fn tolerance(&self) -> f64 {
65 self.tolerance
66 }
67
68 /// The explicit iteration cap, if one was set.
69 pub fn max_iterations(&self) -> Option<usize> {
70 self.max_iterations
71 }
72
73 /// Resolve the effective iteration cap for a system of dimension `n`:
74 /// the explicit cap if set, otherwise `DEFAULT_MAX_ITER_FACTOR · n + 1`.
75 pub fn effective_max_iterations(&self, n: usize) -> usize {
76 self.max_iterations
77 .unwrap_or(DEFAULT_MAX_ITER_FACTOR * n + 1)
78 }
79
80 /// Validate the configuration, returning an error if the tolerance is out
81 /// of range. Called internally before a solve begins.
82 pub(crate) fn validate(&self) -> Result<(), CgError> {
83 if !self.tolerance.is_finite() || self.tolerance <= 0.0 {
84 return Err(CgError::InvalidTolerance(self.tolerance));
85 }
86 Ok(())
87 }
88 }
89
90 impl Default for Config {
91 fn default() -> Self {
92 Self::new()
93 }
94 }
95
/workspace/cgsolve/src/matrix.rs
pub fn from_triplets(n: usize, triplets: &[(usize, usize, f64)]) -> Result<Self, CgError> {
let _ = (n, triplets);
todo!("implement from_triplets (sci-4519)")
} pub fn from_triplets(n: usize, triplets: &[(usize, usize, f64)]) -> Result<Self, CgError> {
if n == 0 {
return Err(CgError::EmptyMatrix);
}
for &(row, col, value) in triplets {
if row >= n || col >= n {
return Err(CgError::IndexOutOfBounds { row, col, dim: n });
}
if !value.is_finite() {
return Err(CgError::NonFiniteEntry { row, col, value });
}
}
// Sort by (row, col) and merge duplicates by summing.
let mut sorted: Vec<(usize, usize, f64)> = triplets.to_vec();
sorted.sort_by_key(|&(r, c, _)| (r, c));
let mut merged: Vec<(usize, usize, f64)> = Vec::with_capacity(sorted.len());
for &(row, col, val) in &sorted {
match merged.last_mut() {
Some(last) if last.0 == row && last.1 == col => last.2 += val,
_ => merged.push((row, col, val)),
}
}
// Build CSR arrays.
let mut row_ptr = vec![0usize; n + 1];
for &(row, _, _) in &merged {
row_ptr[row + 1] += 1;
}
for i in 0..n {
row_ptr[i + 1] += row_ptr[i];
}
let col_idx: Vec<usize> = merged.iter().map(|&(_, c, _)| c).collect();
let values: Vec<f64> = merged.iter().map(|&(_, _, v)| v).collect();
Ok(Self { n, row_ptr, col_idx, values })
}/workspace/cgsolve/src/matrix.rs
pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
let _ = (x, out);
todo!("implement matvec_into (sci-4519)")
} pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
if x.len() != self.n {
return Err(CgError::DimensionMismatch { expected: self.n, got: x.len() });
}
if out.len() != self.n {
return Err(CgError::DimensionMismatch { expected: self.n, got: out.len() });
}
for i in 0..self.n {
let mut sum = 0.0;
for k in self.row_ptr[i]..self.row_ptr[i + 1] {
sum += self.values[k] * x[self.col_idx[k]];
}
out[i] = sum;
}
Ok(())
}/workspace/cgsolve/src/cg.rs
// TODO(sci-4519): implement the conjugate-gradient iteration described in
// the doc comment above. Form the initial residual r = b - A x, iterate the
// alpha/x/r/beta/p updates using one matvec per step, stop on the relative
// residual threshold, return NotPositiveDefinite on non-positive curvature
// and NotConverged if the budget is exhausted. See `tests/integration.rs`
// for the contract under test. The `dot` helper below computes the
// Euclidean dot products (rᵀr, pᵀAp, ...) you will need.
let _ = (a, b, x, tolerance, max_iterations, dot);
todo!("implement cg_iterate (sci-4519)") let n = b.len();
// r = b - A x (x is the zero vector on entry, so r = b)
let ax0 = a.matvec(x)?;
let mut r: Vec<f64> = b.iter().zip(&ax0).map(|(bi, ai)| bi - ai).collect();
// Stopping threshold: relative, falling back to absolute when b = 0.
let b_norm = dot(b, b).sqrt();
let scale = if b_norm > 0.0 { b_norm } else { 1.0 };
let threshold = tolerance * scale;
let mut rr = dot(&r, &r);
let init_residual = rr.sqrt();
if init_residual <= threshold {
return Ok((0, init_residual));
}
let mut p = r.clone();
let mut ap = vec![0.0; n];
for k in 0..max_iterations {
a.matvec_into(&p, &mut ap)?;
let curvature = dot(&p, &ap);
if !curvature.is_finite() || curvature <= 0.0 {
return Err(CgError::NotPositiveDefinite { iteration: k, curvature });
}
let alpha = rr / curvature;
for i in 0..n {
x[i] += alpha * p[i];
r[i] -= alpha * ap[i];
}
let rr_new = dot(&r, &r);
let residual = rr_new.sqrt();
if residual <= threshold {
return Ok((k + 1, residual));
}
let beta = rr_new / rr;
for i in 0..n {
p[i] = r[i] + beta * p[i];
}
rr = rr_new;
}
Err(CgError::NotConverged { max_iterations, residual_norm: rr.sqrt() })Run full test suite
cargo test --release 2>&1
Compiling cgsolve v0.5.1 (/workspace/cgsolve)
Finished `release` profile [optimized] target(s) in 7.84s
Running unittests src/lib.rs (target/release/deps/cgsolve-fea8032f53b8c6f5)
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-a7865db03c0aa40b)
running 27 tests
test already_solved_rhs_terminates_immediately ... ok
test converges_within_n_iterations ... ok
test duplicate_triplets_are_summed ... ok
test distinct_eigenvalues_converge_in_few_steps ... ok
test from_triplets_sums_dups_orders_columns_and_multiplies ... ok
test matvec_matches_dense_definition ... ok
test one_solver_serves_multiple_rhs ... ok
test rejects_dense_shape_mismatch_at_construction ... ok
test rejects_empty_matrix_at_construction ... ok
test rejects_indefinite_matrix ... ok
test rejects_negative_definite_matrix ... ok
test rejects_non_finite_entry_at_construction ... ok
test rejects_non_positive_tolerance ... ok
test rejects_out_of_bounds_triplet ... ok
test rejects_rhs_dimension_mismatch ... ok
test reports_non_convergence_with_tight_iteration_cap ... ok
test residual_is_tiny_for_spd_system ... ok
test residual_tiny_across_many_rhs ... ok
test scaling_rhs_scales_solution ... ok
test solution_is_linear_in_rhs ... ok
test solves_diagonal_system ... ok
test solves_known_2x2_system ... ok
test solves_large_laplacian_via_constructed_rhs ... ok
test solves_one_by_one_system ... ok
test solves_via_constructed_rhs ... ok
test symmetry_check_distinguishes ... ok
test three_distinct_eigenvalues_converge_in_at_most_three ... ok
test result: ok. 27 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Doc-tests cgsolve
running 8 tests
test src/config.rs - config::Config (line 21) ... ok
test src/cg.rs - cg::ConjugateGradient (line 27) ... ok
test src/lib.rs - (line 14) ... ok
test src/lib.rs - (line 30) ... ok
test src/matrix.rs - matrix::SparseMatrix::from_triplets (line 35) ... ok
test src/lib.rs - solve_spd (line 62) ... ok
test src/outcome.rs - outcome::CgOutcome::solution (line 31) ... ok
test src/matrix.rs - matrix::SparseMatrix::matvec (line 155) ... ok
test result: ok. 8 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.47s
[stdout]
Compiling cgsolve v0.5.1 (/workspace/cgsolve)
Finished `release` profile [optimized] target(s) in 7.84s
Running unittests src/lib.rs (target/release/deps/cgsolve-fea8032f53b8c6f5)
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-a7865db03c0aa40b)
running 27 tests
test already_solved_rhs_terminates_immediately ... ok
test converges_within_n_iterations ... ok
test duplicate_triplets_are_summed ... ok
test distinct_eigenvalues_converge_in_few_steps ... ok
test from_triplets_sums_dups_orders_columns_and_multiplies ... ok
test matvec_matches_dense_definition ... ok
test one_solver_serves_multiple_rhs ... ok
test rejects_dense_shape_mismatch_at_construction ... ok
test rejects_empty_matrix_at_construction ... ok
test rejects_indefinite_matrix ... ok
test rejects_negative_definite_matrix ... ok
test rejects_non_finite_entry_at_construction ... ok
test rejects_non_positive_tolerance ... ok
test rejects_out_of_bounds_triplet ... ok
test rejects_rhs_dimension_mismatch ... ok
test reports_non_convergence_with_tight_iteration_cap ... ok
test residual_is_tiny_for_spd_system ... ok
test residual_tiny_across_many_rhs ... ok
test scaling_rhs_scales_solution ... ok
test solution_is_linear_in_rhs ... ok
test solves_diagonal_system ... ok
test solves_known_2x2_system ... ok
test solves_large_laplacian_via_constructed_rhs ... ok
test solves_one_by_one_system ... ok
test solves_via_constructed_rhs ... ok
test symmetry_check_distinguishes ... ok
test three_distinct_eigenvalues_converge_in_at_most_three ... ok
test result: ok. 27 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Doc-tests cgsolve
running 8 tests
test src/config.rs - config::Config (line 21) ... ok
test src/cg.rs - cg::ConjugateGradient (line 27) ... ok
test src/lib.rs - (line 14) ... ok
test src/lib.rs - (line 30) ... ok
test src/matrix.rs - matrix::SparseMatrix::from_triplets (line 35) ... ok
test src/lib.rs - solve_spd (line 62) ... ok
test src/outcome.rs - outcome::CgOutcome::solution (line 31) ... ok
test src/matrix.rs - matrix::SparseMatrix::matvec (line 155) ... ok
test result: ok. 8 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.47sNo file changes, the agent left the workspace untouched.
[test] running cargo test --release (offline)
Finished `release` profile [optimized] target(s) in 0.01s
Running unittests src/lib.rs (target/release/deps/cgsolve-fea8032f53b8c6f5)
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-a7865db03c0aa40b)
running 27 tests
test converges_within_n_iterations ... ok
test distinct_eigenvalues_converge_in_few_steps ... ok
test already_solved_rhs_terminates_immediately ... ok
test duplicate_triplets_are_summed ... ok
test matvec_matches_dense_definition ... ok
test from_triplets_sums_dups_orders_columns_and_multiplies ... ok
test one_solver_serves_multiple_rhs ... ok
test rejects_dense_shape_mismatch_at_construction ... ok
test rejects_empty_matrix_at_construction ... ok
test rejects_indefinite_matrix ... ok
test rejects_negative_definite_matrix ... ok
test rejects_non_finite_entry_at_construction ... ok
test rejects_non_positive_tolerance ... ok
test rejects_out_of_bounds_triplet ... ok
test rejects_rhs_dimension_mismatch ... ok
test reports_non_convergence_with_tight_iteration_cap ... ok
test residual_is_tiny_for_spd_system ... ok
test residual_tiny_across_many_rhs ... ok
test scaling_rhs_scales_solution ... ok
test solves_diagonal_system ... ok
test solution_is_linear_in_rhs ... ok
test solves_known_2x2_system ... ok
test solves_one_by_one_system ... ok
test solves_large_laplacian_via_constructed_rhs ... ok
test solves_via_constructed_rhs ... ok
test symmetry_check_distinguishes ... ok
test three_distinct_eigenvalues_converge_in_at_most_three ... ok
test result: ok. 27 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
Finished `release` profile [optimized] target(s) in 0.01s
Doc-tests cgsolve
running 8 tests
test src/cg.rs - cg::ConjugateGradient (line 27) ... ok
test src/config.rs - config::Config (line 21) ... ok
test src/lib.rs - (line 14) ... ok
test src/lib.rs - (line 30) ... ok
test src/lib.rs - solve_spd (line 62) ... ok
test src/matrix.rs - matrix::SparseMatrix::from_triplets (line 35) ... ok
test src/matrix.rs - matrix::SparseMatrix::matvec (line 155) ... ok
test src/outcome.rs - outcome::CgOutcome::solution (line 31) ... ok
test result: ok. 8 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.51s
[test] all tests passed
[test] reward = 1Reproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_92c34045b5f445c8. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_92c34045b5f445c8 · verifier authoritative; classifier explanatory.