T02: GPU, parallel, and graphics computing

A. Prerequisites and observable goals#

Level: Advanced-entry track. CPU lab tested; optional CUDA source is supplied but not compiled or executed in this environment.

B. The central question: is there enough independent work to pay for coordination?#

A GPU can execute many arithmetic operations across many elements, but offloading requires more than replacing a loop with a kernel. Inputs must be accessible in the appropriate memory, work must be decomposed into independent or synchronized pieces, and outputs must be validated. Small jobs can spend more time launching and transferring than computing. A CPU implementation remains both a fallback and a correctness oracle.

SIMD applies one instruction to multiple data lanes. GPU SIMT execution organizes groups of threads that execute instructions together while exposing thread-level indexing. Divergent control flow can cause different paths to execute under masks rather than giving every branch a separate free processor. A mathematical operation’s parallelism does not establish its achieved hardware utilization.

In the basic CUDA model used by the optional example, a launch creates a grid of blocks, and each block contains threads. Threads in a block can cooperate through shared memory and block-scoped synchronization. Blocks ordinarily must be independently schedulable; do not implement a global barrier by having arbitrary blocks spin for one another. Specialized features exist on particular hardware, but they are outside this introductory kernel. Current NVIDIA documentation was checked for this model and the standard synchronization concepts, not by executing the kernel on a GPU. [CUDA]

C. Memory hierarchy and synchronization scopes#

Registers hold thread-local working values; shared memory provides explicitly managed block-local storage; global device memory stores larger arrays. These names specify access and lifetime concepts, not a promise that every access has a fixed latency. Coalesced accesses let neighboring lanes access nearby global locations efficiently. Random gathers, redundant loads, and low reuse can make data movement dominate arithmetic.

A barrier requires the relevant participating threads to reach it under valid control flow. Returning early from some threads before a block barrier can invalidate the intended algorithm. Instead, out-of-range threads contribute the operation’s identity, such as zero for addition, and participate in the reduction. Host code also needs an appropriate synchronization or event mechanism before consuming a kernel result; a launch can be asynchronous with respect to the host. [CUDA]

D. Derive a tree reduction#

To sum eight numbers, first sum adjacent pairs in parallel, then pairs of those results, then the final two values. The total number of additions is $7$, the same as a sequential sum, while the dependency depth is $3=\log_2 8$. For $n$ a power of two, work is $n-1$ additions and span is $\log_2 n$. Ideal execution on $p$ processors is bounded below by both $W/p$ and $S$, so $T_p\ge\max(W/p,S)$ in an idealized unit-cost model.

On a GPU, each block can reduce a chunk and write one partial sum. A second kernel or a CPU can combine partials. The supplied optional kernel loads up to two input values per thread, places local sums in shared memory, and repeatedly halves the active reduction width with barriers. The CPU combines the much smaller partial array in double precision. This is a teaching implementation, not an optimized production reduction library.

Floating-point addition is not associative. The reduction tree changes rounding order, so exact equality to a serial floating-point sum is not a universal correctness requirement. Compare against a higher-quality reference such as math.fsum for the CPU lab, use absolute and relative tolerances appropriate to magnitude, and test cancellation and dynamic range. A fast but numerically unsuitable reduction is not interchangeable with a stable one.

E. Measure end-to-end and kernel-only costs separately#

A simple offload model is

If CPU time is $2$ ms, launch plus transfer is $3$ ms, and kernel work is $0.2$ ms, end-to-end GPU time is $3.2$ ms: slower despite a ten-times-faster compute phase. Keeping data on the device across several operations can amortize transfer, but then the compared workload must include those operations and their memory use.

Arithmetic intensity is useful operations per byte moved under a stated memory boundary. A reduction has little arithmetic per input byte and is often sensitive to bandwidth and reduction overhead. A matrix multiply can reuse tiles of inputs across many multiply-adds, raising arithmetic intensity when implemented well. Counting theoretical reuse is not proof that an implementation achieves it.

The optional CUDA program reports device-event kernel time and a host-wall interval including its transfer/launch/synchronization work, with allocation policy documented in the source. No CUDA result is included in this book because nvcc and a GPU were not available where its code was tested. The CPU reduction tests did run. This separation prevents a hardware capability assumption from becoming a fabricated benchmark.

F. Graphics is mathematics plus a staged computation#

A simple graphics pipeline transforms object coordinates into world, view, and clip coordinates, assembles primitives, rasterizes covered samples, and evaluates shading and depth tests. Homogeneous coordinates allow translation to be represented by a matrix acting on an augmented vector $(x,y,z,1)$. Perspective division converts clip coordinates to normalized coordinates by dividing by the appropriate homogeneous component; points with problematic denominators require clipping rather than blind division.

Normals do not generally transform like position vectors under nonuniform scaling. To preserve orthogonality to transformed tangents, the normal direction transforms by the inverse transpose of the invertible linear part, then is normalized. This follows from $(A^{-T}n)^T(At)=n^Tt=0$. Singular transforms require separate handling. This is an example of linear algebra determining implementation correctness rather than serving as decorative notation.

Parallel graphics work still needs consistent resource lifetimes, synchronization, color and coordinate conventions, and performance measurement. Rendering a plausible image is not proof that a transform or physical model is correct. Use small geometry with analytically known positions and normals before testing a complex scene.

G. Independent practice#

F foundation C standard course work T transfer S optional stretch L runnable lab or project. Attempt the complete set before opening its solutions.

H. Complete solutions: T02#

Compare the first step where your work differs. A corrected answer is not yet an independent success.

Open every solution

"""T02.05 CPU fallback: changing reduction order changes rounding."""
import math

def pairwise(values: list[float]) -> float:
    work = list(values)
    if any(not math.isfinite(x) for x in work):
        raise ValueError("finite inputs required")
    while len(work) > 1:
        work = [work[i] + work[i + 1] if i + 1 < len(work) else work[i]
                for i in range(0, len(work), 2)]
    return work[0] if work else 0.0

def kahan(values: list[float]) -> float:
    total = correction = 0.0
    for value in values:
        if not math.isfinite(value):
            raise ValueError("finite inputs required")
        y = value - correction
        new = total + y
        correction = (new - total) - y
        total = new
    return total

if __name__ == "__main__":
    values = [0.1] * 10000
    print("loop", sum(values), "tree", pairwise(values),
          "kahan", kahan(values), "fsum", math.fsum(values))

Optional hardware extension: supplied but untested CUDA source#

The following complete program is included so a future GPU run has a concrete starting point. Compile only on a configured CUDA system, validate its result before interpreting timings, and record toolkit, compiler, driver, device, workload, and transfer policy. No speedup or successful compilation is claimed for this source here.

// T02 optional CUDA extension. Not compiled or executed in the reference environment.
#include <cuda_runtime.h>
#include <chrono>
#include <cmath>
#include <iostream>
#include <numeric>
#include <stdexcept>
#include <vector>
void check(cudaError_t status) {
    if(status!=cudaSuccess) throw std::runtime_error(cudaGetErrorString(status));
}
struct DeviceArray {
    float* data=nullptr;
    explicit DeviceArray(std::size_t n) {
        check(cudaMalloc(reinterpret_cast<void**>(&data),n*sizeof(float)));
    }
    ~DeviceArray() { if(data) cudaFree(data); }
    DeviceArray(const DeviceArray&)=delete;
    DeviceArray& operator=(const DeviceArray&)=delete;
};
struct Event {
    cudaEvent_t handle;
    Event() { check(cudaEventCreate(&handle)); }
    ~Event() { cudaEventDestroy(handle); }
    Event(const Event&)=delete;
    Event& operator=(const Event&)=delete;
};
constexpr unsigned block_size=256;
__global__ void reduce_blocks(const float* input,float* output,std::size_t n) {
    __shared__ float partial[block_size];
    const unsigned lane=threadIdx.x;
    const std::size_t i=std::size_t(blockIdx.x)*block_size*2+lane;
    float value=i<n?input[i]:0.0F;
    if(i+block_size<n) value+=input[i+block_size];
    partial[lane]=value;
    __syncthreads();
    for(unsigned stride=block_size/2;stride;stride/=2) {
        if(lane<stride) partial[lane]+=partial[lane+stride];
        __syncthreads(); // all threads participate, including inactive lanes
    }
    if(lane==0) output[blockIdx.x]=partial[0];
}
int main() {
    try {
        constexpr std::size_t n=1U<<20;
        const unsigned blocks=static_cast<unsigned>((n+2*block_size-1)/(2*block_size));
        std::vector<float> input(n,1.0F),partial(blocks);
        const auto begin=std::chrono::steady_clock::now();
        DeviceArray device_input(n),device_output(blocks);
        Event start,stop;
        check(cudaMemcpy(device_input.data,input.data(),n*sizeof(float),cudaMemcpyHostToDevice));
        check(cudaEventRecord(start.handle));
        reduce_blocks<<<blocks,block_size>>>(device_input.data,device_output.data,n);
        check(cudaGetLastError());
        check(cudaEventRecord(stop.handle));
        check(cudaEventSynchronize(stop.handle));
        float kernel_ms=0;
        check(cudaEventElapsedTime(&kernel_ms,start.handle,stop.handle));
        check(cudaMemcpy(partial.data(),device_output.data,blocks*sizeof(float),cudaMemcpyDeviceToHost));
        const double result=std::accumulate(partial.begin(),partial.end(),0.0);
        const double total_ms=std::chrono::duration<double,std::milli>(
            std::chrono::steady_clock::now()-begin).count();
        if(result!=static_cast<double>(n)) throw std::runtime_error("wrong reduction");
        std::cout<<"n="<<n<<" result="<<result<<" kernel_ms="<<kernel_ms
                 <<" allocation_transfer_compute_ms="<<total_ms<<"\n";
    } catch(const std::exception& e) { std::cerr<<e.what()<<"\n"; return 1; }
}

I. Retention and readiness#

Retrieve A01 transformations, A05 rounding, C10 synchronization scopes, and T01 timing boundaries. Readiness requires a correct reduction and a defensible decision about whether offload pays for its costs. Continue with verified hardware tooling, parallel algorithms, graphics pipeline mathematics, and numerical error analysis; profile before replacing a maintained GPU library with a custom kernel.