Writing Your First Kernel: Step-by-Step Introduction to CUDA Programming#

A. Why GPUs for General-Purpose Computing?#

GPUs were originally crafted to manage massive amounts of graphical computations in parallel. Unlike Central Processing Units (CPUs), which handle general-purpose serial tasks, GPUs are designed with a large number of smaller, more specialized cores that excel at doing multiple operations in parallel. This architecture proves beneficial for workloads involving:

• Matrix operations (e.g., linear algebra, deep learning)
• Image processing or computer vision
• Scientific computations (e.g., molecular dynamics, fluid simulations)
• Data analytics and machine learning

B. The Emergence of CUDA#

NVIDIA noticed that scientific and engineering communities, as well as data scientists, were repurposing GPUs for more general calculations. To simplify development, they introduced CUDA—an extension to the C/C++ language (with bindings to other languages too) that provides specific syntax and APIs to program GPUs. With CUDA, you can define special functions (kernels) that run in parallel on thousands of GPU threads to speed up your workloads.

C. High-Level Brainstorm: CPU vs. GPU#

Before diving into syntactical details, it is useful to conceptualize the distinction between CPU and GPU execution flows:

• CPU:

  • Few cores.
  • Optimized for serial tasks.
  • Large cache memory.
  • High clock speed.

• GPU:

  • Thousands of cores in parallel.
  • Optimized for throughput.
  • Smaller cache per core, but huge aggregate memory bandwidth.
  • Ideal for data-parallel tasks.

A. Kernels and Parallel Execution#

The programming model for CUDA is straightforward in concept:

1. You write functions called kernels that describe operations to be done on the GPU.
2. You launch these kernels from your CPU (host) code.
3. The GPU (device) runs these kernels in parallel across many threads.

A kernel is marked in C/C++ code with the __global__ keyword. For example:

1
__global__ void myKernel() {
2
    // Code executed on the GPU
3
}

When you call myKernel<<<blocks, threads>>>() from your host code, you specify how many blocks and threads per block you want. The GPU then launches that many instances of myKernel().

B. Threads, Blocks, and Grids#

CUDA organizes parallel execution hierarchically:

1. A grid consists of one or more blocks.
2. Each block consists of one or more threads.

This structure provides flexibility, allowing your program to scale from smaller GPUs to very large ones without necessarily being hardcoded for a specific number of cores.

C. Memory Model#

The CUDA memory model can be visualized as follows:

• Global memory: Large memory visible to all threads, but slower to access.
• Shared memory: A smaller, high-speed memory shared among threads in the same block.
• Registers: Fastest memory, private to each thread.
• Constant and texture memory: Special memories for read-only or specialized access patterns.

Understanding when and how to use different memory spaces is essential for achieving high performance. We’ll explore each in detail later.

A. Kernel Definition#

1
#include <iostream>
2

3
// Define a simple kernel
4
__global__ void helloKernel() {
5
    // Each thread prints a message.
6
    // This can lead to interleaved prints if using multiple threads.
7
    printf("Hello from thread %d\n", threadIdx.x);
8
}
9

10
int main() {
11
    // Launch the kernel with 1 block of 8 threads
12
    helloKernel<<<1, 8>>>();
13

14
    // Wait for GPU to finish before accessing results (on desktop this is usually automatic,
15
    // but better practice is to use cudaDeviceSynchronize)
16
    cudaDeviceSynchronize();
17

18
    return 0;
19
}

1. We include <iostream> (for C++ output streams) and rely on <cstdio> for printf.
2. The kernel __global__ void helloKernel() will run on the GPU, but it’s launched from the host.
3. helloKernel<<<1, 8>>>(); indicates that we want 8 threads in a single block.
4. threadIdx.x is a built-in CUDA variable indicating the current thread’s index within the block.

B. Compilation and Execution#

• Compile: nvcc hello.cu -o hello
• Run: ./hello

You should see:

1
Hello from thread 0
2
Hello from thread 1
3
...
4
Hello from thread 7

Although this example does not process any data, you’ve successfully compiled and executed your first CUDA kernel!

A. Thread Indexing#

Within a kernel, you can access built-in variables to determine your thread’s position:

• threadIdx.x, blockIdx.x – The index of the current thread or block in the x dimension.
• blockDim.x – The number of threads in each block.
• gridDim.x – The number of blocks in the grid.

If you want a unique global thread ID, you can do:

1
int globalThreadId = blockIdx.x * blockDim.x + threadIdx.x;

B. Multi-Dimensional Indexing#

You can also specify kernels along multiple dimensions (up to 3D for blocks and 3D for grids). For example:

1
dim3 threadsPerBlock(16, 16);
2
dim3 numBlocks(32, 32);
3
myKernel<<<numBlocks, threadsPerBlock>>>();

Inside myKernel, you can reference threadIdx.y, blockIdx.y, and so on. This flexibility is particularly useful for image processing (2D) or volumetric data (3D).

C. Determining Grid and Block Sizes#

Choosing the right threads-per-block (also known as block size) and the number of total threads can significantly impact performance. Typical guidance includes:

• Use a multiple of the warp size (32 threads) whenever possible.
• Consider the hardware constraints: For instance, some GPUs support up to 1024 threads per block.
• Generally, each block runs on one Streaming Multiprocessor (SM), so you often want enough blocks to keep all SMs busy.

A. Example: Allocating and Copying Data#

1
#include <iostream>
2

3
__global__ void incrementKernel(int* d_data) {
4
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
5
    d_data[idx] += 1;
6
}
7

8
int main() {
9
    const int arraySize = 256;
10
    const int bytes = arraySize * sizeof(int);
11

12
    // Allocate host memory
13
    int* h_data = (int*)malloc(bytes);
14
    // Initialize
15
    for(int i = 0; i < arraySize; i++) {
16
        h_data[i] = i;
17
    }
18

19
    // Allocate device memory
20
    int* d_data;
21
    cudaMalloc(&d_data, bytes);
22

23
    // Copy from host to device
24
    cudaMemcpy(d_data, h_data, bytes, cudaMemcpyHostToDevice);
25

26
    // Kernel launch
27
    int threads = 128;
28
    int blocks = arraySize / threads;
29
    incrementKernel<<<blocks, threads>>>(d_data);
30

31
    // Copy results back to host
32
    cudaMemcpy(h_data, d_data, bytes, cudaMemcpyDeviceToHost);
33

34
    // Check results
35
    for(int i = 0; i < 5; i++) {
36
        std::cout << "h_data[" << i << "] = " << h_data[i] << std::endl;
37
    }
38

39
    // Cleanup
40
    cudaFree(d_data);
41
    free(h_data);
42

43
    return 0;
44
}

In this code, each element gets incremented by 1. After copying the data back to the host, you can verify the correctness of the GPU operations.

B. Data Transfer Overheads#

Keep in mind that transferring data across the PCIe bus between the host and device can be a bottleneck. Strategies to minimize or overlap data transfers are essential in high-performance applications.

A. Global Memory#

• Large capacity.
• Slowest to access (compared to shared or registers).
• Persistent across kernel launches if allocated outside the kernel.

B. Shared Memory#

• On-chip, significantly faster than global memory.
• Only accessible by threads within the same block.
• Size is limited (e.g., 48 KB to 100 KB or more per SM, depending on GPU architecture).

Example usage of shared memory:

1
__global__ void sharedMemoryExample(float* d_input, float* d_output) {
2
    extern __shared__ float s_data[];
3

4
    int idx = threadIdx.x + blockDim.x * blockIdx.x;
5
    // Load data from global memory to shared memory
6
    s_data[threadIdx.x] = d_input[idx];
7
    __syncthreads();
8

9
    // Perform some computation
10
    s_data[threadIdx.x] *= 2.0f;
11
    __syncthreads();
12

13
    // Write back to global memory
14
    d_output[idx] = s_data[threadIdx.x];
15
}

extern __shared__ allows you to dynamically specify the size of shared memory at kernel launch:

1
int sharedMemSize = blockSize * sizeof(float);
2
sharedMemoryExample<<<gridSize, blockSize, sharedMemSize>>>(d_input, d_output);

C. Registers#

Registers are the fastest form of memory on the GPU, private to each thread. They store immediate values and variables needed for calculations. While incredibly fast, their limited size means the compiler will spill values to local memory if you use too many.

D. Constant and Texture Memory#

These are specialized forms of memory:

• Constant Memory is read-only from kernels but can be written from the host. Cached for broadcast usage.
• Texture Memory is specialized for 2D spatial locality. May be beneficial for certain image processing tasks.

A. Streams#

A stream in CUDA is essentially a queue of operations that execute in order on the GPU. By using multiple streams, you can overlap:

1. Kernel execution on one stream.
2. Data transfer on another stream.

Example:

1
cudaStream_t stream1, stream2;
2
cudaStreamCreate(&stream1);
3
cudaStreamCreate(&stream2);
4

5
// Launch kernel in stream1
6
myKernel<<<grid, block, 0, stream1>>>(d_dataA);
7
// Copy data in stream2
8
cudaMemcpyAsync(d_dataB, h_dataB, size, cudaMemcpyHostToDevice, stream2);
9

10
// Wait for streams to finish
11
cudaStreamSynchronize(stream1);
12
cudaStreamSynchronize(stream2);
13

14
cudaStreamDestroy(stream1);
15
cudaStreamDestroy(stream2);

B. Overlapping Data Transfers and Computation#

By default, CUDA operations on the default stream (NULL stream) are serialized. Leveraging multiple streams can yield substantial performance improvements when properly orchestrated.

A. Thread Synchronization#

Because GPU threads run in parallel, certain situations require synchronization:

1. Block-level synchronization: Using __syncthreads(), all threads in a block must reach the synchronization point before continuing.
2. Grid-level synchronization: Typically, you must use separate kernel launches or sophisticated techniques.

B. Atomic Operations#

When multiple threads access and modify shared data concurrently, you need atomic operations to avoid race conditions. CUDA provides atomic functions (e.g., atomicAdd, atomicMin, atomicMax) in global and shared memory.

1
__global__ void atomicExample(int* data) {
2
    int idx = threadIdx.x + blockIdx.x * blockDim.x;
3
    atomicAdd(&data[0], idx);
4
}

Use atomic operations judiciously: they can serialize parts of your code and reduce performance if overused.

A. Unified Memory#

Unified Memory simplifies memory management by automatically migrating data between CPU and GPU. You allocate memory with cudaMallocManaged():

1
int* data;
2
cudaMallocManaged(&data, n * sizeof(int));
3
// Use data on CPU or GPU without explicit cudaMemcpy

However, behind the scenes, data movements still occur, so performance may vary. Unified Memory is invaluable for rapid development and simpler code but may not achieve the same performance as carefully optimized explicit transfers.

B. Constant Memory#

Constant memory is beneficial for read-only data that is often accessed by many threads. You declare it with the __constant__ keyword:

1
__constant__ float constData[256];

Then on the host:

1
float hostData[256] = { /* ... initialize ... */ };
2
cudaMemcpyToSymbol(constData, hostData, sizeof(float)*256);

Access in kernel:

1
__global__ void useConstantMemory() {
2
    float value = constData[threadIdx.x];
3
    // ...
4
}

C. Texture Memory#

Texture memory is designed for spatial locality and can provide some caching benefits. It’s often used in graphics or image processing tasks, where threads in close vicinity read from nearby memory addresses.

A. Tools for Debugging#

• cuda-gdb: A GPU-capable version of the GDB debugger.
• NVIDIA Nsight Eclipse Edition: Offers a visual environment for debugging and profiling on Linux and other platforms.
• NVIDIA Nsight Visual Studio Edition: Provides an integrated solution for Windows developers.

B. Profiling with Nsight Systems and Nsight Compute#

• Nsight Systems: Visualize CPU/GPU timelines, identify execution bottlenecks, and see how kernels overlap with memory transfers.
• Nsight Compute: In-depth performance analysis of individual kernels (memory bandwidth, occupancy, etc.).

A. Vector Addition Example#

1
__global__ void vectorAdd(float* A, float* B, float* C, int n) {
2
    int i = blockIdx.x * blockDim.x + threadIdx.x;
3
    if (i < n) {
4
        C[i] = A[i] + B[i];
5
    }
6
}
7

8
int main() {
9
    int n = 1 << 20; // 1 million elements
10
    size_t bytes = n * sizeof(float);
11

12
    // Allocate host memory
13
    float *h_A = (float*)malloc(bytes);
14
    float *h_B = (float*)malloc(bytes);
15
    float *h_C = (float*)malloc(bytes);
16

17
    // Initialize data
18
    for (int i = 0; i < n; i++) {
19
        h_A[i] = 1.0f;
20
        h_B[i] = 2.0f;
21
    }
22

23
    // Allocate device memory
24
    float *d_A, *d_B, *d_C;
25
    cudaMalloc(&d_A, bytes);
26
    cudaMalloc(&d_B, bytes);
27
    cudaMalloc(&d_C, bytes);
28

29
    // Copy data to device
30
    cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice);
31
    cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice);
32

33
    // Launch kernel
34
    int blockSize = 256;
35
    int gridSize = (n + blockSize - 1) / blockSize;
36
    vectorAdd<<<gridSize, blockSize>>>(d_A, d_B, d_C, n);
37

38
    // Copy results back
39
    cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);
40

41
    // Verify
42
    for (int i = 0; i < 5; i++) {
43
        std::cout << "C[" << i << "] = " << h_C[i] << std::endl;
44
    }
45

46
    // Cleanup
47
    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
48
    free(h_A); free(h_B); free(h_C);
49
    return 0;
50
}

B. Matrix Multiplication Example#

Matrix multiplication is a classic example that benefits from GPU parallelism. For instance, we can multiply two square matrices A and B to produce C:

1
__global__ void matMulKernel(float* A, float* B, float* C, int n) {
2
    int row = blockIdx.y * blockDim.y + threadIdx.y;
3
    int col = blockIdx.x * blockDim.x + threadIdx.x;
4

5
    if(row < n && col < n) {
6
        float sum = 0.0f;
7
        for(int k = 0; k < n; k++) {
8
            sum += A[row * n + k] * B[k * n + col];
9
        }
10
        C[row * n + col] = sum;
11
    }
12
}
13

14
int main() {
15
    int n = 512;
16
    size_t bytes = n * n * sizeof(float);
17

18
    // Alloc and init host data
19
    float* h_A = (float*)malloc(bytes);
20
    float* h_B = (float*)malloc(bytes);
21
    float* h_C = (float*)malloc(bytes);
22

23
    for(int i = 0; i < n*n; i++) {
24
        h_A[i] = 1.0f;
25
        h_B[i] = 2.0f;
26
    }
27

28
    // Device
29
    float *d_A, *d_B, *d_C;
30
    cudaMalloc(&d_A, bytes);
31
    cudaMalloc(&d_B, bytes);
32
    cudaMalloc(&d_C, bytes);
33

34
    cudaMemcpy(d_A, h_A, bytes, cudaMemcpyHostToDevice);
35
    cudaMemcpy(d_B, h_B, bytes, cudaMemcpyHostToDevice);
36

37
    // Launch kernel
38
    dim3 threads(16, 16);
39
    dim3 blocks((n + threads.x - 1) / threads.x, (n + threads.y - 1) / threads.y);
40
    matMulKernel<<<blocks, threads>>>(d_A, d_B, d_C, n);
41

42
    cudaMemcpy(h_C, d_C, bytes, cudaMemcpyDeviceToHost);
43

44
    // Check a few values
45
    for(int i = 0; i < 5; i++) {
46
        std::cout << "C[" << i << "] = " << h_C[i] << std::endl;
47
    }
48

49
    // Free
50
    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
51
    free(h_A); free(h_B); free(h_C);
52

53
    return 0;
54
}

This naive matrix multiplication can be optimized using shared memory to reuse sub-blocks of the matrices, drastically reducing memory fetches from global memory.

A. Code Structuring#

1. Separate Device and Host Code

   • Keep kernel definitions in .cu files.
   • Create separate headers for function prototypes to facilitate large-scale projects.

2. Error Checking

   • Always check the return status of CUDA calls.
   • Use macros to simplify error handling. For example:

     1
     #define CUDA_CHECK(ans) { gpuAssert((ans), __FILE__, __LINE__); }
     2
     inline void gpuAssert(cudaError_t code, const char *file, int line) {
     3
         if (code != cudaSuccess) {
     4
             fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
     5
             exit(code);
     6
         }
     7
     }
     

B. Using C++ Features#

While CUDA can be combined with modern C++ features (C++11 or later), there are specific limitations regarding templates, inheritance, and certain class functionalities on the device side. Carefully check NVCC documentation for compatibility.

C. Profiling and Micro-Optimizations#

1. Occupancy

   • The ratio of active warps to the maximum number of warps supported on an SM.
   • Use occupancy calculators (Nsight) to tune block sizes.

2. Loop Unrolling

   • The compiler can unroll small loops to reduce loop overhead.
   • #pragma unroll can be used to guide unrolling.

3. Bank Conflicts in Shared Memory

   • Shared memory is divided into banks. Accessing the same bank by multiple threads concurrently can cause serialization. Reorganize data to avoid conflicts.

D. Multi-GPU and MPI#

For very large-scale problems, you might employ multiple GPUs:

• Use CUDA-aware MPI or peer-to-peer memory access (P2P) on multi-GPU systems.
• Distribute data across nodes in a cluster environment.

E. Libraries and Frameworks#

There is a rich ecosystem of GPU libraries:

• cuBLAS for linear algebra.
• cuFFT for Fast Fourier Transforms.
• Thrust (C++ STL-like library) for vectors, sorting, and more.
• cuDNN for deep learning primitives.

Leveraging these libraries often means you can sidestep writing and optimizing your own GPU kernels for standard operations.