Python API

The CPU calls are synchronous. Numerical breakdown and exhausted iterations return a SolveInfo status; invalid dimensions, layouts, options and Jacobi setup raise an exception. The stopping rule uses the independently evaluated true residual:

||b - A @ x||2 <= max(atol, rtol * ||b||2)

Symmetry and positive definiteness are caller preconditions. non_finite also reports required inner products outside the usable double range, even when the final reported norm is finite.

Convenience and strict functions

cgpy.solve(A, b, x0=None, *, num_threads=1, rtol=1e-08, atol=0.0, maxiter=1000, algorithm='standard', fused=True, residual_refresh=50, preconditioner='none', dense_kernel='native')

Solve a real SPD system and return an owned solution.

Parameters:

  • A (array_like or CsrMatrix) – Nonempty square operator. Dense input is normalized to aligned C-contiguous float64 storage. Symmetry and positive definiteness are caller preconditions.

  • b (array_like) – Real right-hand side of shape (n,), normalized to float64.

  • x0 (array_like, optional) – Real initial guess of shape (n,). Omission selects zeros. The returned solution always owns a separate writable array.

  • num_threads (int, default=1) – Positive native CPU thread count. Multiple threads require OpenMP; the BLAS kernel accepts only one thread.

  • rtol (float, default=1e-8, 0.0) – Finite nonnegative tolerances, with stopping threshold max(atol, rtol * ||b||2).

  • atol (float, default=1e-8, 0.0) – Finite nonnegative tolerances, with stopping threshold max(atol, rtol * ||b||2).

  • maxiter (int, default=1000) – Nonnegative update limit. Zero checks only the initial guess.

  • algorithm ({'standard', 'explicit'}, default='standard') – Residual recurrence or the explicit-residual reference formulation.

  • fused (bool, default=True) – Fuse supported CPU operations. This is independent of algorithm and dense kernel. BLAS retains the selected vector fusion policy.

  • residual_refresh (int, default=50) – Standard-CG true-residual replacement/restart interval. Zero disables periodic refresh, while candidate and final checks remain enabled.

  • preconditioner ({'none', 'jacobi'}, default='none') – Jacobi is prepared for this call and requires a finite positive diagonal/inverse and algorithm='standard'.

  • dense_kernel ({'native', 'blas'}, default='native') – Explicit dense matrix-vector kernel. BLAS requires an enabled build, dense input and num_threads=1. CSR requires 'native'.

Returns:

SolveResult – Owned solution and native status, true residual and operation counts.

Raises:

  • ValueError – Dimensions, real-valued conversion, options, build capabilities or Jacobi setup are invalid. Iteration limits must be integers >= 0.

  • TypeError – A native argument has an unsupported type or integer range.

Notes

Inputs are preserved. Compatible A and b arrays can be borrowed until the synchronous call returns; do not modify them concurrently. Native computation releases the GIL. Numerical failure and iteration limits are reported in status; always inspect converged.

Examples

>>> import numpy as np
>>> import cgpy
>>> result = cgpy.solve(
...     [[4., 1.], [1., 3.]], [1., 2.], preconditioner='jacobi'
... )
>>> result.converged
True
>>> np.testing.assert_allclose(result.x, [1 / 11, 7 / 11])
cgpy.prepare(A, *, num_threads=1, preconditioner='none', dense_kernel='native')

Snapshot a matrix and cache optional Jacobi setup for repeated solves.

Parameters:

  • A (array_like or CsrMatrix) – Nonempty square real SPD operator. Dense inputs are normalized to aligned C-contiguous float64 storage before the owning snapshot.

  • num_threads (int, default=1) – Native CPU threads. Values greater than one require an OpenMP build. The BLAS kernel accepts only one thread.

  • preconditioner ({'none', 'jacobi'}, default='none') – Jacobi requires a finite strictly positive diagonal with a finite positive inverse and supports standard CG only.

  • dense_kernel ({'native', 'blas'}, default='native') – Explicit dense matrix-vector kernel. BLAS requires an enabled build, dense input and num_threads=1. CSR requires 'native'.

Returns:

PreparedSolver – Owning immutable matrix/preconditioner snapshot. Thread count and kernel choice are fixed; each solve selects its tolerances and fusion. Create this type through prepare; it has no public constructor.

Raises:

  • ValueError – Dimensions, options, build capabilities or Jacobi setup are invalid.

  • TypeError – A native integer argument has an unsupported type or range.

Notes

Changes to the original matrix do not affect the prepared solver. Each overlapping solve needs a distinct output and create_workspace() result. Preparation allocates; native solver vectors are reused afterward.

Examples

>>> import numpy as np
>>> import cgpy
>>> plan = cgpy.prepare([[4., 1.], [1., 3.]], preconditioner='jacobi')
>>> x = np.zeros(2)
>>> info = plan.solve_into(
...     np.array([1., 2.]), x, workspace=plan.create_workspace()
... )
>>> info.converged
True
cgpy.solve_into(A, b, x, *, workspace, num_threads=1, rtol=1e-08, atol=0.0, maxiter=1000, algorithm='standard', fused=True, residual_refresh=50, preconditioner='none', dense_kernel='native')

Solve in place with strict arrays and caller-owned workspace.

Parameters:

  • A (numpy.ndarray or CsrMatrix) – Nonempty square real SPD operator. Dense arrays must already be aligned, C-contiguous and float64.

  • b (numpy.ndarray) – Aligned C-contiguous float64 right-hand side of shape (n,).

  • x (numpy.ndarray) – Writable aligned C-contiguous float64 vector of shape (n,). Supplies the initial guess and receives the current approximation.

  • workspace (Workspace) – Reusable vectors with matching size, unused by overlapping calls.

  • num_threads (int, default=1) – Positive native CPU thread count; multiple threads require OpenMP. BLAS requires one thread.

  • rtol (float, default=1e-8, 0.0) – Finite nonnegative tolerances. The true-residual threshold is max(atol, rtol * ||b||2).

  • atol (float, default=1e-8, 0.0) – Finite nonnegative tolerances. The true-residual threshold is max(atol, rtol * ||b||2).

  • maxiter (int, default=1000) – Nonnegative update limit. Zero checks only the initial guess.

  • algorithm ({'standard', 'explicit'}, default='standard') – Residual recurrence or the explicit-residual reference formulation.

  • fused (bool, default=True) – Fuse supported CPU operations. Independent of algorithm and dense kernel; BLAS retains the selected vector update/reduction policy.

  • residual_refresh (int, default=50) – Standard-CG true-residual replacement/restart interval. Zero disables periodic replacement; candidate and final checks remain enabled.

  • preconditioner ({'none', 'jacobi'}, default='none') – Construct Jacobi for this call; requires algorithm='standard'.

  • dense_kernel ({'native', 'blas'}, default='native') – Explicit dense matrix-vector kernel. BLAS requires an enabled build, dense input and num_threads=1; CSR requires 'native'.

Returns:

SolveInfo – Status and counts. Numerical breakdown or an iteration limit may leave an unconverged approximation in x.

Raises:

  • ValueError – Shapes, dtype/layout, aliasing, workspace use, options or Jacobi setup violate the contract. Nonintegral/negative Python iteration limits also raise this exception.

  • TypeError – A native argument has an unsupported type or integer range.

Notes

There is no implicit array conversion. Keep x disjoint from A and b. Borrowed arrays must remain alive and free from conflicting mutation until this synchronous call returns. Native computation releases the GIL. Use prepare to reuse matrix and Jacobi setup across calls.

Examples

>>> import numpy as np
>>> import cgpy
>>> A = np.array([[4., 1.], [1., 3.]])
>>> x = np.zeros(2)
>>> info = cgpy.solve_into(
...     A, np.array([1., 2.]), x, workspace=cgpy.Workspace(2)
... )
>>> info.status
'converged'

solve normalizes real array-like inputs and owns its output. The strict solve_into function and prepared methods require aligned C-contiguous float64 NumPy vectors, with a writable output. Dense strict inputs must also already have that dtype and layout. Output must be disjoint from the matrix and right-hand side. Keep borrowed arrays alive and free from conflicting mutation until the call returns.

The package ships native type stubs and a py.typed marker. CPU solutions and strict array inputs use numpy.typing.NDArray[numpy.float64]; convenience host inputs use ArrayLike. CUDA results contain an owned cuda.Vector. Finite option sets use Literal annotations, and iteration limits use SupportsIndex to retain NumPy integer support. External CUDA inputs use cuda.SupportsCudaArrayInterface; its descriptor remains typed as object for native validation. Array shape, layout, writeability, device, stream and aliasing still require the existing runtime validation; annotations do not establish those conditions.

dense_kernel="native" is the default. The optional "blas" choice is dense-only and accepts num_threads=1. fused independently controls the remaining CPU vector operations. A prepared solver fixes its matrix snapshot, thread count, preconditioner and dense-kernel choice when it is constructed. Dense and CSR inputs return the same PreparedSolver type. CsrMatrix is the public spelling in both CPU and CUDA APIs. Create CPU prepared solvers with cgpy.prepare; cgpy.PreparedSolver is the return type and has no public constructor. CUDA follows the same factory-only contract through cuda.prepare and cuda.prepare_host.

Repeated solves

>>> A = np.array([[4.0, 1.0], [1.0, 3.0]])
>>> b = np.array([1.0, 2.0])
>>> plan = cgpy.prepare(A, preconditioner="jacobi")
>>> workspace = plan.create_workspace()
>>> x = np.zeros_like(b)
>>> info = plan.solve_into(b, x, workspace=workspace, rtol=1e-10)
>>> info.converged
True
>>> np.testing.assert_allclose(x, [1 / 11, 7 / 11])
>>> x.fill(0.0)
>>> plan.solve_into(b, x, workspace=workspace, rtol=1e-10).converged
True

Each overlapping solve needs a distinct workspace and output. Prepared matrix and preconditioner snapshots can be shared. Long native computation releases the GIL; this does not make concurrent mutation of borrowed NumPy arrays safe.

class cgpy.PreparedSolver

An owning matrix snapshot with fixed thread, preconditioner and kernel settings.

Notes

Create through cgpy.prepare to normalize a dense input or wrap a public CsrMatrix. The prepared object owns its matrix and optional Jacobi inverse. Direct construction of this return type is unsupported. Each overlapping solve needs a separate workspace and output. All calls are synchronous, and native computation releases the GIL.

create_workspace(self: cgpy._core.PreparedSolver) → cgpy._core.Workspace

Allocate reusable solver vectors matching this system.

Returns:

Workspace – New independent workspace for one active solve at a time.

property dense_kernel

Explicit ‘native’ or ‘blas’ kernel selected at preparation.

Type:

str

property num_threads

CPU thread count fixed at preparation time.

Type:

int

property preconditioner

‘none’ or ‘jacobi’, fixed at preparation time.

Type:

str

property preconditioner_nbytes

Stored Jacobi inverse bytes; zero without Jacobi.

Type:

int

property preconditioner_setup_ns

Jacobi construction time in nanoseconds; zero without Jacobi.

Type:

int

property shape

Square operator dimensions.

Type:

tuple of int

property size

Matrix dimension and workspace size.

Type:

int

solve_into(self: cgpy._core.PreparedSolver, b: numpy.ndarray, x: numpy.ndarray, *, workspace: cgpy._core.Workspace, rtol: SupportsFloat = 1e-08, atol: SupportsFloat = 0.0, maxiter: object = 1000, algorithm: str = 'standard', fused: bool = True, residual_refresh: object = 50) → cgpy._core.SolveInfo

Solve in place using the prepared operator and caller-owned vectors.

Parameters:

  • b (numpy.ndarray) – Aligned C-contiguous float64 right-hand side of shape (size,).

  • x (numpy.ndarray) – Writable aligned C-contiguous float64 initial guess and output, disjoint from b. Reset it before repeating a solve with the same initial guess.

  • workspace (Workspace) – Matching workspace, unused by overlapping calls.

  • rtol (float, default=1e-8, 0.0) – Finite nonnegative tolerances. True-residual threshold is max(atol, rtol * ||b||2).

  • atol (float, default=1e-8, 0.0) – Finite nonnegative tolerances. True-residual threshold is max(atol, rtol * ||b||2).

  • maxiter (int, default=1000) – Maximum updates; zero checks only the initial guess.

  • algorithm ({'standard', 'explicit'}, default='standard') – Residual recurrence or explicit-residual formulation. Jacobi requires standard.

  • fused (bool, default=True) – Fuse supported CPU operations independently of the fixed dense kernel.

  • residual_refresh (int, default=50) – Standard-CG replacement/restart interval; zero disables periodic refresh.

Returns:

SolveInfo – Convergence/failure status, independently checked residual and counts.

Raises:

  • ValueError – Invalid arrays, options, workspace use, aliasing or algorithm/preconditioner. Iteration limits must implement the integer index protocol and fit a nonnegative native size_t; negative, nonintegral and overflowing counts name the invalid option in this exception.

  • TypeError – An array or workspace argument has an unsupported native type.

Notes

The matrix, thread count, dense kernel and preconditioner are fixed by preparation. Borrowed arrays must remain alive and free from conflicting mutation until the synchronous call returns. The GIL is released for the solve.

class cgpy.Workspace

Own reusable residual, direction, product, scratch and preconditioned vectors.

Parameters:

size (int) – Nonnegative system dimension. A solve requires a matching nonempty matrix.

Raises:

TypeError – Size is negative, out of native integer range or nonintegral.

Notes

One workspace can serve only one active call. Overlapping use raises ValueError; concurrent solves need independent workspaces and outputs. The solver reuses these five vectors instead of allocating them during each prepared solve.

Examples

>>> import cgpy
>>> workspace = cgpy.Workspace(2)
>>> workspace.size
2
property capacity

Allocated vector capacity, at least the current size.

Type:

int

property size

Current vector length, required to match the matrix dimension.

Type:

int

CSR and results

class cgpy.CsrMatrix(indptr, indices, data, shape)

Own a validated canonical compressed sparse row matrix.

Parameters:

  • indptr (array_like of int) – One-dimensional row pointers of length shape[0] + 1. They must start at zero, be nondecreasing, and end at the number of values.

  • indices (array_like of int) – One-dimensional nonnegative column indices, sorted and unique within each row. Each entry must be smaller than shape[1].

  • data (array_like of float) – One-dimensional real values matching indices. Values are copied into float64 storage.

  • shape (tuple of int) – Nonnegative row and column counts. Solving additionally requires a nonempty square matrix.

Raises:

  • ValueError – Index arrays, dimensions, canonical ordering or data shape are invalid.

  • TypeError – Native dimensions or indices cannot be represented by size_t.

Notes

Construction snapshots the arrays. Later changes to them do not alter the operator. Strided, unaligned, read-only and non-native-endian inputs are normalized when needed. Signed/unsigned integer indices are range-checked before native conversion. Canonicalize SciPy CSR with sum_duplicates() and sort_indices() before wrapping arrays that need those operations.

Examples

>>> import cgpy
>>> A = cgpy.CsrMatrix([0, 2, 4], [0, 1, 0, 1], [4., 1., 1., 3.], (2, 2))
>>> A.shape, A.nnz
((2, 2), 4)
property shape

Matrix dimensions (rows, columns).

Type:

tuple of int

property nnz

Number of stored entries, including explicitly stored zeros.

Type:

int

>>> A = cgpy.CsrMatrix([0, 2, 4], [0, 1, 0, 1], [4.0, 1.0, 1.0, 3.0], (2, 2))
>>> A.shape, A.nnz
((2, 2), 4)
>>> result = cgpy.solve(A, np.array([1.0, 2.0]), preconditioner="jacobi")
>>> result.converged
True
>>> result.restarts == result.info.restarts
True
>>> np.testing.assert_allclose(result.x, [1 / 11, 7 / 11])
>>> repeated = cgpy.solve(A, np.array([1.0, 2.0]), preconditioner="jacobi")
>>> result == repeated
False
>>> np.testing.assert_allclose(result.x, repeated.x)

CSR row pointers start at zero, are monotone, and end at nnz. Columns must be in range, sorted and unique within each row. Construction owns a copy of the arrays. Jacobi requires a stored finite strictly positive diagonal and a finite positive inverse. It supports algorithm="standard" only.

SolveResult fields cannot be rebound, but solution contents remain writable. Equality compares result object identity; compare numerical solutions explicitly with a suitable tolerance, as above. The iterations, matvecs and restarts properties expose the corresponding SolveInfo counts directly; they do not keep a second copy of the metadata.

class cgpy.SolveResult(x, info)

Owned solution and independently checked native convergence data.

x is a NumPy array for CPU solves or a CUDA Vector for device solves. Fields cannot be rebound; solution contents remain writable. Results compare by object identity; compare solution values explicitly. Check converged before using an approximation.

property converged

Whether the true residual satisfies the requested tolerance.

property status

Convergence, iteration-limit or numerical-failure status.

property iterations

Number of completed solution updates.

property matvecs

Matrix-vector products, including true-residual checks.

property restarts

Search-direction restarts after true-residual refreshes.

property residual_norm

Independently evaluated final Euclidean residual norm.

class cgpy.SolveInfo

Read-only numerical status and operation counts for a completed native call.

Variables:

  • status (str) – converged, max_iterations, non_positive_curvature or non_finite. Negative curvature has its own status; zero/subnormal required products are reported as outside the usable double range.

  • converged (bool) – Whether the independently evaluated true residual met the tolerance.

  • iterations (int) – Completed solution updates, excluding the initial residual check.

  • matvecs (int) – Matrix-vector products, including initial, refresh and final checks.

  • residual_norm (float) – Independently evaluated final Euclidean true-residual norm.

  • restarts (int) – Direction restarts after a true-residual replacement.

Notes

Numerical breakdown and iteration limits return this record instead of raising an exception. non_finite includes unusable inner-product underflow/overflow; the reported true residual can still be finite. Invalid inputs raise exceptions.

property converged

Whether the true-residual stopping rule was satisfied.

Type:

bool

property iterations

Number of completed solution updates.

Type:

int

property matvecs

Matrix-vector products, including true-residual checks.

Type:

int

property residual_norm

Independently evaluated final Euclidean residual norm.

Type:

float

property restarts

Direction restarts following residual replacement.

Type:

int

property status

Convergence, iteration-limit or numerical-failure status.

Type:

str

Build capabilities

cgpy.openmp_enabled() → bool

Report OpenMP support in the extension imported by this process.

Returns:

bool – Whether native kernel calls can request more than one CPU thread.

cgpy.blas_enabled() → bool

Report optional sequential BLAS support in the imported extension.

Returns:

bool – Whether dense_kernel='blas' is available for dense, one-thread calls.

Notes

This reports the actual imported binary. It does not change kernel selection or inspect another build directory. Native remains the default kernel.

cgpy.build_info() → dict

Describe the native extension imported by this Python process.

Returns:

dict – Compiler, build type, flags and OpenMP support. The blas dictionary records enabled, provider/configuration, integer width, thread mode and symbol prefix. Provider-specific fields are None without BLAS.

Notes

BLAS-enabled builds use the private sequential LP64 OpenBLAS provider. Build information belongs to this loaded binary, not another preset or installed wheel.

These report the extension actually imported by this Python process. They do not describe a different CMake build or a wheel installed elsewhere.

Iteration limits and migration

maxiter and residual_refresh follow the same contract for CPU and CUDA convenience, strict and prepared solves. They accept Python integers, NumPy integer scalars and objects implementing __index__. Python booleans retain their integer meaning of zero or one. Values must be nonnegative and fit the platform’s native size_t; floating-point values are not truncated. Invalid types, negative values and overflow raise ValueError naming the option. Exceptions raised by user-defined __index__ implementations, other than an invalid integer-conversion TypeError, are propagated.

This intentionally changes CPU prepared solves: callers that previously caught TypeError for an invalid iteration limit must catch ValueError. The accepted numerical options, iteration policy and native C++ CgOptions representation are unchanged. Convenience calls check limits before allocating their output/workspace; strict and prepared bindings use the same parser before entering the solver. Integer-protocol callbacks are evaluated once per option and complete before the binding validates and borrows array views; a callback cannot bypass output writeability checks by changing flags after an earlier descriptor was captured.

CSR conversion owns its native snapshot. Foreign dtype, byte order, alignment or strides may require an additional normalization buffer; the native snapshot does not borrow those buffers. NumPy arrays are copied through the buffer interface without converting each scalar through Python. Canonical row/column validation remains in the C++ matrix constructor.