Launching CUDA Kernels for Element-wise Operations

1. Writing a CUDA Kernel

A kernel is a C++ function that is executed by many GPU threads in parallel. To distinguish it from regular C++ functions that run on the host (CPU), we use a special CUDA keyword.

The __global__ Keyword

You declare a function as a CUDA kernel by using the __global__ specifier before its return type. This tells the CUDA compiler (nvcc) two things:

1. This function will be executed on the device (GPU).
2. It can be called from the host (CPU).

Here is the skeleton of our vector addition kernel:

__global__ void vectorAdd(float *A, float *B, float *C, int N) {
    // Kernel code goes here
}

Notice that the parameters are pointers to device memory (which you allocated with cudaMalloc) and a simple integer passed by value.

Unique Identification: The Core of Parallelism

Imagine you've hired thousands of workers to assemble a product, where each product has a unique serial number. To work efficiently, the first thing each worker needs to know is which serial number they are responsible for.

It's the same for GPU threads. When we launch our kernel, thousands of threads will execute the exact same code. To ensure each thread operates on a different element of our vectors, it must calculate its own unique, global ID.

CUDA provides built-in variables for this purpose:

• threadIdx: A dim3 variable containing the index of the thread within its block.
• blockIdx: A dim3 variable containing the index of the block within the grid.
• blockDim: A dim3 variable containing the dimensions of the block (i.e., the number of threads per block).

For a 1D problem like vector addition, we can use the .x component of these variables to compute a unique global index for each thread:

This formula is fundamental to CUDA programming. It works like this: blockIdx.x * blockDim.x calculates the starting index for the current block, and threadIdx.x provides the offset within that block.

For a more detailed visual explanation of this critical concept, please watch the following video.

The vectorAdd Kernel, Complete

Now we can write the full kernel. Each thread will calculate its unique index i, check if that index is within the bounds of our arrays, and if so, perform the addition for that element.

__global__ void vectorAdd(const float *A, const float *B, float *C, int numElements) {
    // Calculate the unique global thread index
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    // Boundary check: Make sure we don't access memory out of bounds
    if (i < numElements) {
        C[i] = A[i] + B[i];
    }
}

The boundary check if (i < numElements) is essential. We often launch slightly more threads than elements for convenience (e.g., rounding up to a full block). This if statement ensures that those extra threads do nothing, preventing them from accessing invalid memory, which would cause a crash or undefined behavior. Your C++ experience with array bounds will make the importance of this check immediately clear.

2. Launching a Kernel

With the kernel written, we now need to launch it from our host code. This is done by calling the kernel function, but with a special execution configuration syntax: the triple angle brackets <<<...>>>.

kernel_name<<<blocksPerGrid, threadsPerBlock>>>(argument1, argument2, ...);

You need to specify two parameters:

1. blocksPerGrid: The number of thread blocks to launch in the grid.
2. threadsPerBlock: The number of threads to launch in each block.

The total number of threads launched will be blocksPerGrid * threadsPerBlock.

How do we choose these values?

• threadsPerBlock: A good starting point is a multiple of 32 (the size of a "warp," the fundamental unit of execution scheduling on NVIDIA GPUs). Common choices are 128, 256, or 512. There is a hardware limit, typically 1024 threads per block. Let's pick 256 for our example.
• blocksPerGrid: This is calculated to ensure we have enough threads to cover all numElements. A standard way to do this is with integer arithmetic that performs a ceiling division:
  (numElements + threadsPerBlock - 1) / threadsPerBlock

So, the host code to launch our vectorAdd kernel would look like this:

// Define the vector size
int numElements = 50000;

// Set the number of threads per block
int threadsPerBlock = 256;

// Calculate the number of blocks in the grid
int blocksPerGrid = (numElements + threadsPerBlock - 1) / threadsPerBlock;

// Launch the kernel!
vectorAdd<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, numElements);

With numElements = 50000 and threadsPerBlock = 256, this would launch (50000 + 255) / 256 = 196 blocks. The total number of threads is 196 * 256 = 50176, which is why the boundary check inside the kernel is so important.

3. The Complete CUDA Program

We've now covered all the individual steps. Let's put them together to see the full picture of a CUDA program, from memory setup to final cleanup.

The following text-based tutorial provides an excellent walkthrough of a complete vector addition program. It reinforces the concepts we've just discussed and places them in the context of a full, working example.

Now, let's tie everything from this lesson and the previous one into a single, cohesive code structure.

The Full Workflow:

1. Host: Allocate memory for vectors h_A, h_B, h_C on the CPU using malloc. Initialize h_A and h_B.
2. Device: Allocate memory for vectors d_A, d_B, d_C on the GPU using cudaMalloc.
3. H -> D Copy: Copy input data from h_A and h_B to d_A and d_B using cudaMemcpy with cudaMemcpyHostToDevice.
4. Execute: Define launch parameters (threadsPerBlock, blocksPerGrid) and launch the vectorAdd kernel: vectorAdd<<<...>>>(d_A, d_B, d_C, N).
5. D -> H Copy: Copy the result from d_C back to h_C using cudaMemcpy with cudaMemcpyDeviceToHost.
6. Cleanup: Free all device memory using cudaFree and all host memory using free.
7. Verify: (Optional but recommended) Check that the results in h_C are correct.

The video below walks through this entire process from start to finish. Watching it will help solidify your understanding of how all these pieces fit together.

For reference, the gpuVectorSum.cu file detailed in the "Summation of Arrays on GPUs" section of the MolSSI Education resource (LINK) is a complete, runnable example that you can also study to see these concepts applied in a full program.

Conclusion

Congratulations! You have now learned how to write and launch a simple parallel algorithm on the GPU. We've bridged the gap between managing data and performing computations, giving you a complete end-to-end framework for CUDA programming.

Key Takeaways:

• Kernels: Functions marked with __global__ run in parallel on the GPU and are called from the host.
• Thread Indexing: The formula blockIdx.x * blockDim.x + threadIdx.x is the standard way for a thread to compute its unique global ID in a 1D grid.
• Kernel Launch: The <<<grid, block>>> syntax specifies the execution configuration, defining the total number of threads to be launched.
• Boundary Checks: Always include a check like if (thread_id < num_elements) inside your kernel to ensure you don't access memory out of bounds.

In this lesson, our kernel accessed memory from a single pool on the GPU called "global memory." However, the GPU provides several different types of memory with different characteristics (size, speed, scope). In our next lesson, we will describe the different GPU memory spaces (global, shared, constant, registers) and their use cases, which is the first step toward optimizing the performance of your kernels.