将共轭梯度求解器移植到 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: profile the original CPU application and identify the dominant kernel, move sparse matrix-vector multiplication (SpMV) to the GPU, migrate the remaining vector operations and reductions, examine why accelerating individual kernels is not enough, keep the CG working set resident on the GPU, and 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 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 } Enter fullscreen mode Exit fullscreen mode 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: KernelTime (s)Fraction of total SpMV 2.819 83.5% WAXPBY 0.342 10.1% DDOT 0.215 6.4% Total 3.378 100% 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 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 > GPU overhead ↓ GPU wins Enter fullscreen mode Exit fullscreen mode 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 Enter fullscreen mode Exit fullscreen mode 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>>( d_rtrans, d_oldrtrans, d_beta ); Enter fullscreen mode Exit fullscreen mode 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 Enter fullscreen mode Exit fullscreen mode 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: MetricMeasurement Achieved occupancy 90.9% Theoretical occupancy 100% DRAM throughput 63.3% of peak SM compute throughput 11.7% Active warps / scheduler ~11.1 Eligible warps / scheduler ~0.06 L1 hit rate 54.0% L2 hit rate 63.1% 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]] Enter fullscreen mode Exit fullscreen mode 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?