tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it.List source files and workspace contents
find /workspace/cgsolve/src -type f | sort && echo "---" && ls /workspace/cgsolve/
/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 --- Cargo.lock Cargo.toml src target tests [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 --- Cargo.lock Cargo.toml src target tests
/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/src/config.rs
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
File does not exist. Note: your current working directory is /workspace/cgsolve.
[error] tool reported failure/workspace/cgsolve/tests/integration.rs
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
File does not exist. Note: your current working directory is /workspace/cgsolve.
[error] tool reported failureList test files
find /workspace/cgsolve/tests -type f | sort
/workspace/cgsolve/tests/smoke.rs [stdout] /workspace/cgsolve/tests/smoke.rs
/workspace/cgsolve/tests/smoke.rs
1 //! Trivial smoke test: the public API is exercised on one small SPD system.
2 //!
3 //! This shows the shape of the crate; the full behavioral conformance suite is
4 //! applied separately by the grader and is not shipped here.
5
6 use cgsolve::{solve_spd, CgError, Config, ConjugateGradient, SparseMatrix};
7
8 #[test]
9 fn solves_a_small_spd_system() {
10 let a = SparseMatrix::from_dense(2, &[4.0, 1.0, 1.0, 3.0]).expect("valid matrix");
11 let out = solve_spd(&a, &[1.0, 2.0]).expect("solve succeeds");
12 let r = a.matvec(out.solution()).unwrap();
13 assert!((r[0] - 1.0).abs() < 1e-9 && (r[1] - 2.0).abs() < 1e-9);
14 }
15
16 #[test]
17 fn reports_rhs_dimension_mismatch() {
18 let a = SparseMatrix::from_dense(2, &[4.0, 1.0, 1.0, 3.0]).unwrap();
19 let cg = ConjugateGradient::new(&a, Config::new());
20 let err = cg.solve(&[1.0, 2.0, 3.0]).unwrap_err();
21 assert!(matches!(err, CgError::DimensionMismatch { expected: 2, got: 3 }));
22 }
23
1 //! The result type returned by a successful conjugate-gradient solve.
2
3 /// The outcome of a converged conjugate-gradient solve.
4 ///
5 /// The crate distinguishes hard failures (reported as
6 /// [`CgError`](crate::CgError)) from a numerically completed solve. This struct
7 /// carries the solution together with diagnostics describing how it was
8 /// obtained.
9 #[derive(Debug, Clone, PartialEq)]
10 #[non_exhaustive]
11 pub struct CgOutcome {
12 /// The computed solution vector `x` of `A x = b`.
13 pub solution: Vec<f64>,
14
15 /// The number of conjugate-gradient iterations performed.
16 pub iterations: usize,
17
18 /// The Euclidean norm of the final residual `b - A x`.
19 pub residual_norm: f64,
20
21 /// Whether the iteration met the relative residual tolerance. Always `true`
22 /// for a value returned by [`solve`](crate::ConjugateGradient::solve)
23 /// (non-convergence is reported as an error instead); retained for callers
24 /// that inspect the diagnostics directly.
25 pub converged: bool,
26 }
27
28 impl CgOutcome {
29 /// Convenience accessor returning the solution vector by reference.
30 ///
31 /// ```
32 /// use cgsolve::{solve_spd, SparseMatrix};
33 /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 4.0]).unwrap();
34 /// let out = solve_spd(&a, &[6.0, 8.0]).unwrap();
35 /// assert_eq!(out.solution().len(), 2);
36 /// ```
37 pub fn solution(&self) -> &[f64] {
38 &self.solution
39 }
40 }
41
/workspace/cgsolve/src/outcome.rs
1 //! Trivial smoke test: the public API is exercised on one small SPD system.
2 //!
3 //! This shows the shape of the crate; the full behavioral conformance suite is
4 //! applied separately by the grader and is not shipped here.
5
6 use cgsolve::{solve_spd, CgError, Config, ConjugateGradient, SparseMatrix};
7
8 #[test]
9 fn solves_a_small_spd_system() {
10 let a = SparseMatrix::from_dense(2, &[4.0, 1.0, 1.0, 3.0]).expect("valid matrix");
11 let out = solve_spd(&a, &[1.0, 2.0]).expect("solve succeeds");
12 let r = a.matvec(out.solution()).unwrap();
13 assert!((r[0] - 1.0).abs() < 1e-9 && (r[1] - 2.0).abs() < 1e-9);
14 }
15
16 #[test]
17 fn reports_rhs_dimension_mismatch() {
18 let a = SparseMatrix::from_dense(2, &[4.0, 1.0, 1.0, 3.0]).unwrap();
19 let cg = ConjugateGradient::new(&a, Config::new());
20 let err = cg.solve(&[1.0, 2.0, 3.0]).unwrap_err();
21 assert!(matches!(err, CgError::DimensionMismatch { expected: 2, got: 3 }));
22 }
23
1 //! The result type returned by a successful conjugate-gradient solve.
2
3 /// The outcome of a converged conjugate-gradient solve.
4 ///
5 /// The crate distinguishes hard failures (reported as
6 /// [`CgError`](crate::CgError)) from a numerically completed solve. This struct
7 /// carries the solution together with diagnostics describing how it was
8 /// obtained.
9 #[derive(Debug, Clone, PartialEq)]
10 #[non_exhaustive]
11 pub struct CgOutcome {
12 /// The computed solution vector `x` of `A x = b`.
13 pub solution: Vec<f64>,
14
15 /// The number of conjugate-gradient iterations performed.
16 pub iterations: usize,
17
18 /// The Euclidean norm of the final residual `b - A x`.
19 pub residual_norm: f64,
20
21 /// Whether the iteration met the relative residual tolerance. Always `true`
22 /// for a value returned by [`solve`](crate::ConjugateGradient::solve)
23 /// (non-convergence is reported as an error instead); retained for callers
24 /// that inspect the diagnostics directly.
25 pub converged: bool,
26 }
27
28 impl CgOutcome {
29 /// Convenience accessor returning the solution vector by reference.
30 ///
31 /// ```
32 /// use cgsolve::{solve_spd, SparseMatrix};
33 /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 4.0]).unwrap();
34 /// let out = solve_spd(&a, &[6.0, 8.0]).unwrap();
35 /// assert_eq!(out.solution().len(), 2);
36 /// ```
37 pub fn solution(&self) -> &[f64] {
38 &self.solution
39 }
40 }
41
/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 });
}
}
// Bucket triplets by row, then sort each row by column and merge duplicates.
let mut row_entries: Vec<Vec<(usize, f64)>> = vec![Vec::new(); n];
for &(row, col, value) in triplets {
row_entries[row].push((col, value));
}
let mut row_ptr = Vec::with_capacity(n + 1);
let mut col_idx = Vec::new();
let mut values = Vec::new();
row_ptr.push(0usize);
for row_data in &mut row_entries {
row_data.sort_by_key(|&(c, _)| c);
let mut j = 0;
while j < row_data.len() {
let col = row_data[j].0;
let mut val = 0.0;
while j < row_data.len() && row_data[j].0 == col {
val += row_data[j].1;
j += 1;
}
col_idx.push(col);
values.push(val);
}
row_ptr.push(col_idx.len());
}
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 = a.dim();
// r = b - A*x (x is the zero vector on entry, so r = b; compute generally).
let mut r: Vec<f64> = {
let ax = a.matvec(x)?;
b.iter().zip(ax.iter()).map(|(&bi, &axi)| bi - axi).collect()
};
let b_norm = dot(b, b).sqrt();
let threshold = if b_norm == 0.0 { tolerance } else { tolerance * b_norm };
let mut rr = dot(&r, &r);
if rr.sqrt() <= threshold {
return Ok((0, rr.sqrt()));
}
let mut p = r.clone();
let mut ap = vec![0.0_f64; n];
for iter in 0..max_iterations {
a.matvec_into(&p, &mut ap)?;
let curvature = dot(&p, &ap);
if curvature <= 0.0 || !curvature.is_finite() {
return Err(CgError::NotPositiveDefinite { iteration: iter + 1, 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 r_norm = rr_new.sqrt();
if r_norm <= threshold {
return Ok((iter + 1, r_norm));
}
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 the test suite
cargo test --release 2>&1
Compiling cgsolve v0.5.1 (/workspace/cgsolve)
Finished `release` profile [optimized] target(s) in 6.91s
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/smoke.rs (target/release/deps/smoke-3b889307a456a1be)
running 2 tests
test reports_rhs_dimension_mismatch ... ok
test solves_a_small_spd_system ... ok
test result: ok. 2 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
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/outcome.rs - outcome::CgOutcome::solution (line 31) ... ok
test src/matrix.rs - matrix::SparseMatrix::matvec (line 160) ... ok
test result: ok. 8 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.48s
[stdout]
Compiling cgsolve v0.5.1 (/workspace/cgsolve)
Finished `release` profile [optimized] target(s) in 6.91s
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/smoke.rs (target/release/deps/smoke-3b889307a456a1be)
running 2 tests
test reports_rhs_dimension_mismatch ... ok
test solves_a_small_spd_system ... ok
test result: ok. 2 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
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/outcome.rs - outcome::CgOutcome::solution (line 31) ... ok
test src/matrix.rs - matrix::SparseMatrix::matvec (line 160) ... ok
test result: ok. 8 passed; 0 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.48sNo file changes, the agent left the workspace untouched.
[test] injected hidden conformance suite
[test] running cargo test --release (offline)
Compiling cgsolve v0.5.1 (/workspace/cgsolve)
Finished `release` profile [optimized] target(s) in 1.13s
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/conformance.rs (target/release/deps/conformance-044e8b6dc5ae63bd)
running 40 tests
test asymmetric_triplets_are_built_and_detected ... ok
test already_solved_rhs_terminates_immediately ... ok
test converges_within_n_iterations ... ok
test distinct_eigenvalues_converge_in_few_steps ... ok
test duplicate_triplets_are_summed ... ok
test duplicate_triplets_can_sum_to_a_stored_zero ... ok
test from_triplets_sums_dups_orders_columns_and_multiplies ... ok
test ill_conditioned_diagonal_still_solves ... ok
test identity_solves_in_one_iteration ... ok
test matvec_into_rejects_wrong_output_length ... ok
test matvec_matches_dense_definition ... ok
test n_distinct_eigenvalues_take_exactly_n ... ok
test large_laplacian_100_converges ... ok
test one_solver_serves_multiple_rhs ... ok
test negative_one_by_one_is_rejected_immediately ... FAILED
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 ... FAILED
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 relative_threshold_handles_huge_rhs ... 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 semidefinite_operator_with_zero_eigenvalue_is_rejected ... ok
test solution_is_linear_in_rhs ... ok
test solves_known_2x2_system ... ok
test solves_diagonal_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 warm_start_at_exact_solution_does_no_work ... ok
test warm_start_from_partial_guess_still_converges ... ok
test zero_matrix_is_not_positive_definite_at_step_zero ... FAILED
failures:
---- negative_one_by_one_is_rejected_immediately stdout ----
thread 'negative_one_by_one_is_rejected_immediately' panicked at tests/conformance.rs:480:13:
assertion `left == right` failed
left: 1
right: 0
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace
---- rejects_negative_definite_matrix stdout ----
thread 'rejects_negative_definite_matrix' panicked at tests/conformance.rs:264:13:
assertion `left == right` failed
left: 1
right: 0
---- zero_matrix_is_not_positive_definite_at_step_zero stdout ----
thread 'zero_matrix_is_not_positive_definite_at_step_zero' panicked at tests/conformance.rs:459:13:
assertion `left == right` failed
left: 1
right: 0
failures:
negative_one_by_one_is_rejected_immediately
rejects_negative_definite_matrix
zero_matrix_is_not_positive_definite_at_step_zero
test result: FAILED. 37 passed; 3 failed; 0 ignored; 0 measured; 0 filtered out; finished in 0.00s
error: test failed, to rerun pass `--test conformance`
[test] tests FAILED
[test] reward = 0Reproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_06be8428670d4663. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.
Trial trial_06be8428670d4663 · verifier authoritative; classifier explanatory.