Module cg 

Conjugate Gradient solver for symmetric positive-definite systems.

Solves Ax = b where A is a symmetric positive-definite (SPD) sparse matrix in CSR format. The algorithm converges in at most n iterations for an n x n system, but in practice converges in O(sqrt(kappa) * log(1/eps)) iterations where kappa = cond(A).

Algorithm

Implements the Hestenes-Stiefel variant of Conjugate Gradient:

r = b - A*x
z = M^{-1} * r        (preconditioner; z = r when disabled)
p = z
rz = r . z

for k in 0..max_iterations:
    Ap = A * p
    alpha = rz / (p . Ap)
    x  = x + alpha * p
    r  = r - alpha * Ap
    if ||r||_2 < tolerance * ||b||_2:
        converged; break
    z  = M^{-1} * r
    rz_new = r . z
    beta = rz_new / rz
    p  = z + beta * p
    rz = rz_new

Preconditioning

When use_preconditioner is true, a diagonal (Jacobi) preconditioner is applied: M = diag(A), so that z_i = r_i / A_{ii}. This reduces the effective condition number for diagonally-dominant systems without adding significant per-iteration cost.

Numerical precision

All dot products and norm computations use f64 accumulation even though the matrix may store f32 values, preventing catastrophic cancellation in the inner products that drive the CG recurrence.

Convergence

Theoretical bound: ||x_k - x*||_A <= 2 * ((sqrt(kappa) - 1)/(sqrt(kappa) + 1))^k * ||x_0 - x*||_A where kappa is the 2-condition number of A.

Structs

ConjugateGradientSolver

Conjugate Gradient solver for symmetric positive-definite sparse systems.

Functions

axpy

Compute y[i] += alpha * x[i] for all i (AXPY operation, f32).

dot_product_f64

Compute the dot product of two f32 slices using f64 accumulation.

norm2

Compute the L2 norm of an f32 slice using f64 accumulation.