A Pedagogical Introduction to Porting a Conjugate Gradient Solver to CUDA

GPU programming tutorials often begin with isolated examples: vector addition, reductions, matrix multiplication, and memory coalescing. These are useful for learning individual CUDA concepts, but there is another question that quickly arises when working with real scientific software:

How do we actually port an existing CPU application to a GPU?

In this article, I will walk through a small but realistic example using HPCCG, a Conjugate Gradient mini-application from the Mantevo project. Rather than rewriting the entire solver for CUDA at once, we will migrate it incrementally:

1. profile the original CPU application and identify the dominant kernel,
2. move sparse matrix-vector multiplication (SpMV) to the GPU,
3. migrate the remaining vector operations and reductions,
4. examine why accelerating individual kernels is not enough,
5. keep the CG working set resident on the GPU, and
6. profile the resulting implementation to see where the bottleneck moved.

As a final experiment, we will return to the original CPU implementation and port it independently to Kokkos, then run essentially the same Kokkos source using both CPU and CUDA backends. This gives us a simple way to compare an explicitly written CUDA implementation with a performance-portability approach.

The purpose is not to produce the fastest possible Conjugate Gradient implementation. Instead, the goal is to make the porting and performance-engineering process itself visible: profile, form a hypothesis, change one part of the application, measure again, and follow the bottleneck.

All of the code and intermediate CUDA-porting commits used in this article are available in my HPCCG GitHub repository, so the individual stages can be checked out and reproduced while following along.

1. The Starting Point: Vanilla HPCCG

HPCCG solves a linear system arising from a 27-point stencil on a structured 3D grid, using an unpreconditioned Conjugate Gradient method. The solver is deliberately simple, which makes it a good porting target: almost all of the work is done by three kernels.

• SpMV (HPC_sparsemv): computes Ap = A * p
• WAXPBY (waxpby): computes w = alpha * x + beta * y
• DDOT (ddot): computes the dot product x · y

One CG iteration strings these together roughly as follows:

for (int k = 1; k < max_iter && normr > tolerance; k++) {
  if (k == 1) {
    waxpby(nrow, 1.0, r, 0.0, r, p);          // p = r
  } else {
    oldrtrans = rtrans;
    rtrans = ddot(nrow, r, r);                // global reduction
    double beta = rtrans / oldrtrans;
    waxpby(nrow, 1.0, r, beta, p, p);         // p = r + beta * p
  }

  normr = sqrt(rtrans);                       // convergence check

  HPC_sparsemv(A, p, Ap);                     // Ap = A * p

  double alpha = ddot(nrow, p, Ap);           // global reduction
  alpha = rtrans / alpha;

  waxpby(nrow, 1.0, x,  alpha, p,  x);        // x = x + alpha * p
  waxpby(nrow, 1.0, r, -alpha, Ap, r);        // r = r - alpha * Ap
}

Keep this loop in mind. Every porting decision later in the article is really a decision about where each line of this loop runs and where its data lives.

Profiling the CPU baseline

Before writing a single line of CUDA, we need to know where the time goes. Conveniently, HPCCG already reports a per-kernel timing breakdown.

For a 100^3 problem, corresponding to one million rows, the original single-process CPU implementation took:

The solver required 149 iterations and reached a final residual of approximately 8.0 × 10⁻²¹.

SpMV is clearly the dominant operation, accounting for more than 80% of the CPU runtime. This is not particularly surprising. Each matrix row accesses neighboring elements through indexed loads, while performing relatively little arithmetic per byte moved.

This gives us our first hypothesis:

If SpMV dominates, moving SpMV to the GPU should be the highest-leverage first step.

Commit 51fdbe5 — Starting off CUDA migration -- sparsemv is hotspot

2. Moving SpMV to the GPU

The data structure problem

The first obstacle is not actually the kernel. It is the matrix representation.

HPCCG stores the matrix as arrays of per-row pointers:

for (int i = 0; i < nrow; i++) {
  double sum = 0.0;

  const double* const cur_vals = A->ptr_to_vals_in_row[i];
  const int* const cur_inds = A->ptr_to_inds_in_row[i];
  const int cur_nnz = A->nnz_in_row[i];

  for (int j = 0; j < cur_nnz; j++)
    sum += cur_vals[j] * x[cur_inds[j]];

  y[i] = sum;
}

This is convenient on the CPU, but copying the top-level structure to device memory would not make its host pointers magically become valid GPU pointers.

I therefore flattened the matrix into Compressed Sparse Row (CSR) form:

• row_offsets[nrow + 1]
• cols[nnz]
• vals[nnz]

Each array is contiguous and can be allocated and copied to the device independently.

Commit 05ace3b — sparsemv device copy and free routine

The first CUDA SpMV

The initial CUDA kernel deliberately uses the simplest possible mapping: one CUDA thread per matrix row.

__global__
void spmv_kernel(DeviceCSRMatrix A,
                 const double* x,
                 double* y)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if (i < A.nrow) {
        double sum = 0.0;

        for (int j = A.row_offsets[i];
             j < A.row_offsets[i + 1]; ++j)
            sum += A.vals[j] * x[A.cols[j]];

        y[i] = sum;
    }
}

At this point I intentionally did not optimize the sparse-memory access pattern. The goal was first to establish a correct GPU baseline.

The rest of CG still ran on the CPU, so each SpMV required the input vector to reach the GPU and the result to return to the CPU.

Commit ef5c262 — Add CUDA CSR SpMV path to HPCCG

This is an important intermediate stage. A GPU kernel can be much faster than its CPU counterpart while the application as a whole remains limited by everything surrounding that kernel.

3. Migrating WAXPBY and DDOT

WAXPBY

WAXPBY is the easiest operation to port. Each output element is independent:

__global__
void waxpby_kernel(int n,
                   double alpha, const double* x,
                   double beta,  const double* y,
                   double* w)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if (i < n)
        w[i] = alpha * x[i] + beta * y[i];
}

During the migration I first compared the CPU and GPU results before switching the actual solver call sites over.

Commits 4b78d30 — Added waxpby CUDA kernel and correctness check
870b4bf — Finished migrating WAXPBY

DDOT: now we need a reduction

DDOT is different. WAXPBY maps one input element to one output element; DDOT collapses an entire vector into a single scalar.

I implemented it as a simple two-stage reduction.

In the first kernel, each thread accumulates a local sum using a grid-stride loop. Threads within each block then reduce those values through shared memory, producing one partial result per block.

Conceptually:

x[i] * y[i]
      ↓
thread-local accumulation
      ↓
shared-memory block reduction
      ↓
one partial / block
      ↓
second reduction
      ↓
device scalar

A second kernel reduces the block-level partial results to the final device scalar.

I deliberately used a custom reduction rather than immediately calling a library routine such as CUB or cuBLAS. For this exercise, seeing the reduction and its synchronization structure explicitly was part of the point.

Commit 1101bf2 — DDOT also is on GPU

But this creates a new problem.

CG needs these scalar results to compute alpha, beta, and the convergence criterion. If the vector operations live on the GPU but every reduction result immediately returns to the CPU, the algorithm repeatedly forces the two processors to synchronize.

We have accelerated all three kernels, but we have not yet fixed the application.

4. Why Accelerating Kernels Is Not Enough

At this point, SpMV, WAXPBY, and DDOT all have CUDA implementations.

It is tempting to declare the port complete.

It isn't.

Nsight Systems showed that the hybrid implementation was repeatedly entering cudaMemcpy. In one intermediate profile I observed:

cudaMemcpy        304 calls
cudaMalloc         10 calls
cudaLaunchKernel 1195 calls

The CUDA API attributed roughly 405 ms to the cudaMemcpy calls, compared with roughly 29 ms of host API time spent launching kernels.

There is an important subtlety here: this does not mean that copying a few bytes over PCIe literally takes milliseconds.

A blocking cudaMemcpy is also a synchronization point. Its API duration can include time spent waiting for previously queued GPU work to complete.

So the important observation was not simply:

PCIe is slow.

It was:

The application is repeatedly forcing the CPU and GPU to synchronize.

This is a different level of performance problem.

We started with a kernel-level bottleneck: SpMV.

After accelerating the kernels, we exposed an application-dataflow bottleneck.

5. Keeping the CG Working Set on the GPU

The next optimization therefore does not change the mathematics of SpMV, DDOT, or WAXPBY at all.

It changes where the data lives.

The main CG working set consists of:

x
r
p
Ap
A

In the GPU-resident version, these vectors and the CSR matrix are allocated on the GPU before entering the iteration loop. The iterative solver then operates directly on device pointers.

Conceptually:

// Allocate and initialize the GPU working set once.

spmv_cuda(A, d_x, d_Ap);
waxpby_cuda(n, 1.0, d_b, -1.0, d_Ap, d_r);
ddot_cuda(n, d_r, d_r, d_partial, d_rtrans);

for (int k = 1; k < max_iter; ++k) {

    // update p on the GPU
    ...

    spmv_cuda(A, d_p, d_Ap);

    ddot_cuda(n, d_p, d_Ap, d_partial, d_pAp);

    // compute alpha and update x and r on the GPU
    ...
}

// Copy final solution back.
cudaMemcpy(x, d_x, ...);

"GPU-resident" here does not mean that the entire CG solver has been turned into one giant CUDA kernel.

The host still orchestrates a sequence of kernels. What changed is that the large working vectors no longer travel between CPU and GPU after every operation.

Commits a4c1be8 — Fully GPU resident CG kernel
bfedb51 — Add GPU-resident native CUDA HPCCG baseline

The result depends strongly on problem size

This produced an interesting result.

For the small 20^3 problem:

The GPU implementation is about 3× slower.

But for 100^3:

Now the GPU implementation is approximately 7.2× faster end-to-end.

This is one of the most useful results of the exercise.

Both problems require the same 149 CG iterations. But 20^3 contains only 8,000 rows, whereas 100^3 contains one million — 125 times as many rows.

For the small problem, there is simply not enough useful work to amortize repeated kernel launches, reductions, and synchronization.

For the larger problem, GPU throughput dominates those fixed costs.

So rather than saying simply that "the GPU is faster," I find it more useful to think in terms of an application-level crossover:

small problem
    useful GPU work < GPU overhead
              ↓
          CPU wins

large problem
    useful GPU work >> GPU overhead
              ↓
          GPU wins

A profiling caveat: CPU timers are not GPU timers

There is another subtle lesson hidden in these results.

HPCCG's original timing infrastructure was written for synchronous CPU functions. When it surrounds a CUDA kernel launch, however, the host can return before the GPU has actually finished executing the kernel.

For example, the GPU-resident 100^3 run reports:

Total:     0.468121 s
DDOT:      0.002686 s
WAXPBY:    0.002234 s
SPARSEMV:  0.002139 s

Those per-operation numbers should not be interpreted as accumulated GPU execution times. They largely reflect asynchronous host-side launch behavior.

I therefore use the application's total wall-clock time for the end-to-end comparison, and Nsight Systems or Nsight Compute when reasoning about individual CUDA operations.

Porting an application to an asynchronous execution model can change not only its performance, but also the meaning of its existing profiler.

6. Following the Bottleneck

Once the large vector transfers were gone, the smaller synchronization points became visible.

Eliminating an unnecessary scalar copy

One particularly simple example was the residual scalar.

Initially, before computing the next residual, I preserved the previous value using a device-to-device copy.

But there was no reason to move the data at all.

I already had two scalar buffers. They could simply exchange roles:

std::swap(d_rtrans, d_oldrtrans);

ddot_cuda(n, d_r, d_r, d_partial, d_rtrans);

compute_division<<<1,1>>>(
    d_rtrans,
    d_oldrtrans,
    d_beta
);

The swap changes two host-side pointer values. No GPU data moves.

That one change reduced the number of cudaMemcpy calls in the profile from:

304 → 156

Exactly 148 copies disappeared.

The interesting part is what happened next.

The total API time attributed to cudaMemcpy barely changed, and end-to-end runtime improved only modestly, from roughly 0.55 s to around 0.49 s in those runs.

Why?

Because the removed copies were cheap device-to-device scalar copies. The expensive calls that remained were primarily host-visible synchronization points.

This is another useful profiling lesson:

Removing many operations is not necessarily the same thing as removing the expensive operations.

Inside the SpMV kernel

Once application dataflow was no longer dominating the discussion, I returned to the original hotspot: SpMV.

Nsight Compute showed the following picture for the 100^3 one-thread-per-row CSR kernel:

At first glance, the occupancy looks excellent.

But occupancy only tells us how many warps can reside on the SM. It does not tell us how many of those warps are actually ready to issue an instruction.

Here, roughly 11 warps were active per scheduler, while only about 0.06 warps were eligible on average. The scheduler had no eligible warp during approximately 98% of sampled cycles.

The source counters provided another clue: Nsight Compute reported approximately 43.7 million excessive memory sectors, about 72% of the total sectors, associated with uncoalesced global accesses.

So increasing occupancy would be the wrong first optimization target. There are already plenty of resident warps.

The more interesting problem is memory access efficiency.

The indirect access

x[A.cols[j]]

is an obvious candidate, but it is not necessarily the only culprit. In a one-thread-per-row CSR mapping, neighboring threads also walk different segments of vals and cols.

Nsight Compute tells us that the current memory access pattern is inefficient. Determining exactly which loads dominate would be the next controlled experiment.

And that is where I intentionally stop optimizing this CUDA implementation for now.

The point of this exercise was not to squeeze every last percent out of CSR SpMV. We now have enough information to ask a different question:

What happens if we solve the same porting problem using a performance-portability abstraction instead of writing CUDA directly?