CUDA 程序的核心思想是异构计算—— CPU（Host）负责控制流与数据搬运，GPU（Device）负责大规模并行计算。 一个 __global__ 函数就是一个 kernel， 它会被成千上万个线程同时执行。

CUDA 的线程组织是一个三级嵌套： 一个 Grid 由多个 Block 构成， 每个 Block 由多个 Thread 构成。 每个 thread 通过内置变量算出自己负责处理数据的哪一个元素。

GPU 的存储层级与 CPU 截然不同—— Global Memory 容量大但延迟高， Shared Memory 延迟低但容量小， 如何编排数据搬运决定了 kernel 性能的上限。

Lab 1 — vector_add.cu

最经典的入门 kernel——两个数组逐元素相加。 目标是理解「全局索引 = blockIdx.x * blockDim.x + threadIdx.x」这一核心公式。

// vector_add.cu — 你的第一个 CUDA kernel
#include <stdio.h>

__global__ void vector_add(float *a, float *b, float *c, int n) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < n) {             // 边界检查：总 thread 数可能大于 n
        c[i] = a[i] + b[i];
    }
}

int main() {
    int n = 1 << 20;        // 1M 个元素
    size_t bytes = n * sizeof(float);

    // 1. Host 内存分配 + 初始化
    float *h_a = (float*)malloc(bytes);
    float *h_b = (float*)malloc(bytes);
    float *h_c = (float*)malloc(bytes);
    for (int i = 0; i < n; i++) {
        h_a[i] = 1.0f;  h_b[i] = 2.0f;
    }

    // 2. Device 内存分配
    float *d_a, *d_b, *d_c;
    cudaMalloc(&d_a, bytes);
    cudaMalloc(&d_b, bytes);
    cudaMalloc(&d_c, bytes);

    // 3. 拷入
    cudaMemcpy(d_a, h_a, bytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, bytes, cudaMemcpyHostToDevice);

    // 4. 启动 Kernel
    int block_size = 256;
    int grid_size  = (n + block_size - 1) / block_size;  // 向上取整
    vector_add<<<grid_size, block_size>>>(d_a, d_b, d_c, n);

    // 5. 拷回 + 验证
    cudaMemcpy(h_c, d_c, bytes, cudaMemcpyDeviceToHost);
    for (int i = 0; i < n; i++) {
        if (h_c[i] != 3.0f) { printf("FAIL at %d\n", i); return 1; }
    }
    printf("PASS: all %d elements correct!\n", n);

    // 6. 清理
    cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
    free(h_a); free(h_b); free(h_c);
    return 0;
}

Lab 2 — naive_gemm.cu

矩阵乘是深度学习的核心计算。 先写一个 naive 版本：每个 thread 算输出矩阵 C 的一个元素， 理解为什么这种写法性能很差（Day 05 再优化）。

// naive_gemm.cu — C[M×N] = A[M×K] × B[K×N]
__global__ void naive_gemm(float *A, float *B, float *C,
                           int M, int N, int K) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < M && col < N) {
        float sum = 0.0f;
        for (int k = 0; k < K; k++) {
            sum += A[row * K + k] * B[k * N + col];
        }
        C[row * N + col] = sum;
    }
}

// 启动配置 (2D grid + 2D block)
int main() {
    int M = 1024, N = 1024, K = 1024;
    // ... malloc, cudaMalloc, cudaMemcpy 省略 ...

    dim3 block(16, 16);              // 16×16 = 256 threads/block
    dim3 grid((N + 15) / 16, (M + 15) / 16);

    naive_gemm<<<grid, block>>>(d_A, d_B, d_C, M, N, K);
    // ... cudaMemcpy D2H, 验证, 计时 ...
}

cuBLAS 基线测量

// 使用 cuBLAS 做同样的矩阵乘
#include <cublas_v2.h>

cublasHandle_t handle;
cublasCreate(&handle);

float alpha = 1.0f, beta = 0.0f;

// cuBLAS 使用列主序 → 参数顺序注意!
cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
             N, M, K,
             &alpha,
             d_B, N,      // B (列主序)
             d_A, K,      // A (列主序)
             &beta,
             d_C, N);     // C (列主序)

// 编译: nvcc -o gemm gemm.cu -lcublas

1. 画出一个 CUDA 程序的五步生命周期（malloc → memcpy → kernel → memcpy → free），标注 Host/Device 的切换点。
2. 给定 n = 10000，blockDim = 256，计算 gridDim 和总 thread 数，指出有多少个 thread 不做有效计算。
3. 写出 2D 情况下的全局索引公式：row = ? , col = ?
4. 解释 Register、Shared Memory、Global Memory 三者的延迟差距，以及 coalesced access 为什么重要。
5. 运行你的 naive GEMM 和 cuBLAS，记录 GFLOPS 差距，分析 naive 版本的主要瓶颈。