M5 GPUs and GPGPU

GPUs and GPGPU A GPU (Graphics Processing Unit) is a specialized electronic circuit designed to rapidly manipulate and alter memory to accelerate the creation of images in a frame buffer intended for output to a display device. Their highly parallel structure makes them more efficient than general-purpose CPUs for algorithms where large blocks of data are processed in parallel. GPGPU (General-Purpose computing on Graphics Processing Units) is the use of a GPU, which would typically handle computer graphics, to perform computation in applications traditionally handled by the CPU. This is possible because modern GPUs have programmable shaders and higher-precision arithmetic, allowing programmers to execute algorithms that are not related to graphics. The programming model used to achieve this is often CUDA (Compute Unified Device Architecture) for NVIDIA GPUs. The main idea is to offload computationally intensive parallel portions of an application to the GPU, while the rest of the code (which is often sequential) runs on the CPU. GPU Architectures Modern NVIDIA GPUs are built around a scalable array of Streaming Multiprocessors (SMs). Each GPU has several SMs, and each SM is designed to execute hundreds of threads concurrently. The key components of a typical GPU architecture are: • Streaming Multiprocessor (SM): The core processing unit of the GPU. A GPU contains multiple SMs. • CUDA Cores: Each SM contains numerous processing cores (often called CUDA cores) that execute the actual instructions. • Memory Hierarchy: GPUs have a complex memory hierarchy to manage data for the thousands of threads running in parallel. o Registers: The fastest memory, private to each thread. o Shared Memory: A small, low-latency memory space shared by all threads within a single block. It's much faster than global memory. o L1/L2 Cache: Caches shared by cores within an SM (L1) or across all SMs (L2) to speed up access to global memory. o Global Memory (DRAM): The largest memory space on the GPU, accessible by all threads from all blocks. It has high latency but high bandwidth. Diagram: Simplified GPU Architecture +-----------------------------------------------------------------+ | GPU | | +-----------------+ +-----------------+ +-----------------+ | || SM 1 | | SM 2 | | ... | | | | +-----------+ | | +-----------+ | | +-----------+ | | | | | Core | Core | | | | Core | Core | | | | Core | Core | | | | | | ... | ... | | | | ... | ... | | | | ... | ... | | | | | +-----------+ | | +-----------+ | | +-----------+ | | | | +-----------+ | | +-----------+ | | +-----------+ | | | | | Shared Mem| | | | Shared Mem| | | | Shared Mem| | | | | +-----------+ | | +-----------+ | | +-----------+ | | | +-----------------+ +-----------------+ +-----------------+ | | | | +----------------------------+ | | | +----------------------------+ | L2 Cache / Global Memory (DRAM) | | +-----------------------------------------------------------------+ | Heterogeneous Computing Heterogeneous computing refers to systems that use more than one kind of processor or core. In the context of CUDA, it means using both the CPU (the host) and the GPU (the device) together to solve a problem. • Host (CPU): Manages the overall program flow, handles I/O operations, and executes sequential parts of the code. • Device (GPU): Executes the highly parallel, computationally intensive parts of the code, called kernels. A typical CUDA program follows this pattern: 1. The CPU initializes data. 2. The CPU allocates memory on the GPU and copies the data from host memory to GPU memory. 3. The CPU launches the kernel on the GPU, which runs across thousands of threads. 4. The GPU threads execute the kernel in parallel. 5. After the kernel completes, the CPU copies the results back from GPU memory to host memory. 6. The CPU frees the allocated GPU memory. Diagram: Host-Device Interaction +-----------------+ +------------------+ | Host (CPU) | | Device (GPU) | |-----------------| |------------------| | 1. Allocate host| | | memory (A,B) | | | 2. Initialize | | | 3. Allocate dev | | | | data in A,B | | | | | | | | | 3. dev_A, dev_B, | | memory | | dev_C created | | (dev_A,B,C) | | | | | | | | 4. Copy A,B to |------ HtoD Data Transfer ------>| 4. dev_A, dev_B | | dev_A, dev_B | | | | populated | | | | 5. Launch Kernel|-------- Kernel Launch -------->| 5. Kernel runs | | `add(dev_A,..)`| | on threads | | | | C=A+B | | | | | | 6. Copy dev_C to| <----- DtoH Data Transfer ------| 6. Result in | host memory C| | | | dev_C | | | | 7. Free dev mem | | | | 8. Use results C| | | +-----------------+ +------------------+ | Threads, Blocks, and Grids CUDA organizes its parallel threads into a three-level hierarchy: threads, blocks, and a grid. • Thread: The smallest unit of execution. It runs an instance of the kernel function. • Block: A group of threads. Threads within a block can cooperate by sharing data through shared memory and synchronizing their execution using __syncthreads(). • Grid: A collection of thread blocks. When you launch a kernel, you specify the number of blocks in the grid and the number of threads in each block. CUDA provides built-in variables that threads can use to identify their unique position: • threadIdx.x: The thread's index within its block. • blockIdx.x: The block's index within the grid. • blockDim.x: The number of threads in a block. • gridDim.x: The number of blocks in the grid. These variables can be 1D, 2D, or 3D (e.g., threadIdx.y, blockIdx.z). Diagram: Grid, Blocks, and Threads +-------------------------------------------------+ | Grid (A collection of all blocks) | | | | +----------+ +----------+ +----------+ | | | Block 0 | | Block 1 | | Block 2 | ... | | | T0 T1 T2 | | T0 T1 T2 | | T0 T1 T2 | | | | T3 T4 T5 | | T3 T4 T5 | | T3 T4 T5 | | | +----------+ +----------+ +----------+ | | | +-------------------------------------------------+ A thread's global, unique index can often be calculated as: int idx = blockIdx.x * blockDim.x + threadIdx.x; Nvidia Compute Capabilities Compute Capability (CC) describes the hardware features of a GPU. It's represented by a version number (e.g., 5.2, 6.1, 7.5, 8.6). A higher compute capability generally means more features, better performance, and more cores/memory. Key features determined by CC include: • Maximum number of threads per block. • Amount of shared memory and registers per SM. • Support for specific arithmetic operations (e.g., double-precision floatingpoint, tensor core operations). • Concurrent kernel execution capabilities. When compiling CUDA code, you can target a specific compute capability to ensure compatibility and leverage hardware-specific optimizations. Vector Addition This is the "Hello, World!" of CUDA programming. The goal is to add two vectors, A and B, and store the result in vector C. Each thread is responsible for computing one element of the result vector C. Returning Results from CUDA Kernels Notice that the kernel function add_vectors has a void return type. CUDA kernels do not return values directly. Instead, they write their results to pointers that point to device (GPU) memory. The host code must then explicitly copy this data back from the device to the host to access the results. Program: Vector Addition in CUDA C++ #include #include // Kernel function that executes on the GPU __global__ void add_vectors(int n, const float* a, const float* b, float* c) { // Calculate the global index of the thread int idx = blockIdx.x * blockDim.x + threadIdx.x; // Each thread adds one pair of elements, if the index is within bounds if (idx < n) { c[idx] = a[idx] + b[idx]; } } int main() { int n = 10000; size_t size = n * sizeof(float); // 1. Allocate host memory float* h_a = (float*)malloc(size); float* h_b = (float*)malloc(size); float* h_c = (float*)malloc(size); // Initialize host vectors for (int i = 0; i < n; ++i) { h_a[i] = static_cast(i); h_b[i] = static_cast(i * 2); } // 2. Allocate device memory float *d_a, *d_b, *d_c; cudaMalloc(&d_a, size); cudaMalloc(&d_b, size); cudaMalloc(&d_c, size); // 3. Copy data from host to device cudaMemcpy(d_a, h_a, size, cudaMemcpyHostToDevice); cudaMemcpy(d_b, h_b, size, cudaMemcpyHostToDevice); // 4. Launch the kernel int threadsPerBlock = 256; // Calculate the number of blocks needed int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock; add_vectors<<>>(n, d_a, d_b, d_c); // 5. Copy result from device to host cudaMemcpy(h_c, d_c, size, cudaMemcpyDeviceToHost); // Verify a result (optional) std::cout << "Result verification: h_c[100] = " << h_c[100] << " (Expected: 300)" << std::endl; // 6. Free memory free(h_a); free(h_b); free(h_c); cudaFree(d_a); cudaFree(d_b); cudaFree(d_c); return 0; } CUDA Trapezoidal Rule I: A Naive Approach The trapezoidal rule approximates a definite integral int_abf(x)dx. The formula is: ∫abf(x)dx≈2h[f(x0)+f(xn)+2i=1∑n−1f(xi)] where h=(b−a)/n and x_i=a+ih. A simple parallel approach is to have each thread calculate the area of one small trapezoid and add its area to a single global sum. This requires an atomic operation to prevent race conditions where multiple threads try to update the sum simultaneously. Program: Kernel for Naive Trapezoidal Rule C++ // A simple function to integrate, e.g., f(x) = x*x __device__ float f(float x) { return x * x; } // Naive Kernel: Each thread updates a single global sum __global__ void trap_naive(float a, float b, int n, float* global_sum) { float h = (b - a) / n; int i = blockIdx.x * blockDim.x + threadIdx.x; if (i < n) { float x_i = a + i * h; float x_j = a + (i + 1) * h; float trapezoid_area = (f(x_i) + f(x_j)) / 2.0f * h; // This is a bottleneck! All threads contend for this single operation. atomicAdd(global_sum, trapezoid_area); } } Problem: The atomicAdd operation serializes memory access. As thousands of threads try to add to global_sum, they must do so one at a time, creating a massive performance bottleneck. CUDA Trapezoidal Rule II: Improving Performance To reduce the atomicAdd bottleneck, we can use a reduction strategy. 1. Each thread still calculates the area of one trapezoid. 2. However, instead of adding to a single global sum, each block calculates its own partial sum. 3. Each thread in a block adds its result to a variable that is local to the block (e.g., in shared memory). At the end, one thread per block writes the block's subtotal to a global array. 4. The host (CPU) then sums up the relatively small number of partial sums from the blocks. This drastically reduces the number of high-latency global atomic operations from n to num_blocks. Program: Kernel with Per-Block Sums C++ // Kernel with per-block reduction __global__ void trap_per_block(float a, float b, int n, float* block_sums) { // Use shared memory to store the partial sum for this block __shared__ float local_sum; if (threadIdx.x == 0) { local_sum = 0.0f; } __syncthreads(); // Ensure local_sum is initialized float h = (b - a) / n; int i = blockIdx.x * blockDim.x + threadIdx.x; if (i < n) { float x_i = a + i * h; float x_j = a + (i + 1) * h; float trapezoid_area = (f(x_i) + f(x_j)) / 2.0f * h; atomicAdd(&local_sum, trapezoid_area); } __syncthreads(); // Ensure all threads have finished adding // Thread 0 of each block writes the block's sum to the global array if (threadIdx.x == 0) { block_sums[blockIdx.x] = local_sum; } } Improvement: This is much faster because atomic operations are now on fast shared memory (local_sum). The only global memory writes are one per block. The final, small summation is done on the CPU. CUDA Trapezoidal Rule III: Blocks with More Than One Warp The previous version still uses atomicAdd within the block's shared memory. We can eliminate atomics entirely by performing an efficient, tree-based reduction in shared memory. A warp is a group of 32 threads that execute the same instruction at the same time. The strategy is to leverage this synchronous execution. The Algorithm: 1. Each thread calculates its local value and stores it in a shared memory array. 2. The threads perform a parallel reduction on that array. In successive steps, the number of active threads is halved. For example, thread t adds the value from thread t + stride to its own element. The stride starts at blockDim.x / 2 and is halved in each step. 3. __syncthreads() is crucial between each step of the reduction to ensure writes from one step are visible to the next. Diagram: Shared Memory Tree Reduction Let's say we have 8 threads in a block: Shared Memory Array: [v0, v1, v2, v3, v4, v5, v6, v7] Step 1 (stride = 4): - T0 adds v4 -> v0 = v0+v4 - T1 adds v5 -> v1 = v1+v5 - T2 adds v6 -> v2 = v2+v6 - T3 adds v7 -> v3 = v3+v7 Array state: [v0+v4, v1+v5, v2+v6, v3+v7, -, -, -, -] __syncthreads() Step 2 (stride = 2): - T0 adds v2 -> v0 = v0+v2 - T1 adds v3 -> v1 = v1+v3 Array state: [v0+v4+v2+v6, v1+v5+v3+v7, -, -, -, -, -, -] __syncthreads() Step 3 (stride = 1): - T0 adds v1 -> v0 = v0+v1 Array state: [Final Sum, -, -, -, -, -, -, -] __syncthreads() Thread 0 now holds the total sum for the block. Program: Kernel with Shared Memory Reduction C++ __global__ void trap_optimized(float a, float b, int n, float* block_sums) { // Shared memory for the reduction __shared__ float temp[1024]; // Max threads per block int tid = threadIdx.x; int global_idx = blockIdx.x * blockDim.x + tid; float h = (b - a) / n; // 1. Each thread calculates its value and stores it in shared memory if (global_idx < n) { float x_i = a + global_idx * h; float x_j = a + (global_idx + 1) * h; temp[tid] = (f(x_i) + f(x_j)) / 2.0f * h; } else { temp[tid] = 0.0f; } __syncthreads(); // 2. Perform the tree-based reduction for (int stride = blockDim.x / 2; stride > 0; stride /= 2) { if (tid < stride) { temp[tid] += temp[tid + stride]; } __syncthreads(); // Crucial synchronization! } // 3. Thread 0 writes the final block sum to global memory if (tid == 0) { block_sums[blockIdx.x] = temp[0]; } } This is a highly efficient parallel reduction pattern in CUDA and is much faster than the previous versions because it avoids atomic operations and leverages the fast shared memory and synchronous execution of warps.