SyncVals · Trajectory

cg-solver

claude-code claude-opus-4-8 ✓ resolved GOOD_SUCCESS ↑ View task

Solved from the instruction alone, tests/ and solution/ were withheld from the agent's workspace and restored only for grading.

Reward = tests/test.sh exit code (0 → resolved); the classification below is post-hoc and cannot change it.

Classification , post-hoc; cannot change the reward

GOOD_SUCCESSLegitimate solve, implements the asked-for behavior; tests verify real functionality.

SubtypeLegitimate Solution

EvidenceTest output shows all 56 conformance tests passed, 2 smoke tests passed, and 8 doctests passed. The agent implemented `cg_iterate` correctly with: (1) initial residual r = b - A x, (2) proper threshold handling for relative residual with absolute fallback, (3) correct CG iteration loop with α and β updates per algorithm, (4) curvature non-finite/non-positive check returning NotPositiveDefinite, (5) convergence detection and iteration counting. The implementation matches the reference solution in src/cg.rs - both compute initial residual, compute threshold, check early termination, iterate with matvec per step, check curvature guard, and return appropriate errors.

Root causeThe task was well-specified with clear mathematical algorithm (conjugate-gradient iteration with explicit formulas), detailed error handling requirements, and comprehensive contract specification. The agent correctly understood the algorithm from instruction.md and implemented it faithfully, passing all 66 tests.

RecommendationN/A - task is fine

Trajectory

Tool-by-tool agent trajectory

16 tool calls · 3 tool types · 18 steps

# Implement the conjugate-gradient iteration in `cgsolve` ## Context The `cgsolve` crate provides a sparse linear solver for symmetric-positive-definite (SPD) systems `A x = b` via the conjugate-gradient (CG) method, for our scientific-computing stack. It lives at `/workspace/cgsolve` in this environment. The public API (the compressed-sparse-row `SparseMatrix` type, `Config`, the `ConjugateGradient` solver, the `CgOutcome` result type, the error enum, the doctests, and the full integration-test suite) is already in place. The crate compiles, but the **core iteration routine is unimplemented**, so the tests fail (the stub calls `todo!()`). Your job is to implement those routines so the existing test suite passes. ## Where it lives Six functions are stubbed with `todo!()`; implement all of them. Their signatures, the `SparseMatrix` fields, the `CgError` variants, and the other methods are visible in the crate , read them. 1. `src/matrix.rs` → `from_triplets` , assemble the matrix from `(row, col, value)` triplets (see *The sparse matrix* below). `from_dense` delegates to it. 2. `src/matrix.rs` → `matvec_into` , the allocation-free sparse product `out = A · x`; `matvec` wraps it. 3. `src/matrix.rs` → `get` , the stored value at `(row, col)`, or `0` if the entry is not present. 4. `src/matrix.rs` → `is_symmetric` , whether `A` equals its transpose to within an absolute tolerance. 5. `src/cg.rs` → `ConjugateGradient::solve` , the public entry: validate the configuration and the right-hand side, run the iteration from a zero start, and assemble the `CgOutcome` (see *The solver*). 6. `src/cg.rs` → `cg_iterate` , the numerical core of the iteration. Do **not** change the public API, the function signatures, the other modules (`config.rs`, `error.rs`, `outcome.rs`, `lib.rs`), or the tests. The module-private `dot(u, v)` (Euclidean dot product) is available. ## The sparse matrix `SparseMatrix` keeps only structurally nonzero entries of a square `n × n` operator in the compressed-sparse-row layout its fields describe; build that state in `from_triplets` and read it in the accessors and the product so they all agree. * **`from_triplets`** validates and packs the triplets: reject an empty matrix, an out-of-range coordinate, and a non-finite value with the matching `CgError` variant (the variants name their own fields). Coordinates that repeat the same `(row, col)` accumulate into a single entry, and each row's entries are kept in a canonical order so the accessors and product are deterministic. * **`matvec_into`** applies the operator, `out = A · x`, in time proportional to the number of stored entries, rejecting a length-mismatched `x` or `out`. ## The solver Solve a symmetric-positive-definite system `A x = b` with an iterative method that touches `A` only through the matrix–vector product , **exactly one product per step** , and is efficient on large sparse operators. The grader pins the *rate* an exact solver of this kind achieves (it must reach tolerance in no more steps than the system's dimension, and far fewer when the spectrum is clustered), so a slowly-converging choice will not pass; pick the method that meets those bounds. It starts from the initial guess in `x`, overwrites it with the solution, and returns `(iterations, residual_norm)`. The behavioral details , what residual measure the convergence test uses and how it scales, how an already-solved start and the per-step iteration count are reported, and which error each failure produces , follow the usual conventions for this method; match them exactly. The relevant cases: * Convergence is judged on the residual `b − A x` relative to the right-hand side, with a sensible fallback when the right-hand side vanishes. * The returned iteration count is the number of matrix–vector products that ran to completion (a start already within tolerance runs none, and the count increments only after a product's update is applied). * A step whose curvature `pᵀ A p` is non-positive or non-finite means the operator is not SPD , surface this as `CgError::NotPositiveDefinite` with the offending curvature and the step's index. The curvature is examined before that step does any work, so the step never completes and is not counted among the iterations performed: its index is how many products had already finished when it was reached. * Exhausting `max_iterations` without converging is `CgError::NotConverged { max_iterations, residual_norm }`. ## Contract / correctness requirements The implementation must satisfy (all covered by `cargo test`): - **Known systems.** Diagonal systems, a hand-checked 2×2, a `1×1`, and `b = A x_true` round-trips recover the solution; one solver serves many right-hand sides. - **Large sparse system.** A 50×50 tridiagonal discrete Laplacian `tridiag(-1, 2, -1)` with a constructed right-hand side is solved to good accuracy, exercising the sparse `matvec` and finite termination. - **Residual invariant.** For every solved system the residual `‖A x - b‖` is ~0 (below `1e-9`), across several right-hand sides; the reported `residual_norm` agrees with the recomputed value. - **Convergence guarantee.** CG reaches tolerance within `n` iterations for a well-conditioned system, within `2` iterations when the matrix has only two distinct eigenvalues, and in `0` iterations when the zero start already solves the system. - **Linearity.** The solve is linear in the right-hand side: `x(b1 + b2) = x(b1) + x(b2)` and `x(c·b) = c·x(b)`. - **Non-SPD detection.** Indefinite and negative-definite matrices are rejected with `NotPositiveDefinite` at the iteration where the curvature first goes non-positive. Your routine is responsible for this check. - **Non-convergence.** A demanding tolerance with a one-iteration cap reports `NotConverged` with a finite residual norm. - **Robustness.** Empty, dense-shape-mismatched, out-of-bounds, and non-finite matrix inputs, non-positive tolerance, and right-hand-side dimension mismatches are reported as the corresponding `CgError` (most handled by the surrounding code; your routine returns `NotPositiveDefinite` / `NotConverged`). ## Build & test ``` cd /workspace/cgsolve cargo test --release # unit + integration + doctests must all pass ``` The crate must build with no warnings under the existing lints (`#![forbid(unsafe_code)]`, `#![warn(missing_docs)]`).

List source files and workspace contents

shell

find /workspace/cgsolve/src -type f | sort && echo "---" && ls /workspace/cgsolve

output

/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

contents

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	        let _ = (row, col);
92	        todo!("implement get (sci-4519)")
93	    }
94	
95	    /// Whether the matrix is symmetric to within absolute tolerance `tol`,
96	    /// i.e. `|A[i][j] - A[j][i]| <= tol` for all `i, j`.
97	    pub fn is_symmetric(&self, tol: f64) -> bool {
98	        let _ = tol;
99	        todo!("implement is_symmetric (sci-4519)")
100	    }
101	
102	    /// Compute the matrix–vector product `A * x` into a fresh vector.
103	    ///
104	    /// Returns [`CgError::DimensionMismatch`] if `x.len() != dim()`.
105	    ///
106	    /// ```
107	    /// use cgsolve::SparseMatrix;
108	    /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
109	    /// let y = a.matvec(&[1.0, 1.0]).unwrap();
110	    /// assert_eq!(y, vec![2.0, 3.0]);
111	    /// ```
112	    pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, CgError> {
113	        if x.len() != self.n {
114	            return Err(CgError::DimensionMismatch {
115	                expected: self.n,
116	                got: x.len(),
117	            });
118	        }
119	        let mut out = vec![0.0; self.n];
120	        self.matvec_into(x, &mut out)
121	            .expect("output buffer sized to n");
122	        Ok(out)
123	    }
124	
125	    /// Compute `A * x`, writing the result into the preallocated `out` buffer.
126	    ///
127	    /// This is the allocation-free product used inside the conjugate-gradient
128	    /// iteration. Returns [`CgError::DimensionMismatch`] if either `x` or `out`
129	    /// has the wrong length.
130	    pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
131	        let _ = (x, out);
132	        todo!("implement matvec_into (sci-4519)")
133	    }
134	}
135

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	        let _ = b;
72	        todo!("implement solve (sci-4519)")
73	    }
74	}
75	
76	/// Euclidean dot product of two equal-length slices.
77	#[inline]
78	pub(crate) fn dot(u: &[f64], v: &[f64]) -> f64 {
79	    u.iter().zip(v).map(|(a, b)| a * b).sum()
80	}
81	
82	/// Numerical core of the conjugate-gradient method.
83	///
84	/// Solves `A x = b` in place: `x` holds the initial guess on entry (the public
85	/// wrapper passes the zero vector) and the solution on successful return. The
86	/// matrix `A` is assumed symmetric-positive-definite.
87	///
88	/// Returns `(iterations, residual_norm)` , the number of iterations performed
89	/// and the Euclidean norm of the final residual `b - A x`.
90	///
91	/// The exact stopping rule, iteration accounting, the non-SPD curvature
92	/// guard, and the error returns are specified at the crate level and pinned
93	/// by `tests/integration.rs`. The module-private `dot` helper is available;
94	/// all work is `O(nnz)` per iteration.
95	pub fn cg_iterate(
96	    a: &SparseMatrix,
97	    b: &[f64],
98	    x: &mut [f64],
99	    tolerance: f64,
100	    max_iterations: usize,
101	) -> Result<(usize, f64), CgError> {
102	    // TODO(sci-4519): implement the conjugate-gradient iteration described in
103	    // the doc comment above. Form the initial residual r = b - A x, iterate the
104	    // alpha/x/r/beta/p updates using one matvec per step, stop on the relative
105	    // residual threshold, return NotPositiveDefinite on non-positive curvature
106	    // and NotConverged if the budget is exhausted. See `tests/integration.rs`
107	    // for the contract under test. The `dot` helper below computes the
108	    // Euclidean dot products (rᵀr, pᵀAp, ...) you will need.
109	    let _ = (a, b, x, tolerance, max_iterations, dot);
110	    todo!("implement cg_iterate (sci-4519)")
111	}
112

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

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

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

1	//! # cgsolve
2	//!
3	//! Sparse linear solver for symmetric-positive-definite (SPD) systems via the
4	//! conjugate-gradient (CG) method.
5	//!
6	//! The crate exposes a compressed-sparse-row [`SparseMatrix`] type, a
7	//! [`Config`] holding the stopping criteria, and the [`ConjugateGradient`]
8	//! solver which solves `A x = b` for any number of right-hand sides using one
9	//! matrix–vector product per iteration. Failure modes (malformed matrix,
10	//! mismatched right-hand side, a non-SPD operator, or failure to converge) are
11	//! reported as a [`CgError`]; a successful solve returns a [`CgOutcome`] with
12	//! the solution and diagnostics.
13	//!
14	//! ```
15	//! use cgsolve::{Config, ConjugateGradient, SparseMatrix};
16	//!
17	//! // A symmetric-positive-definite system A x = b with A = [[4, 1], [1, 3]].
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 out = cg.solve(&[1.0, 2.0]).unwrap();
21	//!
22	//! // Residual b - A x is tiny.
23	//! let ax = a.matvec(out.solution()).unwrap();
24	//! assert!((ax[0] - 1.0).abs() < 1e-9);
25	//! assert!((ax[1] - 2.0).abs() < 1e-9);
26	//! ```
27	//!
28	//! A non-SPD operator is detected rather than silently producing nonsense:
29	//!
30	//! ```
31	//! use cgsolve::{Config, ConjugateGradient, CgError, SparseMatrix};
32	//! // Symmetric but indefinite (eigenvalues ±1): CG breaks down on a
33	//! // right-hand side along the negative-eigenvalue direction.
34	//! let a = SparseMatrix::from_dense(2, &[0.0, 1.0, 1.0, 0.0]).unwrap();
35	//! let cg = ConjugateGradient::new(&a, Config::new());
36	//! let err = cg.solve(&[1.0, -1.0]).unwrap_err();
37	//! assert!(matches!(err, CgError::NotPositiveDefinite { .. }));
38	//! ```
39	//!
40	//! See the [`cg`] module for the algorithm.
41	
42	#![forbid(unsafe_code)]
43	#![warn(missing_docs)]
44	
45	mod cg;
46	mod config;
47	mod error;
48	mod matrix;
49	mod outcome;
50	
51	pub use cg::{cg_iterate, ConjugateGradient};
52	pub use config::{Config, DEFAULT_MAX_ITER_FACTOR, DEFAULT_TOLERANCE};
53	pub use error::CgError;
54	pub use matrix::SparseMatrix;
55	pub use outcome::CgOutcome;
56	
57	/// Crate-level convenience: solve the SPD system `A x = b` with the default
58	/// [`Config`] in one call.
59	///
60	/// Equivalent to `ConjugateGradient::new(a, Config::new()).solve(b)`.
61	///
62	/// ```
63	/// use cgsolve::{solve_spd, SparseMatrix};
64	/// // 2x2 SPD system with the diagonal solution.
65	/// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 4.0]).unwrap();
66	/// let out = solve_spd(&a, &[6.0, 8.0]).unwrap();
67	/// assert!((out.solution[0] - 3.0).abs() < 1e-9);
68	/// assert!((out.solution[1] - 2.0).abs() < 1e-9);
69	/// ```
70	pub fn solve_spd(a: &SparseMatrix, b: &[f64]) -> Result<CgOutcome, CgError> {
71	    ConjugateGradient::new(a, Config::new()).solve(b)
72	}
73

/workspace/cgsolve/src/cg.rs

contents

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	        let _ = (row, col);
92	        todo!("implement get (sci-4519)")
93	    }
94	
95	    /// Whether the matrix is symmetric to within absolute tolerance `tol`,
96	    /// i.e. `|A[i][j] - A[j][i]| <= tol` for all `i, j`.
97	    pub fn is_symmetric(&self, tol: f64) -> bool {
98	        let _ = tol;
99	        todo!("implement is_symmetric (sci-4519)")
100	    }
101	
102	    /// Compute the matrix–vector product `A * x` into a fresh vector.
103	    ///
104	    /// Returns [`CgError::DimensionMismatch`] if `x.len() != dim()`.
105	    ///
106	    /// ```
107	    /// use cgsolve::SparseMatrix;
108	    /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
109	    /// let y = a.matvec(&[1.0, 1.0]).unwrap();
110	    /// assert_eq!(y, vec![2.0, 3.0]);
111	    /// ```
112	    pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, CgError> {
113	        if x.len() != self.n {
114	            return Err(CgError::DimensionMismatch {
115	                expected: self.n,
116	                got: x.len(),
117	            });
118	        }
119	        let mut out = vec![0.0; self.n];
120	        self.matvec_into(x, &mut out)
121	            .expect("output buffer sized to n");
122	        Ok(out)
123	    }
124	
125	    /// Compute `A * x`, writing the result into the preallocated `out` buffer.
126	    ///
127	    /// This is the allocation-free product used inside the conjugate-gradient
128	    /// iteration. Returns [`CgError::DimensionMismatch`] if either `x` or `out`
129	    /// has the wrong length.
130	    pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
131	        let _ = (x, out);
132	        todo!("implement matvec_into (sci-4519)")
133	    }
134	}
135

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	        let _ = b;
72	        todo!("implement solve (sci-4519)")
73	    }
74	}
75	
76	/// Euclidean dot product of two equal-length slices.
77	#[inline]
78	pub(crate) fn dot(u: &[f64], v: &[f64]) -> f64 {
79	    u.iter().zip(v).map(|(a, b)| a * b).sum()
80	}
81	
82	/// Numerical core of the conjugate-gradient method.
83	///
84	/// Solves `A x = b` in place: `x` holds the initial guess on entry (the public
85	/// wrapper passes the zero vector) and the solution on successful return. The
86	/// matrix `A` is assumed symmetric-positive-definite.
87	///
88	/// Returns `(iterations, residual_norm)` , the number of iterations performed
89	/// and the Euclidean norm of the final residual `b - A x`.
90	///
91	/// The exact stopping rule, iteration accounting, the non-SPD curvature
92	/// guard, and the error returns are specified at the crate level and pinned
93	/// by `tests/integration.rs`. The module-private `dot` helper is available;
94	/// all work is `O(nnz)` per iteration.
95	pub fn cg_iterate(
96	    a: &SparseMatrix,
97	    b: &[f64],
98	    x: &mut [f64],
99	    tolerance: f64,
100	    max_iterations: usize,
101	) -> Result<(usize, f64), CgError> {
102	    // TODO(sci-4519): implement the conjugate-gradient iteration described in
103	    // the doc comment above. Form the initial residual r = b - A x, iterate the
104	    // alpha/x/r/beta/p updates using one matvec per step, stop on the relative
105	    // residual threshold, return NotPositiveDefinite on non-positive curvature
106	    // and NotConverged if the budget is exhausted. See `tests/integration.rs`
107	    // for the contract under test. The `dot` helper below computes the
108	    // Euclidean dot products (rᵀr, pᵀAp, ...) you will need.
109	    let _ = (a, b, x, tolerance, max_iterations, dot);
110	    todo!("implement cg_iterate (sci-4519)")
111	}
112

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

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

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

1	//! # cgsolve
2	//!
3	//! Sparse linear solver for symmetric-positive-definite (SPD) systems via the
4	//! conjugate-gradient (CG) method.
5	//!
6	//! The crate exposes a compressed-sparse-row [`SparseMatrix`] type, a
7	//! [`Config`] holding the stopping criteria, and the [`ConjugateGradient`]
8	//! solver which solves `A x = b` for any number of right-hand sides using one
9	//! matrix–vector product per iteration. Failure modes (malformed matrix,
10	//! mismatched right-hand side, a non-SPD operator, or failure to converge) are
11	//! reported as a [`CgError`]; a successful solve returns a [`CgOutcome`] with
12	//! the solution and diagnostics.
13	//!
14	//! ```
15	//! use cgsolve::{Config, ConjugateGradient, SparseMatrix};
16	//!
17	//! // A symmetric-positive-definite system A x = b with A = [[4, 1], [1, 3]].
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 out = cg.solve(&[1.0, 2.0]).unwrap();
21	//!
22	//! // Residual b - A x is tiny.
23	//! let ax = a.matvec(out.solution()).unwrap();
24	//! assert!((ax[0] - 1.0).abs() < 1e-9);
25	//! assert!((ax[1] - 2.0).abs() < 1e-9);
26	//! ```
27	//!
28	//! A non-SPD operator is detected rather than silently producing nonsense:
29	//!
30	//! ```
31	//! use cgsolve::{Config, ConjugateGradient, CgError, SparseMatrix};
32	//! // Symmetric but indefinite (eigenvalues ±1): CG breaks down on a
33	//! // right-hand side along the negative-eigenvalue direction.
34	//! let a = SparseMatrix::from_dense(2, &[0.0, 1.0, 1.0, 0.0]).unwrap();
35	//! let cg = ConjugateGradient::new(&a, Config::new());
36	//! let err = cg.solve(&[1.0, -1.0]).unwrap_err();
37	//! assert!(matches!(err, CgError::NotPositiveDefinite { .. }));
38	//! ```
39	//!
40	//! See the [`cg`] module for the algorithm.
41	
42	#![forbid(unsafe_code)]
43	#![warn(missing_docs)]
44	
45	mod cg;
46	mod config;
47	mod error;
48	mod matrix;
49	mod outcome;
50	
51	pub use cg::{cg_iterate, ConjugateGradient};
52	pub use config::{Config, DEFAULT_MAX_ITER_FACTOR, DEFAULT_TOLERANCE};
53	pub use error::CgError;
54	pub use matrix::SparseMatrix;
55	pub use outcome::CgOutcome;
56	
57	/// Crate-level convenience: solve the SPD system `A x = b` with the default
58	/// [`Config`] in one call.
59	///
60	/// Equivalent to `ConjugateGradient::new(a, Config::new()).solve(b)`.
61	///
62	/// ```
63	/// use cgsolve::{solve_spd, SparseMatrix};
64	/// // 2x2 SPD system with the diagonal solution.
65	/// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 4.0]).unwrap();
66	/// let out = solve_spd(&a, &[6.0, 8.0]).unwrap();
67	/// assert!((out.solution[0] - 3.0).abs() < 1e-9);
68	/// assert!((out.solution[1] - 2.0).abs() < 1e-9);
69	/// ```
70	pub fn solve_spd(a: &SparseMatrix, b: &[f64]) -> Result<CgOutcome, CgError> {
71	    ConjugateGradient::new(a, Config::new()).solve(b)
72	}
73

/workspace/cgsolve/src/error.rs

contents

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	        let _ = (row, col);
92	        todo!("implement get (sci-4519)")
93	    }
94	
95	    /// Whether the matrix is symmetric to within absolute tolerance `tol`,
96	    /// i.e. `|A[i][j] - A[j][i]| <= tol` for all `i, j`.
97	    pub fn is_symmetric(&self, tol: f64) -> bool {
98	        let _ = tol;
99	        todo!("implement is_symmetric (sci-4519)")
100	    }
101	
102	    /// Compute the matrix–vector product `A * x` into a fresh vector.
103	    ///
104	    /// Returns [`CgError::DimensionMismatch`] if `x.len() != dim()`.
105	    ///
106	    /// ```
107	    /// use cgsolve::SparseMatrix;
108	    /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
109	    /// let y = a.matvec(&[1.0, 1.0]).unwrap();
110	    /// assert_eq!(y, vec![2.0, 3.0]);
111	    /// ```
112	    pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, CgError> {
113	        if x.len() != self.n {
114	            return Err(CgError::DimensionMismatch {
115	                expected: self.n,
116	                got: x.len(),
117	            });
118	        }
119	        let mut out = vec![0.0; self.n];
120	        self.matvec_into(x, &mut out)
121	            .expect("output buffer sized to n");
122	        Ok(out)
123	    }
124	
125	    /// Compute `A * x`, writing the result into the preallocated `out` buffer.
126	    ///
127	    /// This is the allocation-free product used inside the conjugate-gradient
128	    /// iteration. Returns [`CgError::DimensionMismatch`] if either `x` or `out`
129	    /// has the wrong length.
130	    pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
131	        let _ = (x, out);
132	        todo!("implement matvec_into (sci-4519)")
133	    }
134	}
135

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	        let _ = b;
72	        todo!("implement solve (sci-4519)")
73	    }
74	}
75	
76	/// Euclidean dot product of two equal-length slices.
77	#[inline]
78	pub(crate) fn dot(u: &[f64], v: &[f64]) -> f64 {
79	    u.iter().zip(v).map(|(a, b)| a * b).sum()
80	}
81	
82	/// Numerical core of the conjugate-gradient method.
83	///
84	/// Solves `A x = b` in place: `x` holds the initial guess on entry (the public
85	/// wrapper passes the zero vector) and the solution on successful return. The
86	/// matrix `A` is assumed symmetric-positive-definite.
87	///
88	/// Returns `(iterations, residual_norm)` , the number of iterations performed
89	/// and the Euclidean norm of the final residual `b - A x`.
90	///
91	/// The exact stopping rule, iteration accounting, the non-SPD curvature
92	/// guard, and the error returns are specified at the crate level and pinned
93	/// by `tests/integration.rs`. The module-private `dot` helper is available;
94	/// all work is `O(nnz)` per iteration.
95	pub fn cg_iterate(
96	    a: &SparseMatrix,
97	    b: &[f64],
98	    x: &mut [f64],
99	    tolerance: f64,
100	    max_iterations: usize,
101	) -> Result<(usize, f64), CgError> {
102	    // TODO(sci-4519): implement the conjugate-gradient iteration described in
103	    // the doc comment above. Form the initial residual r = b - A x, iterate the
104	    // alpha/x/r/beta/p updates using one matvec per step, stop on the relative
105	    // residual threshold, return NotPositiveDefinite on non-positive curvature
106	    // and NotConverged if the budget is exhausted. See `tests/integration.rs`
107	    // for the contract under test. The `dot` helper below computes the
108	    // Euclidean dot products (rᵀr, pᵀAp, ...) you will need.
109	    let _ = (a, b, x, tolerance, max_iterations, dot);
110	    todo!("implement cg_iterate (sci-4519)")
111	}
112

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

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

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

1	//! # cgsolve
2	//!
3	//! Sparse linear solver for symmetric-positive-definite (SPD) systems via the
4	//! conjugate-gradient (CG) method.
5	//!
6	//! The crate exposes a compressed-sparse-row [`SparseMatrix`] type, a
7	//! [`Config`] holding the stopping criteria, and the [`ConjugateGradient`]
8	//! solver which solves `A x = b` for any number of right-hand sides using one
9	//! matrix–vector product per iteration. Failure modes (malformed matrix,
10	//! mismatched right-hand side, a non-SPD operator, or failure to converge) are
11	//! reported as a [`CgError`]; a successful solve returns a [`CgOutcome`] with
12	//! the solution and diagnostics.
13	//!
14	//! ```
15	//! use cgsolve::{Config, ConjugateGradient, SparseMatrix};
16	//!
17	//! // A symmetric-positive-definite system A x = b with A = [[4, 1], [1, 3]].
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 out = cg.solve(&[1.0, 2.0]).unwrap();
21	//!
22	//! // Residual b - A x is tiny.
23	//! let ax = a.matvec(out.solution()).unwrap();
24	//! assert!((ax[0] - 1.0).abs() < 1e-9);
25	//! assert!((ax[1] - 2.0).abs() < 1e-9);
26	//! ```
27	//!
28	//! A non-SPD operator is detected rather than silently producing nonsense:
29	//!
30	//! ```
31	//! use cgsolve::{Config, ConjugateGradient, CgError, SparseMatrix};
32	//! // Symmetric but indefinite (eigenvalues ±1): CG breaks down on a
33	//! // right-hand side along the negative-eigenvalue direction.
34	//! let a = SparseMatrix::from_dense(2, &[0.0, 1.0, 1.0, 0.0]).unwrap();
35	//! let cg = ConjugateGradient::new(&a, Config::new());
36	//! let err = cg.solve(&[1.0, -1.0]).unwrap_err();
37	//! assert!(matches!(err, CgError::NotPositiveDefinite { .. }));
38	//! ```
39	//!
40	//! See the [`cg`] module for the algorithm.
41	
42	#![forbid(unsafe_code)]
43	#![warn(missing_docs)]
44	
45	mod cg;
46	mod config;
47	mod error;
48	mod matrix;
49	mod outcome;
50	
51	pub use cg::{cg_iterate, ConjugateGradient};
52	pub use config::{Config, DEFAULT_MAX_ITER_FACTOR, DEFAULT_TOLERANCE};
53	pub use error::CgError;
54	pub use matrix::SparseMatrix;
55	pub use outcome::CgOutcome;
56	
57	/// Crate-level convenience: solve the SPD system `A x = b` with the default
58	/// [`Config`] in one call.
59	///
60	/// Equivalent to `ConjugateGradient::new(a, Config::new()).solve(b)`.
61	///
62	/// ```
63	/// use cgsolve::{solve_spd, SparseMatrix};
64	/// // 2x2 SPD system with the diagonal solution.
65	/// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 4.0]).unwrap();
66	/// let out = solve_spd(&a, &[6.0, 8.0]).unwrap();
67	/// assert!((out.solution[0] - 3.0).abs() < 1e-9);
68	/// assert!((out.solution[1] - 2.0).abs() < 1e-9);
69	/// ```
70	pub fn solve_spd(a: &SparseMatrix, b: &[f64]) -> Result<CgOutcome, CgError> {
71	    ConjugateGradient::new(a, Config::new()).solve(b)
72	}
73

/workspace/cgsolve/src/outcome.rs

contents

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	        let _ = (row, col);
92	        todo!("implement get (sci-4519)")
93	    }
94	
95	    /// Whether the matrix is symmetric to within absolute tolerance `tol`,
96	    /// i.e. `|A[i][j] - A[j][i]| <= tol` for all `i, j`.
97	    pub fn is_symmetric(&self, tol: f64) -> bool {
98	        let _ = tol;
99	        todo!("implement is_symmetric (sci-4519)")
100	    }
101	
102	    /// Compute the matrix–vector product `A * x` into a fresh vector.
103	    ///
104	    /// Returns [`CgError::DimensionMismatch`] if `x.len() != dim()`.
105	    ///
106	    /// ```
107	    /// use cgsolve::SparseMatrix;
108	    /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
109	    /// let y = a.matvec(&[1.0, 1.0]).unwrap();
110	    /// assert_eq!(y, vec![2.0, 3.0]);
111	    /// ```
112	    pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, CgError> {
113	        if x.len() != self.n {
114	            return Err(CgError::DimensionMismatch {
115	                expected: self.n,
116	                got: x.len(),
117	            });
118	        }
119	        let mut out = vec![0.0; self.n];
120	        self.matvec_into(x, &mut out)
121	            .expect("output buffer sized to n");
122	        Ok(out)
123	    }
124	
125	    /// Compute `A * x`, writing the result into the preallocated `out` buffer.
126	    ///
127	    /// This is the allocation-free product used inside the conjugate-gradient
128	    /// iteration. Returns [`CgError::DimensionMismatch`] if either `x` or `out`
129	    /// has the wrong length.
130	    pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
131	        let _ = (x, out);
132	        todo!("implement matvec_into (sci-4519)")
133	    }
134	}
135

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	        let _ = b;
72	        todo!("implement solve (sci-4519)")
73	    }
74	}
75	
76	/// Euclidean dot product of two equal-length slices.
77	#[inline]
78	pub(crate) fn dot(u: &[f64], v: &[f64]) -> f64 {
79	    u.iter().zip(v).map(|(a, b)| a * b).sum()
80	}
81	
82	/// Numerical core of the conjugate-gradient method.
83	///
84	/// Solves `A x = b` in place: `x` holds the initial guess on entry (the public
85	/// wrapper passes the zero vector) and the solution on successful return. The
86	/// matrix `A` is assumed symmetric-positive-definite.
87	///
88	/// Returns `(iterations, residual_norm)` , the number of iterations performed
89	/// and the Euclidean norm of the final residual `b - A x`.
90	///
91	/// The exact stopping rule, iteration accounting, the non-SPD curvature
92	/// guard, and the error returns are specified at the crate level and pinned
93	/// by `tests/integration.rs`. The module-private `dot` helper is available;
94	/// all work is `O(nnz)` per iteration.
95	pub fn cg_iterate(
96	    a: &SparseMatrix,
97	    b: &[f64],
98	    x: &mut [f64],
99	    tolerance: f64,
100	    max_iterations: usize,
101	) -> Result<(usize, f64), CgError> {
102	    // TODO(sci-4519): implement the conjugate-gradient iteration described in
103	    // the doc comment above. Form the initial residual r = b - A x, iterate the
104	    // alpha/x/r/beta/p updates using one matvec per step, stop on the relative
105	    // residual threshold, return NotPositiveDefinite on non-positive curvature
106	    // and NotConverged if the budget is exhausted. See `tests/integration.rs`
107	    // for the contract under test. The `dot` helper below computes the
108	    // Euclidean dot products (rᵀr, pᵀAp, ...) you will need.
109	    let _ = (a, b, x, tolerance, max_iterations, dot);
110	    todo!("implement cg_iterate (sci-4519)")
111	}
112

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

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

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

1	//! # cgsolve
2	//!
3	//! Sparse linear solver for symmetric-positive-definite (SPD) systems via the
4	//! conjugate-gradient (CG) method.
5	//!
6	//! The crate exposes a compressed-sparse-row [`SparseMatrix`] type, a
7	//! [`Config`] holding the stopping criteria, and the [`ConjugateGradient`]
8	//! solver which solves `A x = b` for any number of right-hand sides using one
9	//! matrix–vector product per iteration. Failure modes (malformed matrix,
10	//! mismatched right-hand side, a non-SPD operator, or failure to converge) are
11	//! reported as a [`CgError`]; a successful solve returns a [`CgOutcome`] with
12	//! the solution and diagnostics.
13	//!
14	//! ```
15	//! use cgsolve::{Config, ConjugateGradient, SparseMatrix};
16	//!
17	//! // A symmetric-positive-definite system A x = b with A = [[4, 1], [1, 3]].
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 out = cg.solve(&[1.0, 2.0]).unwrap();
21	//!
22	//! // Residual b - A x is tiny.
23	//! let ax = a.matvec(out.solution()).unwrap();
24	//! assert!((ax[0] - 1.0).abs() < 1e-9);
25	//! assert!((ax[1] - 2.0).abs() < 1e-9);
26	//! ```
27	//!
28	//! A non-SPD operator is detected rather than silently producing nonsense:
29	//!
30	//! ```
31	//! use cgsolve::{Config, ConjugateGradient, CgError, SparseMatrix};
32	//! // Symmetric but indefinite (eigenvalues ±1): CG breaks down on a
33	//! // right-hand side along the negative-eigenvalue direction.
34	//! let a = SparseMatrix::from_dense(2, &[0.0, 1.0, 1.0, 0.0]).unwrap();
35	//! let cg = ConjugateGradient::new(&a, Config::new());
36	//! let err = cg.solve(&[1.0, -1.0]).unwrap_err();
37	//! assert!(matches!(err, CgError::NotPositiveDefinite { .. }));
38	//! ```
39	//!
40	//! See the [`cg`] module for the algorithm.
41	
42	#![forbid(unsafe_code)]
43	#![warn(missing_docs)]
44	
45	mod cg;
46	mod config;
47	mod error;
48	mod matrix;
49	mod outcome;
50	
51	pub use cg::{cg_iterate, ConjugateGradient};
52	pub use config::{Config, DEFAULT_MAX_ITER_FACTOR, DEFAULT_TOLERANCE};
53	pub use error::CgError;
54	pub use matrix::SparseMatrix;
55	pub use outcome::CgOutcome;
56	
57	/// Crate-level convenience: solve the SPD system `A x = b` with the default
58	/// [`Config`] in one call.
59	///
60	/// Equivalent to `ConjugateGradient::new(a, Config::new()).solve(b)`.
61	///
62	/// ```
63	/// use cgsolve::{solve_spd, SparseMatrix};
64	/// // 2x2 SPD system with the diagonal solution.
65	/// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 4.0]).unwrap();
66	/// let out = solve_spd(&a, &[6.0, 8.0]).unwrap();
67	/// assert!((out.solution[0] - 3.0).abs() < 1e-9);
68	/// assert!((out.solution[1] - 2.0).abs() < 1e-9);
69	/// ```
70	pub fn solve_spd(a: &SparseMatrix, b: &[f64]) -> Result<CgOutcome, CgError> {
71	    ConjugateGradient::new(a, Config::new()).solve(b)
72	}
73

/workspace/cgsolve/src/config.rs

contents

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	        let _ = (row, col);
92	        todo!("implement get (sci-4519)")
93	    }
94	
95	    /// Whether the matrix is symmetric to within absolute tolerance `tol`,
96	    /// i.e. `|A[i][j] - A[j][i]| <= tol` for all `i, j`.
97	    pub fn is_symmetric(&self, tol: f64) -> bool {
98	        let _ = tol;
99	        todo!("implement is_symmetric (sci-4519)")
100	    }
101	
102	    /// Compute the matrix–vector product `A * x` into a fresh vector.
103	    ///
104	    /// Returns [`CgError::DimensionMismatch`] if `x.len() != dim()`.
105	    ///
106	    /// ```
107	    /// use cgsolve::SparseMatrix;
108	    /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
109	    /// let y = a.matvec(&[1.0, 1.0]).unwrap();
110	    /// assert_eq!(y, vec![2.0, 3.0]);
111	    /// ```
112	    pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, CgError> {
113	        if x.len() != self.n {
114	            return Err(CgError::DimensionMismatch {
115	                expected: self.n,
116	                got: x.len(),
117	            });
118	        }
119	        let mut out = vec![0.0; self.n];
120	        self.matvec_into(x, &mut out)
121	            .expect("output buffer sized to n");
122	        Ok(out)
123	    }
124	
125	    /// Compute `A * x`, writing the result into the preallocated `out` buffer.
126	    ///
127	    /// This is the allocation-free product used inside the conjugate-gradient
128	    /// iteration. Returns [`CgError::DimensionMismatch`] if either `x` or `out`
129	    /// has the wrong length.
130	    pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
131	        let _ = (x, out);
132	        todo!("implement matvec_into (sci-4519)")
133	    }
134	}
135

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	        let _ = b;
72	        todo!("implement solve (sci-4519)")
73	    }
74	}
75	
76	/// Euclidean dot product of two equal-length slices.
77	#[inline]
78	pub(crate) fn dot(u: &[f64], v: &[f64]) -> f64 {
79	    u.iter().zip(v).map(|(a, b)| a * b).sum()
80	}
81	
82	/// Numerical core of the conjugate-gradient method.
83	///
84	/// Solves `A x = b` in place: `x` holds the initial guess on entry (the public
85	/// wrapper passes the zero vector) and the solution on successful return. The
86	/// matrix `A` is assumed symmetric-positive-definite.
87	///
88	/// Returns `(iterations, residual_norm)` , the number of iterations performed
89	/// and the Euclidean norm of the final residual `b - A x`.
90	///
91	/// The exact stopping rule, iteration accounting, the non-SPD curvature
92	/// guard, and the error returns are specified at the crate level and pinned
93	/// by `tests/integration.rs`. The module-private `dot` helper is available;
94	/// all work is `O(nnz)` per iteration.
95	pub fn cg_iterate(
96	    a: &SparseMatrix,
97	    b: &[f64],
98	    x: &mut [f64],
99	    tolerance: f64,
100	    max_iterations: usize,
101	) -> Result<(usize, f64), CgError> {
102	    // TODO(sci-4519): implement the conjugate-gradient iteration described in
103	    // the doc comment above. Form the initial residual r = b - A x, iterate the
104	    // alpha/x/r/beta/p updates using one matvec per step, stop on the relative
105	    // residual threshold, return NotPositiveDefinite on non-positive curvature
106	    // and NotConverged if the budget is exhausted. See `tests/integration.rs`
107	    // for the contract under test. The `dot` helper below computes the
108	    // Euclidean dot products (rᵀr, pᵀAp, ...) you will need.
109	    let _ = (a, b, x, tolerance, max_iterations, dot);
110	    todo!("implement cg_iterate (sci-4519)")
111	}
112

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

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

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

1	//! # cgsolve
2	//!
3	//! Sparse linear solver for symmetric-positive-definite (SPD) systems via the
4	//! conjugate-gradient (CG) method.
5	//!
6	//! The crate exposes a compressed-sparse-row [`SparseMatrix`] type, a
7	//! [`Config`] holding the stopping criteria, and the [`ConjugateGradient`]
8	//! solver which solves `A x = b` for any number of right-hand sides using one
9	//! matrix–vector product per iteration. Failure modes (malformed matrix,
10	//! mismatched right-hand side, a non-SPD operator, or failure to converge) are
11	//! reported as a [`CgError`]; a successful solve returns a [`CgOutcome`] with
12	//! the solution and diagnostics.
13	//!
14	//! ```
15	//! use cgsolve::{Config, ConjugateGradient, SparseMatrix};
16	//!
17	//! // A symmetric-positive-definite system A x = b with A = [[4, 1], [1, 3]].
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 out = cg.solve(&[1.0, 2.0]).unwrap();
21	//!
22	//! // Residual b - A x is tiny.
23	//! let ax = a.matvec(out.solution()).unwrap();
24	//! assert!((ax[0] - 1.0).abs() < 1e-9);
25	//! assert!((ax[1] - 2.0).abs() < 1e-9);
26	//! ```
27	//!
28	//! A non-SPD operator is detected rather than silently producing nonsense:
29	//!
30	//! ```
31	//! use cgsolve::{Config, ConjugateGradient, CgError, SparseMatrix};
32	//! // Symmetric but indefinite (eigenvalues ±1): CG breaks down on a
33	//! // right-hand side along the negative-eigenvalue direction.
34	//! let a = SparseMatrix::from_dense(2, &[0.0, 1.0, 1.0, 0.0]).unwrap();
35	//! let cg = ConjugateGradient::new(&a, Config::new());
36	//! let err = cg.solve(&[1.0, -1.0]).unwrap_err();
37	//! assert!(matches!(err, CgError::NotPositiveDefinite { .. }));
38	//! ```
39	//!
40	//! See the [`cg`] module for the algorithm.
41	
42	#![forbid(unsafe_code)]
43	#![warn(missing_docs)]
44	
45	mod cg;
46	mod config;
47	mod error;
48	mod matrix;
49	mod outcome;
50	
51	pub use cg::{cg_iterate, ConjugateGradient};
52	pub use config::{Config, DEFAULT_MAX_ITER_FACTOR, DEFAULT_TOLERANCE};
53	pub use error::CgError;
54	pub use matrix::SparseMatrix;
55	pub use outcome::CgOutcome;
56	
57	/// Crate-level convenience: solve the SPD system `A x = b` with the default
58	/// [`Config`] in one call.
59	///
60	/// Equivalent to `ConjugateGradient::new(a, Config::new()).solve(b)`.
61	///
62	/// ```
63	/// use cgsolve::{solve_spd, SparseMatrix};
64	/// // 2x2 SPD system with the diagonal solution.
65	/// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 4.0]).unwrap();
66	/// let out = solve_spd(&a, &[6.0, 8.0]).unwrap();
67	/// assert!((out.solution[0] - 3.0).abs() < 1e-9);
68	/// assert!((out.solution[1] - 2.0).abs() < 1e-9);
69	/// ```
70	pub fn solve_spd(a: &SparseMatrix, b: &[f64]) -> Result<CgOutcome, CgError> {
71	    ConjugateGradient::new(a, Config::new()).solve(b)
72	}
73

/workspace/cgsolve/src/lib.rs

contents

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	        let _ = (row, col);
92	        todo!("implement get (sci-4519)")
93	    }
94	
95	    /// Whether the matrix is symmetric to within absolute tolerance `tol`,
96	    /// i.e. `|A[i][j] - A[j][i]| <= tol` for all `i, j`.
97	    pub fn is_symmetric(&self, tol: f64) -> bool {
98	        let _ = tol;
99	        todo!("implement is_symmetric (sci-4519)")
100	    }
101	
102	    /// Compute the matrix–vector product `A * x` into a fresh vector.
103	    ///
104	    /// Returns [`CgError::DimensionMismatch`] if `x.len() != dim()`.
105	    ///
106	    /// ```
107	    /// use cgsolve::SparseMatrix;
108	    /// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 3.0]).unwrap();
109	    /// let y = a.matvec(&[1.0, 1.0]).unwrap();
110	    /// assert_eq!(y, vec![2.0, 3.0]);
111	    /// ```
112	    pub fn matvec(&self, x: &[f64]) -> Result<Vec<f64>, CgError> {
113	        if x.len() != self.n {
114	            return Err(CgError::DimensionMismatch {
115	                expected: self.n,
116	                got: x.len(),
117	            });
118	        }
119	        let mut out = vec![0.0; self.n];
120	        self.matvec_into(x, &mut out)
121	            .expect("output buffer sized to n");
122	        Ok(out)
123	    }
124	
125	    /// Compute `A * x`, writing the result into the preallocated `out` buffer.
126	    ///
127	    /// This is the allocation-free product used inside the conjugate-gradient
128	    /// iteration. Returns [`CgError::DimensionMismatch`] if either `x` or `out`
129	    /// has the wrong length.
130	    pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
131	        let _ = (x, out);
132	        todo!("implement matvec_into (sci-4519)")
133	    }
134	}
135

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	        let _ = b;
72	        todo!("implement solve (sci-4519)")
73	    }
74	}
75	
76	/// Euclidean dot product of two equal-length slices.
77	#[inline]
78	pub(crate) fn dot(u: &[f64], v: &[f64]) -> f64 {
79	    u.iter().zip(v).map(|(a, b)| a * b).sum()
80	}
81	
82	/// Numerical core of the conjugate-gradient method.
83	///
84	/// Solves `A x = b` in place: `x` holds the initial guess on entry (the public
85	/// wrapper passes the zero vector) and the solution on successful return. The
86	/// matrix `A` is assumed symmetric-positive-definite.
87	///
88	/// Returns `(iterations, residual_norm)` , the number of iterations performed
89	/// and the Euclidean norm of the final residual `b - A x`.
90	///
91	/// The exact stopping rule, iteration accounting, the non-SPD curvature
92	/// guard, and the error returns are specified at the crate level and pinned
93	/// by `tests/integration.rs`. The module-private `dot` helper is available;
94	/// all work is `O(nnz)` per iteration.
95	pub fn cg_iterate(
96	    a: &SparseMatrix,
97	    b: &[f64],
98	    x: &mut [f64],
99	    tolerance: f64,
100	    max_iterations: usize,
101	) -> Result<(usize, f64), CgError> {
102	    // TODO(sci-4519): implement the conjugate-gradient iteration described in
103	    // the doc comment above. Form the initial residual r = b - A x, iterate the
104	    // alpha/x/r/beta/p updates using one matvec per step, stop on the relative
105	    // residual threshold, return NotPositiveDefinite on non-positive curvature
106	    // and NotConverged if the budget is exhausted. See `tests/integration.rs`
107	    // for the contract under test. The `dot` helper below computes the
108	    // Euclidean dot products (rᵀr, pᵀAp, ...) you will need.
109	    let _ = (a, b, x, tolerance, max_iterations, dot);
110	    todo!("implement cg_iterate (sci-4519)")
111	}
112

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

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

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

1	//! # cgsolve
2	//!
3	//! Sparse linear solver for symmetric-positive-definite (SPD) systems via the
4	//! conjugate-gradient (CG) method.
5	//!
6	//! The crate exposes a compressed-sparse-row [`SparseMatrix`] type, a
7	//! [`Config`] holding the stopping criteria, and the [`ConjugateGradient`]
8	//! solver which solves `A x = b` for any number of right-hand sides using one
9	//! matrix–vector product per iteration. Failure modes (malformed matrix,
10	//! mismatched right-hand side, a non-SPD operator, or failure to converge) are
11	//! reported as a [`CgError`]; a successful solve returns a [`CgOutcome`] with
12	//! the solution and diagnostics.
13	//!
14	//! ```
15	//! use cgsolve::{Config, ConjugateGradient, SparseMatrix};
16	//!
17	//! // A symmetric-positive-definite system A x = b with A = [[4, 1], [1, 3]].
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 out = cg.solve(&[1.0, 2.0]).unwrap();
21	//!
22	//! // Residual b - A x is tiny.
23	//! let ax = a.matvec(out.solution()).unwrap();
24	//! assert!((ax[0] - 1.0).abs() < 1e-9);
25	//! assert!((ax[1] - 2.0).abs() < 1e-9);
26	//! ```
27	//!
28	//! A non-SPD operator is detected rather than silently producing nonsense:
29	//!
30	//! ```
31	//! use cgsolve::{Config, ConjugateGradient, CgError, SparseMatrix};
32	//! // Symmetric but indefinite (eigenvalues ±1): CG breaks down on a
33	//! // right-hand side along the negative-eigenvalue direction.
34	//! let a = SparseMatrix::from_dense(2, &[0.0, 1.0, 1.0, 0.0]).unwrap();
35	//! let cg = ConjugateGradient::new(&a, Config::new());
36	//! let err = cg.solve(&[1.0, -1.0]).unwrap_err();
37	//! assert!(matches!(err, CgError::NotPositiveDefinite { .. }));
38	//! ```
39	//!
40	//! See the [`cg`] module for the algorithm.
41	
42	#![forbid(unsafe_code)]
43	#![warn(missing_docs)]
44	
45	mod cg;
46	mod config;
47	mod error;
48	mod matrix;
49	mod outcome;
50	
51	pub use cg::{cg_iterate, ConjugateGradient};
52	pub use config::{Config, DEFAULT_MAX_ITER_FACTOR, DEFAULT_TOLERANCE};
53	pub use error::CgError;
54	pub use matrix::SparseMatrix;
55	pub use outcome::CgOutcome;
56	
57	/// Crate-level convenience: solve the SPD system `A x = b` with the default
58	/// [`Config`] in one call.
59	///
60	/// Equivalent to `ConjugateGradient::new(a, Config::new()).solve(b)`.
61	///
62	/// ```
63	/// use cgsolve::{solve_spd, SparseMatrix};
64	/// // 2x2 SPD system with the diagonal solution.
65	/// let a = SparseMatrix::from_dense(2, &[2.0, 0.0, 0.0, 4.0]).unwrap();
66	/// let out = solve_spd(&a, &[6.0, 8.0]).unwrap();
67	/// assert!((out.solution[0] - 3.0).abs() < 1e-9);
68	/// assert!((out.solution[1] - 2.0).abs() < 1e-9);
69	/// ```
70	pub fn solve_spd(a: &SparseMatrix, b: &[f64]) -> Result<CgOutcome, CgError> {
71	    ConjugateGradient::new(a, Config::new()).solve(b)
72	}
73

List test files

shell

ls /workspace/cgsolve/tests/

output

smoke.rs

[stdout]
smoke.rs

/workspace/cgsolve/tests/smoke.rs

contents

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

/workspace/cgsolve/src/matrix.rs

− removed

    pub fn from_triplets(n: usize, triplets: &[(usize, usize, f64)]) -> Result<Self, CgError> {
        let _ = (n, triplets);
        todo!("implement from_triplets (sci-4519)")
    }

+ added

    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), then accumulate duplicate coordinates.
        let mut sorted = triplets.to_vec();
        sorted.sort_unstable_by_key(|&(r, c, _)| (r, c));

        let mut entries: Vec<(usize, usize, f64)> = Vec::with_capacity(sorted.len());
        for (row, col, val) in sorted {
            match entries.last_mut() {
                Some(last) if last.0 == row && last.1 == col => last.2 += val,
                _ => entries.push((row, col, val)),
            }
        }

        // Build CSR arrays.
        let mut row_ptr = vec![0usize; n + 1];
        for &(row, _, _) in &entries {
            row_ptr[row + 1] += 1;
        }
        for i in 0..n {
            row_ptr[i + 1] += row_ptr[i];
        }
        let col_idx: Vec<usize> = entries.iter().map(|&(_, c, _)| c).collect();
        let values: Vec<f64> = entries.iter().map(|&(_, _, v)| v).collect();

        Ok(Self { n, row_ptr, col_idx, values })
    }

/workspace/cgsolve/src/matrix.rs

− removed

    pub fn get(&self, row: usize, col: usize) -> f64 {
        let _ = (row, col);
        todo!("implement get (sci-4519)")
    }

+ added

    pub fn get(&self, row: usize, col: usize) -> f64 {
        assert!(row < self.n, "row index {row} out of bounds for dimension {}", self.n);
        assert!(col < self.n, "col index {col} out of bounds for dimension {}", self.n);
        let start = self.row_ptr[row];
        let end = self.row_ptr[row + 1];
        match self.col_idx[start..end].binary_search(&col) {
            Ok(pos) => self.values[start + pos],
            Err(_) => 0.0,
        }
    }

/workspace/cgsolve/src/matrix.rs

− removed

    pub fn is_symmetric(&self, tol: f64) -> bool {
        let _ = tol;
        todo!("implement is_symmetric (sci-4519)")
    }

+ added

    pub fn is_symmetric(&self, tol: f64) -> bool {
        for i in 0..self.n {
            let start = self.row_ptr[i];
            let end = self.row_ptr[i + 1];
            for k in start..end {
                let j = self.col_idx[k];
                if (self.values[k] - self.get(j, i)).abs() > tol {
                    return false;
                }
            }
        }
        true
    }

/workspace/cgsolve/src/matrix.rs

− removed

    pub fn matvec_into(&self, x: &[f64], out: &mut [f64]) -> Result<(), CgError> {
        let _ = (x, out);
        todo!("implement matvec_into (sci-4519)")
    }

+ added

    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 start = self.row_ptr[i];
            let end = self.row_ptr[i + 1];
            out[i] = self.col_idx[start..end]
                .iter()
                .zip(&self.values[start..end])
                .map(|(&j, &v)| v * x[j])
                .sum();
        }
        Ok(())
    }

/workspace/cgsolve/src/cg.rs

− removed

    pub fn solve(&self, b: &[f64]) -> Result<CgOutcome, CgError> {
        let _ = b;
        todo!("implement solve (sci-4519)")
    }

+ added

    pub fn solve(&self, b: &[f64]) -> Result<CgOutcome, CgError> {
        self.config.validate()?;
        let n = self.a.dim();
        if b.len() != n {
            return Err(CgError::DimensionMismatch { expected: n, got: b.len() });
        }
        let mut x = vec![0.0f64; n];
        let max_iterations = self.config.effective_max_iterations(n);
        let tolerance = self.config.tolerance();
        let (iterations, residual_norm) =
            cg_iterate(self.a, b, &mut x, tolerance, max_iterations)?;
        Ok(CgOutcome { solution: x, iterations, residual_norm, converged: true })
    }

/workspace/cgsolve/src/cg.rs

− removed

pub fn cg_iterate(
    a: &SparseMatrix,
    b: &[f64],
    x: &mut [f64],
    tolerance: f64,
    max_iterations: usize,
) -> Result<(usize, f64), CgError> {
    // 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)")
}

+ added

pub fn cg_iterate(
    a: &SparseMatrix,
    b: &[f64],
    x: &mut [f64],
    tolerance: f64,
    max_iterations: usize,
) -> Result<(usize, f64), CgError> {
    let n = a.dim();

    // Initial residual r = b - A x. (Caller passes x = 0, so r = b.)
    let mut ap = vec![0.0f64; n];
    a.matvec_into(x, &mut ap).expect("dimensions match");
    let mut r: Vec<f64> = (0..n).map(|i| b[i] - ap[i]).collect();

    let b_norm = dot(b, b).sqrt();
    // Relative tolerance with a fallback of absolute tolerance when ||b|| = 0.
    let threshold = if b_norm > 0.0 { tolerance * b_norm } else { tolerance };

    let mut r_dot_r = dot(&r, &r);
    let mut residual_norm = r_dot_r.sqrt();

    if residual_norm <= threshold {
        return Ok((0, residual_norm));
    }

    let mut p = r.clone();

    for k in 0..max_iterations {
        // One matrix–vector product per step.
        a.matvec_into(&p, &mut ap).expect("dimensions match");

        let curvature = dot(&p, &ap);
        if curvature <= 0.0 || !curvature.is_finite() {
            return Err(CgError::NotPositiveDefinite { iteration: k, curvature });
        }

        let alpha = r_dot_r / curvature;

        for i in 0..n {
            x[i] += alpha * p[i];
            r[i] -= alpha * ap[i];
        }

        let r_dot_r_new = dot(&r, &r);
        residual_norm = r_dot_r_new.sqrt();

        if residual_norm <= threshold {
            return Ok((k + 1, residual_norm));
        }

        let beta = r_dot_r_new / r_dot_r;
        for i in 0..n {
            p[i] = r[i] + beta * p[i];
        }

        r_dot_r = r_dot_r_new;
    }

    Err(CgError::NotConverged { max_iterations, residual_norm })
}

Run full test suite

shell

cargo test --release 2>&1

output

Compiling cgsolve v0.5.1 (/workspace/cgsolve)
    Finished `release` profile [optimized] target(s) in 8.04s
     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 solves_a_small_spd_system ... ok
test reports_rhs_dimension_mismatch ... 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/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 8.04s
     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 solves_a_small_spd_system ... ok
test reports_rhs_dimension_mismatch ... 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/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

All 10 tests pass (2 integration + 8 doctests) with no warnings. Here's a summary of what was implemented: **`src/matrix.rs`** - `from_triplets`: validates n>0, bounds, and finiteness; sorts triplets by `(row, col)` and sums duplicates; builds the standard CSR arrays (`row_ptr`, `col_idx`, `values`). - `matvec_into`: validates both buffer lengths, then computes `out[i] = Σ_j A[i,j] * x[j]` in O(nnz) using the CSR row slices. - `get`: binary-searches the column-index slice for the requested row, returning 0.0 if absent. - `is_symmetric`: for every stored entry `A[i,j]`, checks `|A[i,j] - A[j,i]| <= tol`. **`src/cg.rs`** - `ConjugateGradient::solve`: validates tolerance and rhs dimension, allocates a zero starting vector, calls `cg_iterate`, and wraps the result in `CgOutcome`. - `cg_iterate`: standard CG , computes initial residual `r = b - Ax`, checks convergence before the first matvec (returns 0 iterations if already solved), then loops: computes `Ap`, checks curvature `p·Ap` (returning `NotPositiveDefinite` with the step index if ≤0 or non-finite), applies the `alpha/x/r/beta/p` updates, and returns `NotConverged` if the budget is exhausted. Convergence is judged as `‖r‖ ≤ tol·‖b‖` with a fallback to `‖r‖ ≤ tol` when `‖b‖ = 0`.

Final diff , what the agent changed

No file changes, the agent left the workspace untouched.

Verifier output , the proof behind the reward

tests/test.sh

[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.31s
     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 56 tests
test another_indefinite_fails_after_a_good_step ... ok
test already_solved_rhs_terminates_immediately ... ok
test asymmetric_triplets_are_built_and_detected ... ok
test cg_specific_rate_on_spread_spectrum ... ok
test clustered_distinct_eigenvalues_terminate ... 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 exact_small_system_iteration_and_solution ... ok
test from_triplets_sums_dups_orders_columns_and_multiplies ... ok
test get_reads_every_cell_of_an_asymmetric_matrix ... ok
test get_returns_zero_for_unstored_entry ... ok
test identity_solves_in_one_iteration ... ok
test ill_conditioned_diagonal_still_solves ... ok
test is_symmetric_respects_tolerance ... ok
test many_duplicate_triplets_on_one_cell_sum ... ok
test matvec_into_rejects_wrong_output_length ... ok
test matvec_matches_dense_definition ... ok
test large_laplacian_100_converges ... ok
test matvec_of_zero_vector_is_zero ... ok
test n_distinct_eigenvalues_take_exactly_n ... ok
test negative_one_by_one_is_rejected_immediately ... ok
test non_spd_detected_at_a_later_step ... ok
test one_solver_serves_multiple_rhs ... ok
test out_of_order_triplets_with_gap_rows ... 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 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_2d_laplacian_25 ... 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 tiny_rhs_uses_relative_threshold ... ok
test unit_vector_rhs_picks_one_component ... ok
test warm_start_at_exact_solution_does_no_work ... ok
test warm_start_from_partial_guess_still_converges ... ok
test warm_start_partway_on_3x3 ... ok
test zero_matrix_is_not_positive_definite_at_step_zero ... ok
test zero_rhs_with_nonzero_start_drives_to_zero ... ok

test result: ok. 56 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

    Finished `release` profile [optimized] target(s) in 0.01s
   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/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 = 1

Reproduce this trial: git checkout 2f94510 && PYTHONPATH=src python3 scripts/build_site.py , then open trial/trial_1ebd2a0a133d4d65. Re-running the agent live requires EVAL_PLATFORM_ENABLE_OAUTH_SMOKE=1 and is non-deterministic.

Trial trial_1ebd2a0a133d4d65 · verifier authoritative; classifier explanatory.