simple vector addition code for CUDA. with rich comments

vector_add.cu

#include <chrono>
#include "../common/book.h"

#define N (1 << 16)
// element count: 65536

// max blocks -> 65535
// max threads per block -> 512

using namespace std;
using namespace std::chrono;

__global__
void add(const int* a, const int* b, int*c)
{
    // block layout
    // +----+----+----+----+----+
    // | B0 | T0 | T1 | T2 | T3 |
    // +----+----+----+----+----+
    // | B1 | T0 | T1 | T2 | T3 |
    // +----+----+----+----+----+
    // | B2 | T0 | ██ | T2 | T3 |
    // +----+----+----+----+----+
    // | B3 | T0 | T1 | T2 | T3 |
    // +----+----+----+----+----+
    //
    // blockDim.x = 4
    // so ██ will be executed when
    // threadIdx.x == 1 & blockIdx.x == 2
    //
    // our vectors are N elements long
    //     +---+---+---+---+
    //  A  | 0 | 1 | 2 | 3 |
    //     +---+---+---+---+
    //
    //     +---+---+---+---+
    //  B  | 0 | 1 | 4 | 9 |
    //     +---+---+---+---+
    //
    //     +----+----+----+----+
    //  C  | 0  | 2  | 6  | 12 |
    //     +----+----+----+----+

    unsigned int idx = threadIdx.x + blockIdx.x * blockDim.x;  // on the CPU we would start at 0
    while (idx < N) {
        c[idx] = a[idx] + b[idx];
        idx += blockDim.x * gridDim.x;  // and increment by 1,
        // but we are on the GPU, so we increment by the number of threads per block *
        // the number of blocks in the grid. This is simply the number of threads per
        // block multiplied by the number of blocks in the grid, or blockDim.x * gridDim.x
    }
}


int main()
{
    auto ts_start = high_resolution_clock::now();

    int a[N], b[N], c[N];
    int *dev_a, *dev_b, *dev_c;

    HANDLE_ERROR(cudaMalloc((void**)&dev_a, N * sizeof(int)));
    HANDLE_ERROR(cudaMalloc((void**)&dev_b, N * sizeof(int)));
    HANDLE_ERROR(cudaMalloc((void**)&dev_c, N * sizeof(int)));

    for (int i = 0; i < N; i++) {
        a[i] = i;
        b[i] = i * i;
    }

    HANDLE_ERROR(cudaMemcpy(dev_a,
                            a,
                            N * sizeof(int),
                            cudaMemcpyHostToDevice ) );

    HANDLE_ERROR(cudaMemcpy(dev_b,
                            b,
                            N * sizeof(int),
                            cudaMemcpyHostToDevice ) );

    // math.floor division;
    // 512 blocks and 512 threads per block
    add<<<128, 128>>>(dev_a, dev_b, dev_c);

    HANDLE_ERROR(cudaMemcpy(c,
                            dev_c,
                            N * sizeof(int),
                            cudaMemcpyDeviceToHost ) );

    // verify that the GPU did the work we requested
    bool success = true;
    for (int i = 0; i < N; i++) {
        if ((a[i] + b[i]) != c[i]) {
            printf("Error: %d + %d != %d", a[i], b[i], c[i]);
            success = false;
        }
    }

    if (success) {
        printf("We did it!\n");
    }

    auto end = duration<double, milli>(high_resolution_clock::now() - ts_start).count();
    printf("time taken to add %d elements = %f ms\n", N, end);

    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);

    return EXIT_SUCCESS;
}