CUDA (compute unified device architecture) 是 NVIDIA 所推出的一种并行计算平台和并行计算 api。

CUDA 在并行计算上可以大显神威，因此，我们先要找到一个可并行的问题。一个很简单的可并行问题就是计算无穷级数(infinite series)。圆周率 pi 可以通过一个著名的无穷级数(leibniz formula)进行计算，具体可查看 wiki。我们可以用此公式来逼近圆周率，c 代码如下。

c 实现莱布尼兹级数版 pi

#define LOOP_N 8192000000

double pi() {
    double pi_qv = 1.0;
    int flag = -1;
    for (long i = 1; i < LOOP_N; i++) {
        pi_qv += flag * (1./(2 * i + 1));
        flag = -flag;
    }
    return pi_qv * 4;
}

下面我们编译并运行此程序，并用 time 来测量程序运行时间，在我的机器上差不多需要 8s 吧。

gcc pi_number.c -O3 -opi
time ./pi

cuda 实现莱布尼兹级数版 pi

我们再试试用 cuda 来编写并运行代码

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

#define BLOCK_NUM 32   //块数量
#define THREAD_NUM 256 // 每个块中的线程数
#define LOOP_N BLOCK_NUM * THREAD_NUM * 1000000

__global__ void leib_pi(double* g_sum) {
    const int tid = threadIdx.x;
    const int bid = blockIdx.x;
    double tmp = 0;
    int flag = -1;
    int idx = bid * THREAD_NUM + tid;
    int start = idx * 100000 + 1;
    int end = start + 100000;
    for (int i = start; i < end; i++) {
        tmp += flag * (1./(2 * i + 1));
        flag = -flag;
    }
    g_sum[bid*THREAD_NUM+tid] = tmp;
}

int main(){
    double *h_sum, *g_sum;
    double pi_v = 1;

    // allocate host memory
    h_sum = (double*) malloc(sizeof(double) * BLOCK_NUM * THREAD_NUM);

    // Allocate device memory
    cudaMalloc((void **)&g_sum, sizeof(double) * BLOCK_NUM * THREAD_NUM);

    // Execute kernels
    leib_pi<<<BLOCK_NUM,THREAD_NUM>>>(g_sum);

    // Transfer output from device memory to host
    cudaMemcpy(h_sum, g_sum, sizeof(double)*BLOCK_NUM*THREAD_NUM, cudaMemcpyDeviceToHost);

    for (int i = 0; i < BLOCK_NUM * THREAD_NUM; i++) {
        pi_v += h_sum[i];
    }

    printf("calculate %.10f\n", pi_v*4);

    cudaFree(g_sum);
    free(h_sum);
}

下面我们编译并运行代码，结果显示，程序只用了 0.07s 就跑完了。

nvcc pi_number.cu -o pi
time ./pi

解释 cuda 代码

• 上面的代码中首先我们注意到函数 leib_pi 前面加上了 global 标志，这告诉编译器，此方法是在 gpu 上运行的。而下面的 main 方法则还是在 cpu 上运行。

• cpu 是不能直接访问 gpu 内存的，反过来也不行。在 cuda 中，内存分为 host memory 和 device memory (也称作 gpu memory)，而指向 host memory 的指针被称为 host pointer，指向 device memory 的指针被称为 device memory。

• 在上面的代码中我们首先定义了两个指针 h_sum 指向 host memory ，而 g_sum 指向 gpu memory，然后用 malloc 系统调用给 h_sum 分配内存，用 cudamalloc 给 g_sum 分配 gpu 内存。

• 接下来执行计算

leib_pi<<<BLOCK_NUM,THREAD_NUM>>>(g_sum);

cuda 通过将一组线程组成一个 thread block，并一次 launch 若干个 thread block 的形式来进行并行计算，形成一种网格式的并行结构。在这个例子中，我们 launch 了 32 个 block，每个 block 包含 256 个 thread。你可以理解为有 8192 个线程一起被 launch 出去，并行地开始计算他们各自那部分的内容。每个线程在执行中，可以拿它其对应的 blockId 和 threadId，比如第一个 block 的第一个线程对应的两个 id 就是(0,0)。

const int tid = threadIdx.x;
const int bid = blockIdx.x;

我们将每个线程计算出来的部分结果存入到 gpu memory 中的对应位置。

gpu 计算结束后，返回到 cpu 并调用 cudaMemcpy，将 gpu memory 拷贝到 host memory，然后在 cpu 上完成剩余操作。

用多少 blocks 和 threads?

支持 CUDA 的 gpu 都被分成若干个流处理器(streaming multiprocessors)，每个流处理器有若干个 32 位寄存器，共享内存， a maximum number of active blocks, AND a maximum number of active threads。这些信息可以在 wiki 上找到。比如说你正在使用 2080Ti 的 gpu，你可以看到你的 CC (compute capability) 是 7.5，而对应的 max number of threads per block 是 1024，等等。

注意，number of threads per block 应该始终是 32 的倍数，因为 gpu kernel 总是以 32 个线程一组来启动线程，也就是说，如果你将上面程序中 THREAD_NUM 配置为 50，kernel 会启动 64 个线程(per block)。

此外，我们应当尽可能多的 launch 线程，这就需要考虑你在使用什么 gpu，并根据 gpu 的 CC 来决定具体参数了。

除此之外，还需要考虑寄存器和共享内存这些资源（这些超出了本文的范围）。

PS: 源码见 github。