Understanding CUDA Programming

When leveraging GPU processing to solve computational problems, we must first transfer data from CPU memory to GPU memory via high-speed interconnects like PCIe or NVLink. The second step involves executing the actual computation, typically through CUDA kernel functions. Finally, after computation completes, results are transferred back from GPU to CPU memory, potentially using network or disk storage for persistence.

Basic CUDA Syntax

CUDA introduces special function qualifiers to distinguish code that runs on the GPU. The __global__ qualifier indicates a function that executes on the GPU and can be called from the host (CPU):

__global__ void compute_kernel(void) {
    // GPU computation code here
}

The __global__ qualifier signals to the CUDA compiler (NVCC) that this function should be compiled for GPU execution. NVCC orchestrates the compilation process, potentially invoking multiple compilers and tools to transform the code into separate host and device components. Host code typical compiles with standard GCC/G++, while __global__ functions compile to GPU-executable code.

Kernel Execution

To invoke a kernel function, we use special syntax that specifies execution configuration parameters:

compute_kernel<<<grid_size block_size="">>>( );
</grid_size>

The first parameter (grid_size) indicates how many thread blocks to launch, while the second (block_size) specifies the number of threads per block. For example:

compute_kernel<<<1,1>>>( );  // Launch 1 block with 1 thread

Memory Management in CUDA

CUDA provides specialized functions for memory allocation and data transfer between host and device:

• cudaMalloc() - Allocate memory on the GPU
• cudaFree() - Free GPU memory
• cudaMemcpy() - Copy data between host and device

These APIs use pointers to reference memory spaces. It's crucial to remember that CPU pointers should not be dereferenced in device code, and vice versa, as they refer to distinct memory spaces managed by different processors.

Thread Hierarchy in CUDA

CUDA organizes threads into a three-level hierarchy: grids contain blocks, and blocks contain threads. This hierarchy enables massive parallelism:

• Grid: Collection of blocks that can execute in parallel
• Block: Group of threads that can cooperate and share memory
• Thread: Basic execution unit with its own private memory

Blocks are identified by blockIdx, and threads within blocks by threadIdx. For example, in a 1D configuration:

__global__ void vector_add(float *input_a, float *input_b, float *result) {
    int idx = blockIdx.x;  // Current block index
    result[idx] = input_a[idx] + input_b[idx];
}

Host Code Implementation

A complete CUDA program involves both host and device code. Here's a basic template:

#include <cuda_runtime.h>

__global__ void kernel_name(float *input_a, float *input_b, float *result, int n) {
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    if (index < n)
        result[index] = input_a[index] + input_b[index];
}

int main() {
    // Set GPU device
    cudaSetDevice(0);
    
    // Allocate host and device memory
    float *h_a, *h_b, *h_result;
    float *d_a, *d_b, *d_result;
    int size = N * sizeof(float);
    
    h_a = (float *)malloc(size);
    h_b = (float *)malloc(size);
    h_result = (float *)malloc(size);
    
    cudaMalloc((void **)&d_a, size);
    cudaMalloc((void **)&d_b, size);
    cudaMalloc((void **)&d_result, size);
    
    // Initialize host data
    initialize_arrays(h_a, h_b, N);
    
    // Copy data from host to device
    cudaMemcpy(d_a, h_a, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, size, cudaMemcpyHostToDevice);
    
    // Execute kernel on device
    int threadsPerBlock = 256;
    int blocks = (N + threadsPerBlock - 1) / threadsPerBlock;
    kernel_name<<<blocks threadsperblock="">>>(d_a, d_b, d_result, N);
    
    // Copy result from device to host
    cudaMemcpy(h_result, d_result, size, cudaMemcpyDeviceToHost);
    
    // Release memory
    free(h_a);
    free(h_b);
    free(h_result);
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_result);
    
    return 0;
}
</blocks></cuda_runtime.h>

Handling Index Calculations

When working with CUDA, proper index calculation is essential. For 1D grids and blocks:

int global_index = threadIdx.x + blockIdx.x * blockDim.x;

For 2D configurations:

int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
int global_index = row * grid_width + col;

Handling Edge Cases

When the number of elements isn't evenly divisible by the number of threads, we need to handle edge cases:

__global__ void safe_kernel(float *input_a, float *input_b, float *result, int n) {
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    if (index < n)  // Check bounds to prevent overflow
        result[index] = input_a[index] + input_b[index];
}

When launching the kernel, we calculate the number of blocks needed:

int threads_per_block = 256;
int blocks = (N + threads_per_block - 1) / threads_per_block;
safe_kernel<<<blocks threads_per_block="">>>(d_a, d_b, d_result, N);
</blocks>

Matrix Multiplication Example

Matrix multiplication is a common operation that benefits from GPU acceleration. First, let's examine the CPU implementation:

// CPU matrix multiplication
// A: h x k matrix, B: k x w matrix, C: h x w matrix (result)
void cpu_matmul(float *A, float *B, float *C, int h, int k, int w) {
    for (int row = 0; row < h; row++) {
        for (int col = 0; col < w; col++) {
            float sum = 0.0f;
            for (int i = 0; i < k; i++) {
                sum += A[row * k + i] * B[i * w + col];
            }
            C[row * w + col] = sum;
        }
    }
}

Now, let's implement the same operation on GPU:

__global__ void gpu_matmul(float *matrix_a, float *matrix_b, float *result_matrix, 
                          int height, int width, int inner_dim) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    
    // Check if within matrix bounds
    if (row >= height || col >= width) return;
    
    float sum = 0.0f;
    for (int i = 0; i < inner_dim; i++) {
        sum += matrix_a[row * inner_dim + i] * matrix_b[i * width + col];
    }
    result_matrix[row * width + col] = sum;
}

To launch this kernel:

// Launch parameters
int threads_per_block = 16;  // 16x16 threads per block
dim3 blocks((width + threads_per_block - 1) / threads_per_block,
            (height + threads_per_block - 1) / threads_per_block);

// Execute kernel
gpu_matmul<<<blocks dim3="" threads_per_block="">>>(
    d_matrix_a, d_matrix_b, d_result_matrix, height, width, inner_dim);
</blocks>