Your First CUDA Kernel

# Check the CUDA compiler version
$ nvcc --version
nvcc: NVIDIA (R) Cuda compiler driver
Cuda compilation tools, release 12.6, V12.6.77

# Check GPU driver & device info
$ nvidia-smi
+-----------------------------------------------------------------------------+
| NVIDIA-SMI 560.35    Driver Version: 560.35    CUDA Version: 12.6     |
|-------------------------------+----------------------+----------------------+
| GPU  Name        Persistence-M| Bus-Id        Disp.A | Volatile Uncorr. ECC |
| Fan  Temp   Perf  Pwr:Usage/Cap|         Memory-Usage | GPU-Util  Compute M. |
|===============================+======================+======================|
|   0  NVIDIA GeForce ...  Off  | 00000000:01:00.0  On |                  N/A |
+-------------------------------+----------------------+----------------------+

// Runs on GPU, called from CPU
__global__ void myKernel(float* data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) data[idx] *= 2.0f;
}

// Runs on GPU, called from GPU only
__device__ float square(float x) {
    return x * x;
}

// Runs on CPU (default — __host__ is optional)
__host__ void printResult(float* data, int n) {
    for (int i = 0; i < n; i++) printf("%f\n", data[i]);
}

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>

// ── The kernel: each thread computes one element ──
__global__ void vectorAdd(const float* a, const float* b, float* c, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        c[idx] = a[idx] + b[idx];
    }
}

int main(void) {
    int n = 1000000;                     // 1 million elements
    size_t size = n * sizeof(float);

    // ── Step 0: Allocate and initialise host memory ──
    float* h_a = (float*)malloc(size);
    float* h_b = (float*)malloc(size);
    float* h_c = (float*)malloc(size);

    for (int i = 0; i < n; i++) {
        h_a[i] = 1.0f;
        h_b[i] = 2.0f;
    }

    // ── Step 1: Allocate device memory ──
    float* d_a;
    float* d_b;
    float* d_c;
    cudaMalloc((void**)&d_a, size);
    cudaMalloc((void**)&d_b, size);
    cudaMalloc((void**)&d_c, size);

    // ── Step 2: Copy input data Host → Device ──
    cudaMemcpy(d_a, h_a, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, size, cudaMemcpyHostToDevice);

    // ── Step 3: Launch the kernel ──
    int blockSize = 256;                  // threads per block
    int gridSize = (n + blockSize - 1) / blockSize;  // enough blocks to cover n
    vectorAdd<<<gridSize, blockSize>>>(d_a, d_b, d_c, n);

    // ── Step 4: Copy result Device → Host ──
    cudaMemcpy(h_c, d_c, size, cudaMemcpyDeviceToHost);

    // ── Verify the result ──
    int errors = 0;
    for (int i = 0; i < n; i++) {
        if (fabs(h_c[i] - 3.0f) > 1e-5) {
            errors++;
        }
    }
    printf("Vector addition: %s (%d errors)\n",
           errors == 0 ? "PASSED" : "FAILED", errors);

    // ── Step 5: Free device and host memory ──
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
    free(h_a);
    free(h_b);
    free(h_c);

    return 0;
}

// The standard formula: ceiling division
int blockSize = 256;
int gridSize  = (n + blockSize - 1) / blockSize;

// Example: n = 1000
//   gridSize = (1000 + 255) / 256 = 1255 / 256 = 4
//   Total threads = 4 × 256 = 1024 (covers all 1000 elements)
//   Threads 1000–1023 will hit the `if (idx < n)` guard and exit

// 2D grid example for matrix operations:
dim3 block(16, 16);       // 16×16 = 256 threads per block
dim3 grid(
    (width  + block.x - 1) / block.x,
    (height + block.y - 1) / block.y
);
matrixKernel<<<grid, block>>>(d_matrix, width, height);

#ifndef CUDA_CHECK_CUH
#define CUDA_CHECK_CUH

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = (call); \
        if (err != cudaSuccess) { \
            fprintf(stderr, "CUDA error at %s:%d — %s\n", \
                    __FILE__, __LINE__, \
                    cudaGetErrorString(err)); \
            exit(EXIT_FAILURE); \
        } \
    } while (0)

#endif

#include "cuda_check.cuh"

int main(void) {
    float* d_data;
    CUDA_CHECK( cudaMalloc((void**)&d_data, 1024 * sizeof(float)) );

    // Kernel launches don't return cudaError_t, so check after:
    myKernel<<<4, 256>>>(d_data, 1024);
    CUDA_CHECK( cudaGetLastError() );
    CUDA_CHECK( cudaDeviceSynchronize() );

    CUDA_CHECK( cudaMemcpy(h_data, d_data, 1024 * sizeof(float),
                           cudaMemcpyDeviceToHost) );
    CUDA_CHECK( cudaFree(d_data) );

    return 0;
}

# Compile and run the vector addition program
$ nvcc -o vector_add vector_add.cu
$ ./vector_add
Vector addition: PASSED (0 errors)

# Compile device code to relocatable device code
$ nvcc -dc -arch=sm_89 -o kernels.o kernels.cu

# Compile host code
$ nvcc -dc -arch=sm_89 -o main.o main.cu

# Link everything together
$ nvcc -arch=sm_89 -o program kernels.o main.o

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <cuda_runtime.h>

__global__ void vectorMul(const float* a, const float* b, float* c, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        c[idx] = a[idx] * b[idx];
    }
}

int main(void) {
    int n = 1000000;
    size_t size = n * sizeof(float);

    float* h_a = (float*)malloc(size);
    float* h_b = (float*)malloc(size);
    float* h_c = (float*)malloc(size);

    for (int i = 0; i < n; i++) {
        h_a[i] = 1.0f;
        h_b[i] = 2.0f;
    }

    float* d_a; float* d_b; float* d_c;
    cudaMalloc((void**)&d_a, size);
    cudaMalloc((void**)&d_b, size);
    cudaMalloc((void**)&d_c, size);

    cudaMemcpy(d_a, h_a, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, size, cudaMemcpyHostToDevice);

    int blockSize = 256;
    int gridSize = (n + blockSize - 1) / blockSize;
    vectorMul<<<gridSize, blockSize>>>(d_a, d_b, d_c, n);

    cudaMemcpy(h_c, d_c, size, cudaMemcpyDeviceToHost);

    int errors = 0;
    for (int i = 0; i < n; i++) {
        if (fabs(h_c[i] - 2.0f) > 1e-5) errors++;
    }
    printf("Vector multiply: %s (%d errors)\n",
           errors == 0 ? "PASSED" : "FAILED", errors);

    cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
    free(h_a); free(h_b); free(h_c);

    return 0;
}

#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>

__global__ void vectorAdd(const float* a, const float* b, float* c, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        c[idx] = a[idx] + b[idx];
    }
}

int main(void) {
    int n = 10000000;  // 10 million elements
    size_t size = n * sizeof(float);

    float* h_a = (float*)malloc(size);
    float* h_b = (float*)malloc(size);
    float* h_c = (float*)malloc(size);

    for (int i = 0; i < n; i++) {
        h_a[i] = 1.0f;
        h_b[i] = 2.0f;
    }

    float* d_a; float* d_b; float* d_c;
    cudaMalloc((void**)&d_a, size);
    cudaMalloc((void**)&d_b, size);
    cudaMalloc((void**)&d_c, size);

    cudaMemcpy(d_a, h_a, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, size, cudaMemcpyHostToDevice);

    int blockSizes[] = {32, 64, 128, 256, 512, 1024};
    int numTests = 6;

    for (int t = 0; t < numTests; t++) {
        int bs = blockSizes[t];
        int gs = (n + bs - 1) / bs;

        cudaEvent_t start, stop;
        cudaEventCreate(&start);
        cudaEventCreate(&stop);

        cudaEventRecord(start);
        vectorAdd<<<gs, bs>>>(d_a, d_b, d_c, n);
        cudaEventRecord(stop);
        cudaEventSynchronize(stop);

        float ms = 0.0f;
        cudaEventElapsedTime(&ms, start, stop);
        printf("Block size %4d → %8d blocks → %.4f ms\n", bs, gs, ms);

        cudaEventDestroy(start);
        cudaEventDestroy(stop);
    }

    cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
    free(h_a); free(h_b); free(h_c);

    return 0;
}