C++ Basics to GPU: From Hello World to CUDA

These are my working notes from learning C++ and then CUDA, collected into one reference. The first half covers the language basics and object orientation; the second half moves onto the GPU: kernels, the execution configuration, unified memory, error handling, multiprocessors, streams, and a few small projects.

Hello World and more

Preprocessor directives run before compilation. Use #include <iostream>, with <> for standard libraries and "" for user-defined ones. using namespace std; lets you use standard-library objects without the std:: prefix.

int main() returns an int and takes no arguments except command-line ones:

int main(int argc, char *argv[]) { }
// or
int main(int argc, char **argv) { }

argc is the number of arguments (default 1, the path to the executable), and argv holds them as C-strings, with argv[0] being the executable itself.

cout << x prints x and cin >> x reads into x:

int x, y;
std::cin >> x >> y;
std::cout << x << y;

const <type> var = val; behaves like const in JavaScript. The core datatypes and sizes:

int: 4 bytes, range $[-2^{31}, 2^{31}]$
double: 8 bytes, 15 decimal digits
char: 1 byte
bool: 1 byte
std::string: needs <string>; has .length() and indexes with []

Cast with (type), for example (double) x. Operators come in the usual families: logical &&, ||, !; bitwise << (left shift), >> (right shift), ~ (NOT), & (AND), | (OR), ^ (XOR); and the ternary ? :.

Control flow includes switch, do-while, and a range-based for:

switch (value) {
  case 1:
    cout << "Case 1";
    break;
  case 2:
    cout << "Case 2";
  default:
    cout << "Default";
}
 
do {
  cout << "Did";   // runs at least once
} while (false);
 
int fibonacci[5] = {0, 1, 1, 2, 3};
for (int number : fibonacci) {
  std::cout << number;
}

Auto and functions

auto lets the compiler deduce a variable's type:

int fibonacci[5] = {0, 1, 1, 2, 3};
for (auto number : fibonacci) {
  std::cout << number;
}

Functions declare a return type, arguments, and optional defaults (which must follow non-default parameters). Overloading lets several functions share a name as long as they differ in parameter types or count:

int add(int a, int b) {
  return a + b;
}
 
double add(double a, double b) {
  return a + b;
}
 
int main() {
  cout << add(3, 2);     // calls add(int, int)
  cout << add(5.3, 1.4); // calls add(double, double)
}

Data structures

Fixed-size arrays, including 2D:

char game[3][3] = {
  {'x', 'o', 'o'},
  {'o', 'x', 'x'},
  {'o', 'o', 'x'}
};

Vectors are arrays with a mutable size:

#include <vector>
 
std::vector<int> weights;
weights.push_back(25);
weights.push_back(45);
weights.pop_back();   // weights now holds {25}
// also: .front, .back, .size, .empty

The standard library also ships stacks, queues, sets, unordered sets, and hashmaps.

References

The address-of operator is &:

cout << message << endl;   // value: Hello World!
cout << &message << endl;  // address: 0x7ffee9b21af0

A reference is an alias for an existing variable, declared with &. References must be initialized and should not be reassigned:

int& length;         // invalid: must be initialized
int x = 5, y = 2;
int& xref = x;
xref = y;            // does not rebind: assigns y's value into x

Their main use is function arguments, since a function can modify a value passed by reference:

void change(int& x) {
  x = 3;
}
int main() {
  int x = 5;
  cout << x;   // 5
  change(x);
  cout << x;   // 3
}

Memory: stack and heap

C++ uses two main regions of memory, the stack and the heap.

Feature	Stack	Heap
Scope	Local variables within functions	Dynamically allocated variables
Lifetime	Automatic (freed when the function exits)	Manual (`free()` or `delete`)
Allocation	Automatic (compile-time)	Dynamic (runtime)
Size	Limited, often smaller	Much larger, prone to fragmentation
Speed	Faster	Slower
Example	`int x = 3;`	`int* ptr = malloc(sizeof(int));`

The stack is a highly ordered region that handles function calls and local variables, and it is released automatically when a function returns. The heap is an unordered region of dynamically allocated memory that you manage yourself with new and delete; it persists until you explicitly release it.

Allocate on the heap when memory must outlive a function's scope. With C-style malloc and free:

int* ptr = (int*)malloc(sizeof(int));  // uninitialized memory on the heap
*ptr = 3;
free(ptr);                             // prevents a leak

With C++ new and delete (which zero-initialize basic types):

int* ptr = new int;
*ptr = 3;
delete ptr;

And for arrays, new[] and delete[]:

int* arr = new int[10];
arr[0] = 3;
delete[] arr;

Object orientation

By default all class members are private. Mark them public to access them from outside:

class School {
public:
  std::string name;
  int age;
  void getName();
};
 
void School::getName() {
  std::cout << "Name";
}
 
School yeet;
yeet.name = "YeetCamp";
yeet.age = 21;
yeet.getName();

Access specifiers control visibility:

Access	`public`	`protected`	`private`
Inside the class	yes	yes	yes
Inside derived classes	yes	yes	no
Outside the class	yes	no	no

A constructor runs when an object is created; a destructor runs when it goes out of scope. A member initializer list is the way to initialize const members (which cannot use =):

class House {
private:
  std::string location;
  int rooms;
public:
  House(std::string loc = "New York", int num = 5) {
    location = loc;
    rooms = num;
  }
  void summary() {
    std::cout << location << " house with " << rooms << " rooms";
  }
  ~House() {
    std::cout << "Moved away from " << location;
  }
};
 
class Book {
private:
  const std::string title;
  const int pages;
public:
  Book() : title("Diary"), pages(100) {}   // member initializer list
};

Encapsulation

Getters and setters keep fields private while exposing controlled access:

class Clock {
private:
  int time = 1200;
public:
  int getTime() {
    return time;
  }
  void setTime(int new_time) {
    time = new_time;
  }
};
 
int main() {
  Clock alarm;
  alarm.setTime(930);
  std::cout << alarm.getTime();   // 930
}

Inheritance

A derived class calls its base constructor through its own initializer list:

class Animal {
private:
  std::string gender;
  int age;
public:
  Animal(std::string new_gender, int new_age)
    : gender(new_gender), age(new_age) {}
};
 
class Dog : public Animal {
private:
  std::string breed;
public:
  Dog(std::string new_gender, int new_age, std::string new_breed)
    : Animal(new_gender, new_age), breed(new_breed) {}
  void sound() {
    std::cout << "Woof\n";
  }
};
 
int main() {
  Dog buddy("male", 8, "Husky");
  buddy.sound();   // Woof
}

Inheritance chains (multilevel) construct base to derived:

class A { public: A() { std::cout << "Constructing A\n"; } };
class B : public A { public: B() { std::cout << "Constructing B\n"; } };
class C : public B { public: C() { std::cout << "Constructing C\n"; } };
 
int main() {
  C example;   // prints A, then B, then C
}

Polymorphism

A derived class can override a base method of the same name:

class Animal {
public:
  void action() { std::cout << "The animal does something.\n"; }
};
 
class Fish : public Animal {
public:
  void action() { std::cout << "Fish swim.\n"; }
};
 
class Bird : public Animal {
public:
  void action() { std::cout << "Birds fly.\n"; }
};
 
int main() {
  Animal newAnimal;
  Fish newFish;
  Bird newBird;
  newAnimal.action();
  newFish.action();
  newBird.action();
}

CUDA: the first kernel

On the GPU, functions are called kernels. Because moving data between host (CPU) and device (GPU) is expensive, it only pays off when the compute done per transferred element is high. A matrix addition transfers $3N^2$ elements for $N^2$ additions, a ratio of $O(1)$ , so it barely benefits; a matrix multiplication does $N^3$ work for the same transfer, a ratio of $O(N)$ , so larger matrices win more.

A __global__ function runs on the GPU and returns void. It can be launched from the host with an execution configuration <<<blocks, threads>>>, and the host waits for it with cudaDeviceSynchronize() since launches are asynchronous:

void CPUFunction() {
  printf("This function is defined to run on the CPU.\n");
}
 
__global__ void GPUFunction() {
  printf("This function is defined to run on the GPU.\n");
}
 
int main() {
  CPUFunction();
  GPUFunction<<<1, 1>>>();   // launch on the GPU
  cudaDeviceSynchronize();   // wait for it to finish
}

The configuration <<<A, B>>> launches A blocks of B threads each (max 1024 threads per block), and a grid is a collection of blocks. Inside a kernel you navigate the hierarchy with gridDim.x (blocks per grid), blockIdx.x (this block's index), blockDim.x (threads per block), and threadIdx.x (this thread's index). Production code should check each API call and call cudaGetLastError() after launches.

Compiling with nvcc

nvcc will feel familiar to gcc users. To compile and run a .cu file:

nvcc -o out some-CUDA.cu -run

The -o flag names the output binary, and the convenience flag -run executes it after a successful compile.

Unified memory and accelerated loops

A kernel can loop by using its thread index as an iteration number:

__global__ void loop(int N) {
  printf("This is iteration number %d\n", threadIdx.x);
}
 
int main() {
  int N = 10;
  loop<<<1, N>>>(N);
  cudaDeviceSynchronize();
}

The standard way to map threads across blocks to a flat data index is threadIdx.x + blockIdx.x * blockDim.x. Allocate unified memory with cudaMallocManaged and release it with cudaFree:

__global__ void initializeElementsTo(int initialValue, int *a, int N) {
  int i = threadIdx.x + blockIdx.x * blockDim.x;
  if (i < N) {
    a[i] = initialValue;
  }
}
 
int main() {
  int N = 1000;
  int *a;
  size_t size = N * sizeof(int);
  cudaMallocManaged(&a, size);
 
  initializeElementsTo<<<4, 256>>>(6, a, N);
  cudaDeviceSynchronize();
 
  cudaFree(a);
}

When the data is larger than the grid, a grid-stride loop lets each thread handle multiple elements, striding by the total number of threads gridDim.x * blockDim.x:

__global__ void doubleElements(int *a, int N) {
  int idx = blockIdx.x * blockDim.x + threadIdx.x;
  int gridStride = gridDim.x * blockDim.x;
  for (int i = idx; i < N; i += gridStride) {
    a[i] *= 2;
  }
}

Error handling

Most CUDA functions return a cudaError_t, which you should check against cudaSuccess:

cudaError_t err;
err = cudaMallocManaged(&a, N);
if (err != cudaSuccess) {
  printf("Error: %s\n", cudaGetErrorString(err));
}

A small macro keeps this tidy across many calls:

#include <stdio.h>
#include <assert.h>
 
inline cudaError_t checkCuda(cudaError_t result) {
  if (result != cudaSuccess) {
    fprintf(stderr, "CUDA Runtime Error: %s\n", cudaGetErrorString(result));
    assert(result == cudaSuccess);
  }
  return result;
}
 
int main() {
  checkCuda(cudaDeviceSynchronize());
}

Multi-dimensional data

Grids and blocks can have up to three dimensions using CUDA's dim3 type, which is handy for 2D matrices without any performance penalty:

dim3 threads_per_block(16, 16, 1);
dim3 number_of_blocks(16, 16, 1);
someKernel<<<number_of_blocks, threads_per_block>>>();

Here gridDim.x, gridDim.y, blockDim.x, and blockDim.y would all equal 16.

Project: thermal conductivity

A 2D heat-diffusion stencil, accelerated by mapping each grid point to a thread and comparing the GPU result against a CPU reference:

#include <stdio.h>
#include <math.h>
 
#define I2D(num, c, r) ((r)*(num)+(c))
 
__global__ void step_kernel_mod(int ni, int nj, float fact, float* temp_in, float* temp_out) {
  int i00, im10, ip10, i0m1, i0p1;
  float d2tdx2, d2tdy2;
 
  int row = threadIdx.x + blockIdx.x * blockDim.x;
  int col = threadIdx.y + blockIdx.y * blockDim.y;
 
  if (row > 0 && col > 0 && row < ni-1 && col < nj-1) {
    i00 = I2D(ni, row, col);
    im10 = I2D(ni, row-1, col);
    ip10 = I2D(ni, row+1, col);
    i0m1 = I2D(ni, row, col-1);
    i0p1 = I2D(ni, row, col+1);
 
    d2tdx2 = temp_in[im10] - 2*temp_in[i00] + temp_in[ip10];
    d2tdy2 = temp_in[i0m1] - 2*temp_in[i00] + temp_in[i0p1];
 
    temp_out[i00] = temp_in[i00] + fact*(d2tdx2 + d2tdy2);
  }
}
 
int main() {
  const int ni = 200, nj = 100, nstep = 200;
  float tfac = 8.418e-5;   // thermal diffusivity of silver
  const int size = ni * nj * sizeof(float);
 
  float *temp1, *temp2, *temp_tmp;
  cudaMallocManaged(&temp1, size);
  cudaMallocManaged(&temp2, size);
 
  dim3 threads(32, 16, 1);
  dim3 blocks((ni/threads.x)+1, (nj/threads.y)+1, 1);
 
  for (int istep = 0; istep < nstep; istep++) {
    step_kernel_mod<<<blocks, threads>>>(ni, nj, tfac, temp1, temp2);
    cudaDeviceSynchronize();
    temp_tmp = temp1; temp1 = temp2; temp2 = temp_tmp;   // swap buffers
  }
 
  cudaFree(temp1);
  cudaFree(temp2);
}

Profiling with nsys

nsys profile generates a report you can open in Nsight Systems. Passing --stats=true prints summary tables, including the CUDA API summary, the kernel summary, and the memory-time and memory-size operation summaries.

Streaming multiprocessors and warps

GPUs execute on streaming multiprocessors (SMs). Two portability-friendly heuristics help performance: choose a number of blocks that is a multiple of the number of SMs, and a block size that is a multiple of 32 (the warp size, the grouping of threads an SM schedules together). Rather than hard-code the SM count, query it at runtime:

int deviceId;
cudaGetDevice(&deviceId);
 
cudaDeviceProp props;
cudaGetDeviceProperties(&props, deviceId);   // props.major, props.minor, ...
 
int numberOfSMs;
cudaDeviceGetAttribute(&numberOfSMs, cudaDevAttrMultiProcessorCount, deviceId);

Project: optimized vector addition

Sizing the grid to the SM count and prefetching the unified memory turns a naive vector add into a much faster one:

#include <stdio.h>
 
__global__ void initWith(float num, float *a, int N) {
  int index = threadIdx.x + blockIdx.x * blockDim.x;
  int stride = blockDim.x * gridDim.x;
  for (int i = index; i < N; i += stride) a[i] = num;
}
 
__global__ void addVectorsInto(float *result, float *a, float *b, int N) {
  int index = threadIdx.x + blockIdx.x * blockDim.x;
  int stride = blockDim.x * gridDim.x;
  for (int i = index; i < N; i += stride) result[i] = a[i] + b[i];
}
 
int main() {
  int deviceId, numberOfSMs;
  cudaGetDevice(&deviceId);
  cudaDeviceGetAttribute(&numberOfSMs, cudaDevAttrMultiProcessorCount, deviceId);
 
  const int N = 2<<24;
  size_t size = N * sizeof(float);
 
  float *a, *b, *c;
  cudaMallocManaged(&a, size);
  cudaMallocManaged(&b, size);
  cudaMallocManaged(&c, size);
 
  cudaMemPrefetchAsync(a, size, deviceId);
  cudaMemPrefetchAsync(b, size, deviceId);
  cudaMemPrefetchAsync(c, size, deviceId);
 
  size_t threadsPerBlock = 256;
  size_t numberOfBlocks = 32 * numberOfSMs;
 
  initWith<<<numberOfBlocks, threadsPerBlock>>>(3, a, N);
  initWith<<<numberOfBlocks, threadsPerBlock>>>(4, b, N);
  initWith<<<numberOfBlocks, threadsPerBlock>>>(0, c, N);
 
  addVectorsInto<<<numberOfBlocks, threadsPerBlock>>>(c, a, b, N);
  cudaDeviceSynchronize();
 
  cudaMemPrefetchAsync(c, size, cudaCpuDeviceId);
 
  cudaFree(a);
  cudaFree(b);
  cudaFree(c);
}

Asynchronous memory prefetching

Prefetching migrates unified memory to a CPU or GPU in the background, before the code needs it, cutting page-fault and on-demand-migration overhead. It also moves data in larger chunks, so it fits best when access patterns are known ahead of time and not sparse:

int deviceId;
cudaGetDevice(&deviceId);
 
cudaMemPrefetchAsync(pointerToSomeUMData, size, deviceId);          // prefetch to GPU
cudaMemPrefetchAsync(pointerToSomeUMData, size, cudaCpuDeviceId);   // prefetch to host

Project: SAXPY

A single-precision a*x + b over large vectors, again sized to the SMs:

#include <stdio.h>
 
#define N 2048 * 2048
 
__global__ void saxpy(int *a, int *b, int *c) {
  int tid = blockIdx.x * blockDim.x + threadIdx.x;
  int stride = gridDim.x * blockDim.x;
  for (int i = tid; i < N; i += stride)
    c[i] = 2 * a[i] + b[i];
}
 
int main() {
  int deviceId, SMs;
  cudaGetDevice(&deviceId);
  cudaDeviceGetAttribute(&SMs, cudaDevAttrMultiProcessorCount, deviceId);
 
  int *a, *b, *c;
  int size = N * sizeof(int);
  cudaMallocManaged(&a, size);
  cudaMallocManaged(&b, size);
  cudaMallocManaged(&c, size);
 
  for (int i = 0; i < N; ++i) { a[i] = 2; b[i] = 1; c[i] = 0; }
 
  cudaMemPrefetchAsync(a, size, deviceId);
  cudaMemPrefetchAsync(b, size, deviceId);
  cudaMemPrefetchAsync(c, size, deviceId);
 
  saxpy<<<32 * SMs, 256>>>(a, b, c);
  cudaDeviceSynchronize();
 
  cudaFree(a); cudaFree(b); cudaFree(c);
}

Concurrent streams

Operations within a stream run in order, operations in different non-default streams have no ordering relative to each other, and the default stream blocks (it waits for all others and holds them up until it finishes). Creating your own streams lets independent work overlap:

cudaStream_t stream;
cudaStreamCreate(&stream);
 
someKernel<<<number_of_blocks, threads_per_block, 0, stream>>>();   // stream is the 4th argument
 
cudaStreamDestroy(stream);

Applied to the vector-add initializations, each initWith runs on its own stream:

cudaStream_t stream1, stream2, stream3;
cudaStreamCreate(&stream1);
cudaStreamCreate(&stream2);
cudaStreamCreate(&stream3);
 
initWith<<<numberOfBlocks, threadsPerBlock, 0, stream1>>>(3, a, N);
initWith<<<numberOfBlocks, threadsPerBlock, 0, stream2>>>(4, b, N);
initWith<<<numberOfBlocks, threadsPerBlock, 0, stream3>>>(0, c, N);
 
addVectorsInto<<<numberOfBlocks, threadsPerBlock>>>(c, a, b, N);
cudaDeviceSynchronize();
 
cudaStreamDestroy(stream1);
cudaStreamDestroy(stream2);
cudaStreamDestroy(stream3);

Project: n-body simulator

The capstone: an all-pairs gravitational n-body step. Each thread accumulates the force on one body from every other, then a second kernel integrates positions. Sizing blocks to the SMs and prefetching pushes the interactions-per-second number up:

#include <math.h>
#include <stdio.h>
#include <assert.h>
 
#define SOFTENING 1e-9f
 
typedef struct { float x, y, z, vx, vy, vz; } Body;
 
__global__ void bodyForce(Body *p, float dt, int n) {
  int index = threadIdx.x + blockIdx.x * blockDim.x;
  int stride = blockDim.x * gridDim.x;
 
  for (int i = index; i < n; i += stride) {
    float Fx = 0.0f, Fy = 0.0f, Fz = 0.0f;
    for (int j = 0; j < n; j++) {
      float dx = p[j].x - p[i].x;
      float dy = p[j].y - p[i].y;
      float dz = p[j].z - p[i].z;
      float distSqr = dx*dx + dy*dy + dz*dz + SOFTENING;
      float invDist = rsqrtf(distSqr);
      float invDist3 = invDist * invDist * invDist;
      Fx += dx * invDist3; Fy += dy * invDist3; Fz += dz * invDist3;
    }
    p[i].vx += dt*Fx; p[i].vy += dt*Fy; p[i].vz += dt*Fz;
  }
}
 
__global__ void integratePosition(Body *p, float dt, int n) {
  int index = threadIdx.x + blockIdx.x * blockDim.x;
  int stride = blockDim.x * gridDim.x;
  for (int i = index; i < n; i += stride) {
    p[i].x += p[i].vx*dt;
    p[i].y += p[i].vy*dt;
    p[i].z += p[i].vz*dt;
  }
}
 
int main(const int argc, const char** argv) {
  int device_id, num_of_sms;
  cudaGetDevice(&device_id);
  cudaDeviceGetAttribute(&num_of_sms, cudaDevAttrMultiProcessorCount, device_id);
 
  int nBodies = 2<<11;
  if (argc > 1) nBodies = 2<<atoi(argv[1]);
 
  const float dt = 0.01f;
  const int nIters = 10;
  int bytes = nBodies * sizeof(Body);
 
  float *buf;
  cudaMallocManaged(&buf, bytes);
  Body *p = (Body*)buf;
 
  size_t threadsPerBlock = 128;
  size_t numberOfBlocks = 32 * num_of_sms;
  cudaMemPrefetchAsync(buf, bytes, device_id);
 
  for (int iter = 0; iter < nIters; iter++) {
    bodyForce<<<numberOfBlocks, threadsPerBlock>>>(p, dt, nBodies);
    integratePosition<<<numberOfBlocks, threadsPerBlock>>>(p, dt, nBodies);
    cudaDeviceSynchronize();
  }
 
  cudaFree(buf);
}

That is the whole arc: from printf on the CPU to thousands of GPU threads integrating a gravitational system in parallel.