GPUs

Introduction

Continuing with the other posts¹ about ~how computers work at an ELI24 level, this post aims to understand GPU devices, how their architecture differs from CPUs (and why), as well as some rudimentary lessons in writing software for these devices.

I'll begin with an architectural comparison against the maybe-more-familiar CPU:

First of all, we can see that a GPU consists of a larger number of cores and a shallower cache hierarchy. It's worth noting that CPU cores, in general, are much more capable than those present in GPU devices, and the presence of additional caches vastly reduces memory access latency. The GPU paradigm of parallelization is evidenced at the hardware design level by the presence of singular control units (purple) responsible for entire thread groups of cores, implying that each core in a thread group must execute the same instruction.

The differences in architectures serve different purposes. CPUs are designed for low-latency, low-throughput serial operations, whereas GPUs are optimized for high-latency, high-throughput parallel operations.

Sample Benchmarks

We can gain some intuition about the ramifications of these hardware differences via benchmarks on modern devices^2,3:

† per Streaming Multiprocessor (a row of connected cores)

Parallel Operations

An example class of a parallel operation for which GPUs outperform CPUs is vector arithmetic:

An implementation of vector addition designed for execution on a CPU might look something like:

void add(int n, float* x, float* y, float* z) {
	for (int i = 0; i < n; i++)
		z[i] = x[i] + y[i];
}

this code iterates over every single element in each input vector x, y which will be $O(n)$. However, since each element in z is independent of any other element in z, we can see that this approach can be improved by parallelizing the summation and computing each entry simultaneously. The GPU code (which we call a kernel) which achieves this looks like:

__global__ void add(int n, float* x, float* y, float* z) {
	int i = blockIdx.x * blockDim.x + threadIdx.x;
	if (i < n)
		z[i] = x[i] + y[i];
}

Here, we use the __global__ keyword, which indicates to the compiler that this function should be executed on the GPU. The next bit of GPU-spice is the index variable i which uses "thread local" variables which would be contextually available in each instance of this function being run in parallel. The way CUDA achieves this is with the concept of thread blocks.⁴ When we launch our GPU kernel, we must specify the number blocks to launch and how many threads each block contains – this pattern informs how we index into all of our data structures to be accessed by the GPU. For example, if we launch $n/2$ blocks each containing two threads, the resulting thread organization will resemble:

Each thread executes the same code and is now able to compute the appropriate vector index via its own block- and thread-index.

The last notable difference between the two implementation is the existence of a bounds check on the computed index i. This is necessary in the event that the dimension of our input vectors does not align with our block dimension e.g. for an odd value of $n$, our block alignment would be:

With a GPU kernel in hand, we do still need actually delegate execution from the CPU. This is a multistep process following the procedure:

1. Allocate memory on the GPU
2. Copy the data to the GPU (since it doesn't have direct access itself)
3. Run the kernel
4. Copy the results back to the CPU
5. Free the memory

1. Allocate memory on the GPU

// we use the _d suffix to denote that this will be sent to the device (GPU)
float* a_d;
float* b_d;
float* c_d;
int N = 4096;
cudaMalloc((void**) &a_d, N * sizeof(float));
cudaMalloc((void**) &b_d, N * sizeof(float));
cudaMalloc((void**) &c_d, N * sizeof(float));

2. Copy the data to the GPU

cudaMemcpy(a_d, a, N * sizeof(float), cudaMemcpyHostToDevice);
cudaMemcpy(b_d, a, N * sizeof(float), cudaMemcpyHostToDevice);

this function is nearly identical to the C standard counterpart, with the added argument denoting the type or direction of the copy. Read more about it here.⁵

3. Run the kernel

to run the kernel we specify how many blocks, and how many threads in each block. CUDA exposes these configurables with some funky syntax,⁶ but it's basically just division.

As an example, a function declared as __global__ void Func(float* parameter); must be called like this: Func<<< Dg, Db, Ns >>>(parameter);

where Dg, Db, Ns are the dims of the grid, each block in the grid, and the number of bites in shared memory per block (defaulted to 0).

int BLOCK_SZ = 2048;
add<<<ceil(N / (float) BLOCK_SIZE), BLOCK_SIZE>>>(N, a_d, b_d, c_d);

4. Copy the results back to the CPU

Same as step 2, we just reverse the direction:

cudaMemcpy(c, c_d, N * sizeof(float), cudaMemcpyDeviceToHost);

5. Free the memory

Similar to standard memory management APIs, CUDA provides a method to allow us to mark those device-variables for re-allocation via cudaFree.

cudaFree(a_d);
cudaFree(b_d);
cudaFree(c_d);

Comparison

We can observe how the two implementations compare on each respective processing unit, noting that the GPU pays a penalty up front for fewer, slower caches, but remains leveled off for far-longer than the CPU, and also that for sufficiently large $n$, we start to pay exponential time penalties for CPU-based computation.

Comparing the ratio of execution time in $n$, this tradeoff becomes more evident:

Here we can see that that even our naive GPU kernel is nearly 200x faster for larger inputs. We can drive this number even higher by reducing memory accesses per operation.

Kernel Grids

Returning to a fixed instance of our vector addition kernel operating for a vector of size 6:

The code to delegate execution of our kernel to the GPU would resemble:

int N = 6;
int BLOCK_SIZE = 2;
add<<<ceil(N / (float) BLOCK_SIZE), BLOCK_SIZE>>>(N, a_d, b_d, c_d);

and our resulting kernel grid would be:

Not that the blockDim along the x dimension is constant for all blocks, so the only values which would change within our kernel code from earlier would be the block index (valued 0, 1, or 2) and the block's thread index (valued 0 or 1). We can observe then that our index value is a bijection onto the indices of our input vectors:

__global__ void add(int n, float* x, float* y, float* z) {
	int i = blockIdx.x * blockDim.x + threadIdx.x;
	if (i < n)
		z[i] = x[i] + y[i];
}

e.g. the computed index of the first thread in the second block would correspond to

int i = 1 * 2 + 0;

computing:

Our example kernel currently only references the .x dimension of our thread- and block-index variables, but CUDA allows us to configure multidimensional kernel grids. In some scenarios, we might prefer a two or even three dimensional kernel grid:

For $n > 2$-dimensional kernel grids, we specify the kernel parameters via dim3 types like so:

dim3 dimGrid(3,3,3);
dim3 dimbBlock(3,3,3);
add<<<dimGrid, dimBlock>>>(N, a_d, b_d, b_c);

This abstraction is mainly syntactic sugar for index-readability. Since arrays of arbitrary dimension are still stored as flattened contiguous chunks of memory, we still need to make use of the dimension variables in order compute our row/column/... indices. Indexing into a two-dimensional array via dimension variables is given by the formula:

A naive (oh, so terribly, terribly naive⁷) matrix multiplication kernel on a 2-dimensional kernel grid might then resemble:

__global__ void matmul(int n, float* a, float* b, float* c) {
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	if (row < n && col < n) {
		float dot_prod = 0.f;
		for (int i =0; i < n; i++) {
			dot_prod += a[row*n + i] * b[i*n + col];
		}
		c[row * n + col] = dot_prod;
	}
}

We could similarly achieve the same result with a 1-dimensional kernel grid like so:

__global__ void matmul(int n, float* a, float* b, float* c) {
	int idx = blockIdx.x * blockDim.x + threadIdx.x;
	int row = idx / n;
	int col = idx % n;
	if (row < n && col < n) {
		float dot_prod = 0.f;
		for (int i =0; i < n; i++) {
			dot_prod += a[row*n + i] * b[i*n + col];
		}
		c[row * n + col] = dot_prod;
	}
}

This adds a negligible amount of space overhead, as well as a bit more cognitive overhead, hence the preference for leveraging kernel grid syntactic sugar wherever possible.

This formula generalizes to higher dimensions as well like so:

Practical GPU usage: writing kernels for neural nets

Feed Forward

The feed forward step of an individual neural network $y_k$ neuron boils down to computing the dot product of an input vector $\vec{x}$ and network weights $\vec{w}$, plus some scalar bias $b_k$ associated with the neuron:

We can generalize the computations for all neurons by treating both the inputs and the weights as matrices, and the bias values as a vector:

The kernel for this computation is similar to the matrix multiplication kernel, except our matrices are no longer uniform/square, so we need to take additional params to describe their dimensions:

__global__ void forward(
	int batch_sz, int n, int out_w,
	float* inputs, float* weights, float* biases, float* outputs
) {
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	if (row < batch_sz && col < out_w) {
		// initialize the outputs to be the corresponding weights
		outputs[row * out_w + col] = biases[col];
		for (int i =0; i < n; i++) {
			outputs[row * out_w + col] +=
				weights[i * out_w + col] * inputs[row * n + i];
		}
	}
}

Activation

In order to prevent the network from collapsing to a single linear layer, we introduce a non-linearity such as ReLU:⁸

The $\text{ReLU}$ kernel is pretty straightforward:

__global__ void relu(int w, int h, float* inputs, float* outputs) {
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	if (row < h && col < w) {
		float activation = inputs[row * w + col];
		outputs[row * w + col] = activation > 0.f ? activation : 0.f;
	}
}

we apply this activation to every layer in the network barring the final output representing e.g. the probability that a given input corresponds to a given output.

Clamping

We need to translate a $\text{ReLU}$'d activation to a probability, such as:

The caveat of using an exponential translation like $\text{softmax}$ is that, for many inputs, we run the risk of numerical overflow since we're summing over many potentially-huge numbers. The workaround is to subtract the largest element of the input vector $\vec{x}$ that we're applying the $\text{softmax}$ to from the exponent of each element within:

The corresponding kernel (hopefully starting to look like a familiar pattern) is:

__global__ void softmax(int w, int h, float* inputs, float* outputs) {
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	if (row < h && col < w) {
		// find the maximum
		float max_val = intput[row * w];
		for (int i = 1; i < w; i++) {
			max_val = max(max_val, inputs[row * w + i])l
		}

		float divisor = 1e-9f; // prevent div by zero
		for (int i = 0; i. < w; i ++) {
			divisor += exp(inputs[row * w + i] - max_val);
		}
		outsputs[row * w + col] =
			exp(inputs[row * w + col] - max_val) / (divisor);
	}
}

Loss

The final function we might be interested in for a feed forward pass would be a loss function to measure the network's performance. A cross-entropy loss⁹ function takes two vectors and iterates over all $n$ possible output values, comparing the predicted results $\vec{p}$ of each to the actual output $\vec{q}$, assuming the outputs of the model are probabilities:

The lower our predicted probability is for a target score, the higher the loss:

Our loss kernel can be unidimensional since our output is a single dimension:

__global__ void cross_entropy_loss(
	int w, int h, float* preds, float* reals, float* outputs
) {
	int ix = blockIdx.x * blockDim.x + threadIdx.x;
	if (idx < h) {
		float loss = 0.f;
		for (int i = 0; i < w; i++) {
			loss -= real[idx * w + i] * log(max(1e-6, preds[idx * w + i]));
		}
		outputs[idx] = loss;
	}
}

Backpropagation

Backpropagation is the process of making incremental updates to the initially-random weights and biases of a neural network based on partial derivative of the loss with respect to the weight or bias of each neuron:

We can compute the partials at each layer by leveraging the chain rule. Starting with the last layer and the loss with respect the the activations at connecting it to the output, and working our way backwards multiplying each interior partial using the partials computed for the layer in front of it:

For our $\text{softmax}$ activation and cross entropy loss function $H$, we can compute the partials with respect to the inputs:

Since each element of the softmax depends on every other element in the input, the derivative will be a matrix. To simplify all the constituent partials, we make use of properties of logarithms to rearrange the inner expressions like so:

This lets us compute the derivative of a $\text{softmax}$ element into that element times the partial of the logarithm of that element. This is useful since it also lets us remove the division operator as well as one of the exponentiations:

Now, taking the partial derivative:

The first term simplifies to 1 if the element index matches the derivative index, otherwise it disappears:

as the more complicated second term, we can again apply the chain rule to simplify the computation for ourselves:

The term inside the summation depends on $x$ for only a single element: $x_k$, the rest are constants which will become zero, so we can further simplify:

So the complete equation for the partial derivative of $s_k$ becomes:

and we can plug this back into our matrix of softmax back propagation updates:

Now for the derivatives of our loss function w.r.t. the activations:

which is a bit simpler:

Multiplying the two collections, we attain the derivative of our loss w.r.t. the inputs:

and since $p$ is a probability distribution, the sum over all elements $(p_1 + p_2 + \cdots + p_n)$ will be 1, so we can simplify the resultant expression:

so our loss w.r.t. inputs is just the vector difference between our softmax'd outputs and the true distribution. Accordingly, the kernel to compute this is much friendlier than all the algebra above:

__global__ void cross_entropy_backwards(
	int w, int h, float* preds, float* reals, float* outputs
) {
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	if (row < h && col < w)
		outputs[row * w + col] = preds[row * w + col] - reals[row * w + col];
}

With the derivative of our loss w.r.t. the inputs from the output layers, we can compute each preceding derivative as well via the chain rule by multiplying the expression by the partial derivative of the inputs with respect to the activations of the preceding layer:

Recalling that our equation for the inputs to the previous layers is the matrix multiplication of the activations and the weights plus the biases:

Now here I've made an egregious mess of notation – the superscripts are indices not exponents. Taking the derivative of this expression yields just the weights of that previous layer

The corresponding kernel then just needs to perform matrix multiplication between the weights and the error from the next layer:

__global__ void backward(
	int batch_size, int n, int out_w,
	float* weights, float* biases, float* d_loss, float* out_d_loss
) {
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	if (row < batch_size && col < out_w) {
		float dloss = 0.f;
		for (int i = 0; i < n; i++) {
			float w = weights[i * out_w + col];
			dloss += w * d_loss[row * n + i];
		}
		out_d_loss[row * out_w + col] = dloss;
	}
}

The derivative of our chosen non-linearity is straightforward:

To apply backpropagation, we again apply the chain rule: multiplying by the partial derivative of our loss function w.r.t. the previous activation:

in kernel form:

__global__ void relu_backwards(
	int w, int h, int ns, float* a, float* d_loss, float* b
) {
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	if (row < w && col < h) {
		float activation = a[row * w + col]
		b[row * w + col] = activation > 0.f ? d_loss[row * w + col] : 0.f;
	}
}

Loss of weights and biases

To compute the error of our weights and biases, we take the partial derivative with respect to the desired quantity:

To update these values, we augment the existing weight or bias according to the derivative times some learning rate divided by our batch size:

and the last kernel to perform this gradient descent:

__global__ void update_layer(
	int w, int h, int batch_size, float η,
	float* weights, float* biases, float* activations, float* d_loss
) {
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	if (row < h && col < w) {
		float d_w = 0.f;
		float d_b = 0.f;
		for (int i = 0; i < batch_size; i++) {
			float activation = activations[i * h + row];
			float d_l = d_loss[i * w + col];
			dw += activation + d_l;
			db += d_l;
		}
		weights[row * w + col] -= η * d_w / batch_size;
		biases[col] -= η * d_b / batch_size;
	}
}

Performance for Measuring GPU Kernels

When we run code designed for the GPU, there's three primary devices that we're concerned with:

• the GPU itself
• the GPU's high bandwidth memory (DRAM)
• the CPU responsible for scheduling kernels and copying data to/from the HBM

Each of these devices has corresponding penalties referring to as compute time, memory access latency, and overhead respectively.

Continuing with the arbitrary example of executing an MNIST handwritten numerical digit, suppose we profiled our application and found that it takes 2s to load the data and 9.5s to train our neural net for 10 epochs:

We might be able to bring down the CPU overhead by using more sophisticated methods of optimized string parsing to ingest the MNIST csv data,¹⁰ bringing its overall time down to 1s, but we're still left with ~1s epochs:

The next step might be to observe that we don't need all of the data during the first epoch, just the data required for the first batch, so we can run these two tasks nearly simultaneously:

However, notice that this does not eliminate the overhead time, so even if we were (hypothetically *cough* *cough*) able optimize our GPU-based training computations, our first epoch would still be bottlenecked by the time required to load the first batch of training data:

CUDA API

A brief note about the CUDA API which helps understand precisely how we're able to parallelize these first two tasks. By design, CUDA is not a synchronous API, so when we run a kernel with the method<<< ... >>>(params) code, the CPU does not await the results of those tasks. Calling additional kernels adds them to a queue of kernels to be scheduled by the CPU which is managed by CUDA drivers for us. Making a synchronous/CPU-bound call after a kernel has been scheduled does wait for the scheduled kernels to complete before executing.

E.g. using the neural net methods implemented above, our training loop might look something like this:

load_mnist();

for (int epoch = 0; epoch < N_EPOCHS; epoch++) {
	for (int batch = 0; batch < train_length / BATCH_SIZE; batch++) {
		forward<<dimGrid, dimBlock>>>(...);
		relu<<dimGrid, dimBlock>>>(...);
		forward<<dimGrid, dimBlock>>>(...);
		relu<<dimGrid, dimBlock>>>(...);
		forward<<dimGrid, dimBlock>>>(...);
		softmax<<dimGrid, dimBlock>>>(...);
		cross_entropy_loss<<dimGrid, dimBlock>>>(...);

		cudaMemcpy(out_h, out_d, out_size, cudaMemcpyDeviceToHost);
		compute_accuracy(...);

		cross_entropy_backwards<<<<dimGrid, dimBlock>>>(...);
		backward<<dimGrid, dimBlock>>>(...);
		relu_backwards<<dimGrid, dimBlock>>>(...);
		backward<<dimGrid, dimBlock>>>(...);
		relu_backwards<<dimGrid, dimBlock>>>(...);

		update_layer<<dimGrid, dimBlock>>>(...);
		update_layer<<dimGrid, dimBlock>>>(...);
		update_layer<<dimGrid, dimBlock>>>(...);
	}
}

So, in order to parallelize the training and loading tasks, we can first schedule all of our kernels up front:

load_mnist();

for (int epoch = 0; epoch < N_EPOCHS; epoch++) {
	for (int batch = 0; batch < train_length / BATCH_SIZE; batch++) {
		forward<<dimGrid, dimBlock>>>(...);
		relu<<dimGrid, dimBlock>>>(...);
		forward<<dimGrid, dimBlock>>>(...);
		relu<<dimGrid, dimBlock>>>(...);
		forward<<dimGrid, dimBlock>>>(...);
		softmax<<dimGrid, dimBlock>>>(...);
		cross_entropy_loss<<dimGrid, dimBlock>>>(...);

		cross_entropy_backwards<<<<dimGrid, dimBlock>>>(...);
		backward<<dimGrid, dimBlock>>>(...);
		relu_backwards<<dimGrid, dimBlock>>>(...);
		backward<<dimGrid, dimBlock>>>(...);
		relu_backwards<<dimGrid, dimBlock>>>(...);

		update_layer<<dimGrid, dimBlock>>>(...);
		update_layer<<dimGrid, dimBlock>>>(...);
		update_layer<<dimGrid, dimBlock>>>(...);

		cudaMemcpy(out_h, out_d, out_size, cudaMemcpyDeviceToHost);
		compute_accuracy(...);
	}
}

Then, instead of loading the full data set, we can instead define a batch-loader to be called once before scheduling our kernels, and then call it again before our synchronization point:

load_next_batch();

for (int epoch = 0; epoch < N_EPOCHS; epoch++) {
	for (int batch = 0; batch < train_length / BATCH_SIZE; batch++) {
		forward<<dimGrid, dimBlock>>>(...);
		relu<<dimGrid, dimBlock>>>(...);
		forward<<dimGrid, dimBlock>>>(...);
		relu<<dimGrid, dimBlock>>>(...);
		forward<<dimGrid, dimBlock>>>(...);
		softmax<<dimGrid, dimBlock>>>(...);
		cross_entropy_loss<<dimGrid, dimBlock>>>(...);

		cross_entropy_backwards<<<<dimGrid, dimBlock>>>(...);
		backward<<dimGrid, dimBlock>>>(...);
		relu_backwards<<dimGrid, dimBlock>>>(...);
		backward<<dimGrid, dimBlock>>>(...);
		relu_backwards<<dimGrid, dimBlock>>>(...);

		update_layer<<dimGrid, dimBlock>>>(...);
		update_layer<<dimGrid, dimBlock>>>(...);
		update_layer<<dimGrid, dimBlock>>>(...);

		if (epoch == 0) load_next_batch();

		cudaMemcpy(out_h, out_d, out_size, cudaMemcpyDeviceToHost);
		compute_accuracy(...);
	}
}

GPU perf

When evaluating and improving code targeting the GPU, we're interested in the two other metrics besides CPU overhead. Our GPU compute time is bound by the throughput of the datatype being used on the device, typically measured in floating point operations per second (or FLOPS). The 3090 Ti³ we've been using as a benchmark has a theoretical limit of 40 TFLOPS on 32-bit floating point numbers, and a memory bandwidth of 1.01 TB/s. Therefore, we can perform computations much faster than we can actually load data from memory, so we might introduce a third meta-metric of computation intensity to measure how many FLOPS we can execute per byte loaded from memory.

To find the optimal performance ratio, we can plot the theoretical FLOPs against computational intensity to find the desired intersection:

If we're on the right side, we're compute bound, the GPU is working at full power. This is good. Here, we gain performance improvements by using datatypes with better OPS, optimizing our algorithms, or using more expensive/sophisticated hardware such as tensor cores for a higher throughput limit.

The left is bad, the GPU is being underworked, waiting for data from memory. To move from left to right, we can use different kinds of memory, or do some kernel fusion.

Kernel fusion

Consider one of the numerous bottlenecks in our relatively naive CUDA implementation of a neural network:

__global__ void relu(int w, int h, float* inputs, float* outputs) {
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	if (row < h && col < w) {
		float activation = inputs[row * w + col];
		outputs[row * w + col] = activation > 0.f ? activation : 0.f;
	}
}

This function has two memory access (reading a float from inputs and writing a float to outputs) both of size 4 bytes. We're also only doing two operations: the comparison and assignment of activation. This gives us a computation intensity of $\frac{2 \text{ FLOPs}}{8 \text{ B}}$ (not very intense imo).

Since our relu function is exclusively used after invoking the forward function, we can fuse these two functions together, completely eliminating the need for the memory accesses and operations added by relu:

__global__ void forward_relu(
	int batch_sz, int n, int out_w,
	float* inputs, float* weights, float* biases, float* outputs
) {
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	if (row < batch_sz && col < out_w) {
		float out = biases[col];
		outputs[row * out_w + col] = biases[col];
		for (int i =0; i < n; i++) {
			out += weights[i * out_w + col] * inputs[row * n + i];
		}
		outputs[row * out_w + col] = out > 0.f ? out : 0.f;
	}
}

We can similarly fuse the backwards pass. This achieves the "hypothetical" improvements alluded to earlier, reducing out total execution time by a factor of 2. At the moment, the rinky-dink example is hardly utilizing the GPU. Since our data in network is very small, the kernels aren't even filling all possible threads. If we increase the batch_size from 16 to 64, we can achieve an even higher rate of parallelization, bringing our per-epoch training time down from ~half a second to a tenth of a second, yielding a 5x improvement over the initial implementation.

GPU Memory Hierarchy

With respect to programming for a GPU with CUDA, there are two kinds of memory that we're concerned with: on- and off-chip. There are some segments of memory connected to the actual GPU chip, as well as smaller amounts of memory within the chip itself. Each thread within a block is able to access the global VRAM which resides off-chip:

Memory allocated via e.g. cudaMalloc(...) or statically initialized globally (with the __device__ keyword) lives in global memory which is the largest and slowest memory pool the GPU can access. Each thread also has access to a local set of on-chip registers. Variables created within our kernels are stored in these registers. Of course, if we initialize too many variables within our kernels (more than ~255 for modern commercial hardware), we induce kernel register occupancy, which will cause values in our registers to be evicted to "local memory" which is actually off-chip ("local" indicates that the variables are local to that thread).

We also have a small (on the order of 64 KB) section of VRAM called "constant memory" which is off-chip, cached, and read-only. Access to different addresses within this special segment are serialized, meaning that if we accesses the same address from multiple threads, we can achieve faster reads than if we made identical accesses to global memory.

To use constant memory, we use the __constant__ keyword when declaring a variable, and then use the cudaMemcpyToSymbol function to move data from the CPU to constant memory.

Finally, we have shared memory which is block-local memory:

Shared memory is on-chip, and therefore preferable to local, constant, and global memory, especially when the data be stored would be used by multiple threads. We specify __shared__ when declaring kernel variables in order to store them in this region.

Optimizing w.r.t. memory layout

With deeper knowledge of the memory hierarchy on our GPU, we can return to one of our previous kernels to see how we might take advantage of the layout in order to milk out some performance gains.

Recall the naive matrix multiplication kernel from earlier:

__global__ void matmul(int n, float* a, float* b, float* c) {
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	if (row < n && col < n) {
		float dot_prod = 0.f;
		for (int i =0; i < n; i++) {
			dot_prod += a[row*n + i] * b[i*n + col];
		}
		c[row * n + col] = dot_prod;
	}
}

For two $n \times n$ square matrices, this code will access each row of a $n$ times and each column of b ntimes, storing the dot product between corresponding elements in the thread's registers.