Conjugate Gradient Method Library

cgpy solves real symmetric positive-definite systems A x = b with CG and Jacobi PCG. Its C++20 core and CPython 3.12 interface support float64 dense and canonical CSR matrices on CPU or an explicit optional CUDA backend.

First solve

>>> import numpy as np
>>> import cgpy
>>> A = np.array([[4.0, 1.0], [1.0, 3.0]])
>>> b = np.array([1.0, 2.0])
>>> result = cgpy.solve(A, b, rtol=1e-10, preconditioner="jacobi")
>>> result.status
'converged'
>>> np.testing.assert_allclose(result.x, [1 / 11, 7 / 11])
>>> bool(result.residual_norm <= 1e-10 * np.linalg.norm(b))
True

The convenience call preserves its inputs and returns an owned solution. See Installation for package installation.

API choices

CPU and CUDA provide corresponding solve, strict and prepared APIs. Choose based on what should be owned or reused:

API

Use case

Ownership and reuse

solve

One solve

Fresh output, workspace and any Jacobi setup per call.

solve_into

Caller-managed buffers

Borrow the matrix and RHS; use caller-owned output and workspace.

prepare(...).solve_into

Fixed matrix, repeated RHS

Own a matrix snapshot and cache optional Jacobi; use caller-owned output and workspace.

cuda.Session.solve

CUDA matrices that may change between calls

Reuse one fixed-size workspace; borrow each matrix and return a fresh output, with Jacobi rebuilt per call.

Each overlapping solve needs its own workspace and output. See Python API for strict CPU array contracts and CUDA for device data, preparation and stream lifetimes.

Architecture and limits

  • Storage, backend kernels, solver coordination and bindings have separate responsibilities; the native core builds without Python.

  • RAII owners manage storage and execution resources. Lightweight views borrow them, while prepared solvers and sessions expose explicit reuse choices.

  • Host and device execution share small numerical policies, including true-residual convergence checks, while retaining their own scheduling.

Symmetry and positive definiteness are caller preconditions; the library does not verify them through factorization. Arithmetic is float64. See Architecture for design decisions and Benchmarks for representative performance, measurement scope and remaining limitations.

Guides