Conjugate Gradient

Conjugate Gradient

Iterative CG solver for a symmetric positive-definite sparse linear system. Demonstrates GPU-resident buffers — once data lands on the device it never leaves until the algorithm converges.

Source: src/3.Maths/ConjugateGradient/Program.cs

Resident arrays

FloatResidentArray / IntResidentArray keep memory device-side. Explicit RefreshDevice() / RefreshHost() calls are the only transfers:

public static void ConjugateGradient(FloatResidentArray X, SparseMatrix A, FloatResidentArray B, int maxiter, float eps)
{
    A.RefreshDevice(); X.RefreshDevice(); B.RefreshDevice();
    // … iterate fully on device …
    X.RefreshHost();
}

Main loop

wrapper.Fmsub(R, B, A, X, N);              // R = B - A*X
wrapper.Copy(P, R, N);
int k = 0;
while (k < maxiter)
{
    wrapper.Multiply(AP, A, P, N);         // AP = A*P
    float r = ScalarProd(R, R, N);
    float alpha = r / ScalarProd(P, AP, N);
    wrapper.Saxpy(X, X, alpha, P, N);      // X += alpha*P
    wrapper.Saxpy(R, R, -alpha, AP, N);    // R -= alpha*AP
    float rr = ScalarProd(R, R, N);
    if (rr < eps * eps) break;

    float beta = rr / r;
    wrapper.Saxpy(P, R, beta, P, N);       // P = R + beta*P
    ++k;
}

Scalar products call into a small native helper (CUB under the hood) via P/Invoke. Within Hybridizer Essentials, custom CUB-backed reductions aren't exposed natively — see Reduction, GenericReduction, and LambdaReduction for hand-rolled alternatives that stay in pure C#.

The driver runs ~10 000-row Laplacians by default; convergence without a preconditioner is slow but the example focuses on data movement, not numerical performance.