“Straightforward”
You create an int variable which is basically the thread index, since we are launching enough threads, each thread does the computation for index, index + stride, index + 2*stride and so on.
The stride is essentially the number of total number of threads launched which can be calculated by blockDim.x * gridDim.x where blockDim.x tells you how many threads are in the block and gridDim.x tells you how many blocks are in the grid.
In a strict sense every kernel launch is essentially a grid created that handles computation. More on the architecture of GPU and how it relates to grid, blocks and threads.
__global__ void vector_add(const float* A, const float* B, float* C, int N) {
int i = blockDim.x * blockIdx.x + threadIdx.x;
int stride = blockDim.x * gridDim.x;
//the parallelism is in the many threads; the loop is just a way to “tile”
//the work across them.
//the for loop is for a particular thread to handle computation
//for example if you dont use the stride, and increment it simply by
// 1 then essentially
//your entire computation is being overlapped by just each thread ie
//every thread is doing the addition computation for all N elements
for(int idx = i; idx<N; idx += stride){
C[idx] = A[idx] + B[idx];
}
}
“I didn’t want to do matmul next”
ReLU is essentially max(0, x), and given a 1D array of elements, all you have to do is write the same formula in the code. The rest of the kernel looks pretty much the same as the vector addition kernel. Since we launch enough threads -
thread 0 → 0, 0 + stride, 0 + 2*stride
thread 1 → 1, 1 + stride, 1 + 2*stride and so on
Similarly up until the nth thread. Like mentioned previously, the “for” loop in the kernel is to iterate through the elements that are computed per thread where the increment is equal to the stride. The parallelism is due to the fact that many threads are launched in parallel.
The initial kernel -
#include <cuda_runtime.h>
__global__ void relu_kernel(const float* input, float* output, int N) {
int i = blockDim.x * blockIdx.x + threadIdx.x;
int stride = blockDim.x * gridDim.x;
for(int idx = i; idx < N; idx += stride){
output[idx] = fmaxf(0.0f, input[idx]);
}
}
// input, output are device pointers (i.e. pointers to memory on the GPU)
extern "C" void solve(const float* input, float* output, int N) {
int threadsPerBlock = 256;
int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;
relu_kernel<<<blocksPerGrid, threadsPerBlock>>>(input, output, N);
cudaDeviceSynchronize();
}
// relu is simply max(0, x)
We use fmaxf for floating point values and fmax for double values. max is part of the cmath library in C++ and will not work inside a CUDA kernel. We may use fmax also for floating point values, I made that error initially and did not notice since the kernel passed the test results. Turns out that fmax will convert your float value to a double, perform the operation and then convert it back.