国产加速器海光DCU&GPGPU深算处理器异构编程实战（中）

#include <cstdio>
#include <cstdlib>
int main()
{
    int N = 256 * 1024 * 1024;    // 数组长度256M
    int size = N * sizeof(float); // 内存空间1GB
    // 分配内存空间
    float *x = (float *) malloc(size);
    float *y = (float *) malloc(size);
    // 初始化a,b,c
    float a = 2.5, b = 2.0, c = 1.0;
    // 循环8次完成
    for (int j = 0; j < 8; j++) {
        // 初始化x
        for (int i = 0; i < N; i++) {
            x[i] = 2.0;
        }
        // 计算y = ax^2 +  bx + c;
        for (int i = 0; i < N; i++) {
            y[i] =  a * x[i] * x[i] + b * x[i] + c;
        }
        // 检查结果
        bool correct = true;
        for (int i = 0; i < N; i++) {
            if (y[i] != 15.0) {
                correct = false;
                break;
            }
        }
        if (correct) {
            printf("Loop %d: Correct results.\n", j+1);
        } else {
            printf("Loop %d: Wrong results.\n", j+1);
        }
    }
    // 释放空间
    free(x);
    free(y);
    return 0;
}

在上面的代码中，为了避免使用过大的内存，我们将从2G个x计算2G个y的流程分8次完成，每次处理256M个。我们首先通过malloc函数分配两个长度为256M的数组x和y（内存空间为2×256M×4B=2GB2×256M×4B=2GB），之后初始化参数a、b、c。之后就进入大循环，在每个大循环中，首先初始化数组x，这在实际问题中可能是从文件中读取，或是以某种方式计算，但简单起见在这里我们让x中的256M个数全都等于2.0，之后计算这些数的二次函数，并将结果保存在数组y中，我们知道y中的所有数都应该等于15.0并据此来检验计算机是否给出了正确的结果。我们将这个大循环执行8次。最后我们释放了数组所占的内存。这是一个串行程序，计算机会使用一个CPU核顺序的执行上述过程。

我们用GNU编译器编译并运行，为了让程序执行得更快，让编译器以O3的级别优化：

g++ cpu_origin.cc -o cpu_origin -O3
./cpu_origin

在我的计算机上，逐行打印出了以下结果，每打印一行凭感觉用时不到一秒：

Loop 1: Correct results.
Loop 2: Correct results.
Loop 3: Correct results.
Loop 4: Correct results.
Loop 5: Correct results.
Loop 6: Correct results.
Loop 7: Correct results.
Loop 8: Correct results.

看来我们确实给了计算机施加了一定的压力。现在我们想知道计算机完成上面的工作具体用了多长时间，为此我们需要编写一个计时器。

题外话：优化代码对于现代编译器来说是非常重要的工作，想体会这一点，请尝试将编译命令中的-O3选项去掉（GNU编译器默认不优化）重新编译并运行。

3.1.2 为程序计时

我们创建一个名为timer.h的头文件，输入以下代码：

#ifndef TIMER_H
#define TIMER_H
#include <sys/time.h>
struct my_timer
{
    timeval ts, te; //起始时刻，终止时刻
    float dt; // 时间间隔,单位毫秒(ms)
    void start(){
        gettimeofday(&ts, NULL);
    }
    void stop(){
        gettimeofday(&te, NULL);
        long int dt_sec  = te.tv_sec - ts.tv_sec;
        long int dt_usec = te.tv_usec - ts.tv_usec;
        dt = dt_sec * 1.0e3 + dt_usec / 1.0e3;
    }
};
#endif

在上面的代码中，我们定义了一个简单的类my_timer，通过成员函数start()记录起始时刻，stop()记录终止时刻并计算时间间隔。代码中利用了sys/time.h头文件中的timeval类，其中包含以秒为单位成员变量的tv_sec和以微秒（μs）为单位的成员变量tv_usec。时间间隔dt是以终止时刻减去起始时刻，并以毫秒（ms）为单位。

现在我们可以用my_timer进行计时了，include头文件timer.h，定义一个计时器timer：

    // 计时器
    my_timer timer;
    double t_init = 0.0;
    double t_calc = 0.0;
    double t_chck = 0.0;

对初始化，计算和结果检验分别计时：

    for (int j = 0; j < 8; j++) {
        // 初始化x
        timer.start();
        for (int i = 0; i < N; i++) {
            x[i] = 2.0;
        }
        timer.stop();
        t_init += timer.dt;
        // 计算y = ax^2 +  bx + c;
        timer.start();
        for (int i = 0; i < N; i++) {
            y[i] =  a * x[i] * x[i] + b * x[i] + c;
        }
        timer.stop();
        t_calc += timer.dt;
        // 检查结果
        timer.start();
        bool correct = true;
        for (int i = 0; i < N; i++) {
            if (y[i] != 15.0) {
                correct = false;
                break;
            }
        }
        if (correct) {
            printf("Loop %d: Correct results.\n", j+1);
        } else {
            printf("Loop %d: Wrong results.\n", j+1);
        }
        timer.stop();
        t_chck += timer.dt;
    }

并在大循环结束后输出这些时间：

    printf("Initialization took %8.3f ms.\n", t_init);
    printf("Calculation    took %8.3f ms.\n", t_calc);
    printf("Validation     took %8.3f ms.\n", t_chck);

编译并运行(cpu_timer.cc)：

g++ cpu_timer.cc -o cpu_timer -O3
./cpu_timer

这次屏幕上显示了程序各部分的计时：

Initialization took  868.988 ms.
Calculation    took 1266.479 ms.
Validation     took 1909.695 ms.

这次我们有了比较精确的计时，程序的大循环部分总计耗时约4秒。这是个不错的结果，毕竟我们让计算机处理了2G个数据，但如果我们还想让程序运行地更快，需要怎样去做呢？

3.1.3 使用多个CPU核

现在CPU早已进入了多核时代，你所用的PC或服务器的配置单上一定会标注你的CPU是多个核的，而在上面的例子中我们只用了一个核心，现在要问的是，多个核能加速上面的程序吗？

观察我们的代码，很快我们发现对于初始化和计算二次函数这两个循环来说，虽然在刚才我们使用一个CPU核心时，循环需要一次接一次的做，但其实每次循环都是独立的，即第i次循环只会负责x[i]和y[i]，而不用处理比如说x[i+1]。在这种情况下，我们用n个CPU核心，每个核负责处理对i循环中的N/n次循环，因为这些核心彼此不必等待，可以让它们同时工作起来，期望速度会提升n倍。

打个比方，好比1千米的路，我们要在路边每隔10米种一颗树，我们让一位工人从头到尾种这100棵树，但这太慢了，所以我们决定多找几位工人来帮助他，比如说再找3位工人。我们给工人们指派任务后他们就开工了，工人A走到10米的位置种第一棵树、工人B走到20米的位置种第二棵树，工人C和D也分别在30和40米的位置开始工作；工人A种完第一棵树之后，走到50米的位置开始种下一颗树，差不多同时，工人B也走到60米的位置继续工作，以此类推。从每个工人的角度，他的同伴正在做什么对他来说并不影响，他们各自开着自己运树苗的车，拿着自己的铁锹和水桶。假设这几位工人的工作效率差不多，我们可以预期他们将以4倍于一个人的速度种完这100棵树。这种分工与使用多CPU核进行并行计算是同一思路，对于我们的程序来说，我们用多个线程（thread，相当于上面比喻中的工人），每个线程使用一个CPU核心（工人自身的劳动力以及工具）同时做计算（种树）。

思考：这种分工方式有什么好处（上面的工作，还可以让每个工人各负责一整段250米）。

提示：如果这段路变成了1.15千米会怎么样。

有了这个思路后，现在的问题是怎样让我们的程序使用多个CPU核，怎样告诉它们去分工。这里我们可以通过OpenMP简单地实现多核并行，OpenMP是一套基于共享内存模型的多线程编程接口，在后面我们还会专门介绍它，现在的使用方法是在初始化和计算两个for循环上面添加一句（以初始化的for循环为例）：

        #pragma omp parallel for
        for (int i = 0; i < N; i++) {
            x[i] = 2.0;
        }

这是一句OpenMP指导语句，用来告诉编译器将下面的for循环并行执行，我们重新编译一下(cpu_omp.cc)，这次我们需要加上选项-fopenmp，让编译器编译时考虑这些指导语句：

g++ cpu_omp.cc -o cpu_omp -O3 -fopenmp

我们配置一个环境变量OMP_NUM_THREADS为2，告诉计算机用两个线程执行OpenMP的并行部分，现在运行，如果这时你有两个空闲的CPU核心，就会让它们都开始工作：

export OMP_NUM_THREADS=2
./cpu_omp

这次的计时是：

Initialization took  429.644 ms.
Calculation    took  640.189 ms.
Validation     took 1915.407 ms.

可以看到初始化和计算的耗时缩短了一半。在我的计算机上有足够的CPU核心，可以进一步设置OMP_NUM_THREADS=4，再次运行：

Initialization took  278.735 ms.
Calculation    took  435.696 ms.
Validation     took 1939.763 ms.

这次耗时是单个核心的三分之一，虽然没达到4倍的加速，但看来多个CPU核心并行计算确实有用。这时你可能想，如果我有更多的核心（设想一下成千上万个）会不会更快呢？

思考：现在花在检查结果的时间最长了，请尝试对它也使用并行加速。

3.1.4 为DCU编写代码

DCU是专为并行计算设计的，此刻你可以认为它是有几千个核心的处理器。但是目前的程序不会自动使用DCU，我们需要对程序做一些修改。在本章，我们介绍使用HIP编程接口为DCU编程，我们也可以通过OpenMP来使用DCU，在后文有单独章节介绍。

现在开工，首先DCU有自己的内存，和CPU使用的内存（下文称主存）是不一样的，所以需要对DCU的内存进行管理，为此我们在代码中添加以下内容：

    float *dx, *dy;
    hipMalloc(&dx, size);
    hipMalloc(&dy, size);

dx和dy分别是数组x和y在DCU内存中的对应，我们给它们分配主存上同样大小的内存空间，但这一次我们使用hipMalloc()而非malloc()。

为了给DCU编程，我们需要在代码中使用一些DCU专用的函数（专用名词：运行时应用编程接口，runtime API，得名于这些编程接口的作用是与运行时环境进行交互），它们以hip前缀开头，包含在一个名为hip_runtime.h的头文件中，hipMalloc是其中之一，下面我们还会用到其他的。记得在文件开头加一行：

#include <hip/hip_runtime.h>

现在到初始化数组x的for循环了，但是DCU是无法直接执行这个for循环的，我们得为它专门写一个函数：

__global__ void init(float *x, int N)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) {
        x[i] = 2.0;
    }
}

这个函数在DCU编程中被称为核函数，它前面有一个修饰符__global__，表示这是可以被CPU调用，并在DCU上执行的函数。

观察一下这个函数，你会发现上面的for循环不见了，只保留了循环终止的判断i<N，i也不再从0开始递增，变成了用几个没有定义的变量计算（blockIdx.x、blockDim.x以及threadIdx.x）得到。回忆一下，刚才我们说过使用多核CPU做并行计算时，是通过多个线程，每个线程使用一个CPU核心的方式同时计算。但我们编程的时候并没有遇到线程这个概念，只是在运行之前设置了一个OMP_NUM_THREADS的线程数，这是因为OpenMP让编译器帮助你做了很多工作，现在我们得自己来做这些。

DCU上能用更多的线程，这些线程是分块来管理的，称为线程块，blockIdx.x是线程块的索引（即编号），blockDim.x是线程块的大小，threadIdx.x是每个线程块内线程的索引，这些都是DCU的内置变量，因此你不用定义就可以使用（看到后面的.x你应该就能猜到线程块可以是多维的，我们会在后面提到），DCU在执行核函数之前，会为每个线程唯一确定这些变量。这时，我们看第一个线程块（blockIdx.x=0）的线程，对于每个线程来说都有一个私有的i，就等于线程索引threadIdx.x，比如索引是0的线程，它的i是0；索引是10的线程，它的i是10。假设线程块的大小是256（blockDim.x=256），对于第二个线程块（blockIdx.x=1）中索引为0的线程，i是256；索引是15的线程，它的i是271，以此类推。每个线程负责把数组x的一个元素x[i]写成2.0。所以方式还是一样，每个线程使用一个DCU核心完成这N次循环中的一次，一大堆线程同时开始工作，直到把整个数组初始化完。当然我们保留了判断（i<N），以此避免有线程对超过数组边界的内存写入数据。

现在我们调用这个核函数，这也需要特殊的写法：

        int blockSize = 256;
        int numBlocks = (N + blockSize - 1) / blockSize;
        // 初始化x
        timer.start();
        init<<<numBlocks, blockSize>>>(dx,N);
        timer.stop();
        t_init += timer.dt;

我们按上面的想法，把线程块blockSize定为256，线程块的个数numBlocks通过代码中的公式确定，用以保证在N无法被blockSize整除时总的线程数大于N。调用函数时init时，在函数后加上由三对尖括号扩起来的线程块个数和线程块大小，最后再跟上函数的参数。这是在程序中调用在DCU上运行的核函数的固定写法。

不必担心总线程数256M是否过大，我们在下一节会介绍DCU启动线程的方式。但注意线程块大小（blockSize）是不能大于1024的。

按照相同的方式，我们再写一个负责计算二次函数的核函数，将另一个for循环替换掉：

__global__ void calculate(float *x, float *y, float a,
                          float b, float c, int N)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) {
        y[i] = a * x[i] * x[i] + b * x[i] + c;
    }
}

并调用：

        calculate<<<numBlocks, blockSize>>>(dx,dy,a,b,c,N);

至此，我们可以检查结果了，我们暂时还是让CPU串行的来做，但问题是刚才提到了，DCU内存和计算机主存是不同的，虽然我们在DCU有了结果，但CPU却还不知道，因此我们需要将DCU内存上的数组dy的内容传递到主存上的数组y，通过以下运行时编程接口：

        hipMemcpy(y, dy, size, hipMemcpyDeviceToHost);

我们把它包含在检查结果的计时范围之内（毕竟如果不需要检查结果也就用不着传输了）。再检视一下整个代码，我们还需要在程序最后释放掉DCU内存中的dx和dy（虽然在这里程序马上就结束了，释放并不重要，但用好的编程习惯规范地编程会让我们少犯很多错误）:

    hipFree(dx);
    hipFree(dy);

另外我们还发现其实有了数组dx之后，x也不再需要了，我们将对x的malloc和free都删去，大功告成！

方便起见，将DCU版代码完整展示如下（dcu.cc）：

#include <cstdio>
#include <cstdlib>
#include "timer.h"
#include <hip/hip_runtime.h>
__global__ void init(float *x, int N)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) {
        x[i] = 2.0;
    }
}

__global__ void calculate(float *x, float *y, float a,
                          float b, float c, int N)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) {
        y[i] = a * x[i] * x[i] + b * x[i] + c;
    }
}

int main()
{
    int N = 256 * 1024 * 1024;    // 数组长度256M
    int size = N * sizeof(float); // 内存空间1GB
    // 分配内存空间
    float *y = (float *) malloc(size);
    float *dx, *dy;
    hipMalloc(&dx, size);
    hipMalloc(&dy, size);
    // 计时器
    my_timer timer;
    double t_init = 0.0;
    double t_calc = 0.0;
    double t_chck = 0.0;
    // 初始化a,b,c
    float a = 2.5, b = 2.0, c = 1.0;
    // 循环8次完成
    for (int j = 0; j < 8; j++) {
        // 初始化x
        int blockSize = 256;
        int numBlocks = (N + blockSize - 1) / blockSize;
        timer.start();
        init<<<numBlocks, blockSize>>>(dx,N);
        timer.stop();
        t_init += timer.dt;
        // 计算y = ax^2 +  bx + c;
        timer.start();
        calculate<<<numBlocks, blockSize>>>(dx,dy,a,b,c,N);
        timer.stop();
        t_calc += timer.dt;
        // 检查结果
        timer.start();
        bool correct = true;
        hipMemcpy(y, dy, size, hipMemcpyDeviceToHost);
        for (int i = 0; i < N; i++) {
            if (y[i] != 15.0) {
                correct = false;
                break;
            }
        }
        if (correct) {
            printf("Loop %d: Correct results.\n", j+1);
        } else {
            printf("Loop %d: Wrong results.\n", j+1);
        }
        timer.stop();
        t_chck += timer.dt;
    }
    printf("Initialization took %8.3f ms.\n", t_init);
    printf("Calculation    took %8.3f ms.\n", t_calc);
    printf("Validation     took %8.3f ms.\n", t_chck);
    // 释放空间
    hipFree(dx);
    hipFree(dy);
    free(y);
    return 0;
}

3.1.5 编译运行DCU程序

现在该编译我们的DCU版代码了（dcu.cc），GNU编译器是无法完成的，我们需要DTK软件栈中的hipcc工具，这次不需要-O3因为对hipcc这是默认的：

hipcc dcu.cc -o dcu

我们运行程序，得到以下结果：

Initialization took    1.571 ms.
Calculation    took    0.070 ms.
Validation     took 3387.858 ms.

这个结果看起来很棒，我们的计时器显示除了验证结果的耗时增长了不少之外，初始化和计算的部分几乎都不耗时了。但我必须告诉你，这里我们的计时器显示的结果是有问题的。我们用DTK软件栈中的性能分析工具来看看真实的性能，使用命令：

hipprof ./dcu

这次屏幕的最后多了两张表格（图3-1）表格里记录了我们调用的HIP编程接口的耗时以及核函数的耗时，结果是以纳秒（ns）为单位的，我们换算一下，数据传输hipMemcpy总耗时1.8秒，核函数init和calculate耗时都是约34毫秒。根据这些真实的时间，单从核函数和for循环对比，相比单核CPU分别提速了大约25倍和37倍。

图3-1 hipprof统计的耗时

你可能觉得传输的时间也对不上，但这是因为编译器不同造成的，hipcc同样可以编译CPU版的代码（使用clang编译器），编译后运行的检查结果耗时大约1.6秒，DCU版多了传输时间1.8秒，总耗时约3.4秒。

3.1.6 核函数执行是异步的

现在你可能会问，我们的计时器为什么不准呢？这是因为DCU执行核函数与CPU执行其他部分是异步的。回到我们的代码：

        int blockSize = 256;
        int numBlocks = (N + blockSize - 1) / blockSize;
        // 初始化x
        timer.start();
        init<<<numBlocks, blockSize>>>(dx,N);
        timer.stop();
        t_init += timer.dt;

上面这些代码其实都是由单个CPU核串行执行的，首先设置两个变量存储线程块的大小和线程块的个数，让计时器记下初始时刻。随后调用核函数init，CPU会执行一系列操作，把init函数的相关信息打包发送给DCU。此时CPU就不再管init函数了，完全由DCU独立执行，CPU会继续往下，让计时器记录终止时刻并计算时间间隔。

说的更清楚一点，如果CPU执行和DCU核函数执行是同步的，CPU会等待DCU核函数执行完毕后，再继续执行timer.stop()；但实际情况是异步的，CPU把任务交给DCU之后，不待DCU运行完毕就直接继续执行timer.stop()了。这就解释了为什么我们的计时器不准，因为它实际记录的是CPU把核函数init打包交给DCU的时间，这段耗时相比DCU执行核函数是很短的。

为了让计时器记录准确的核函数执行时间，我们需要在timer.stop()之前让CPU与DCU同步，通过调用HIP运行时编程接口hipDeviceSynchronize()：

        timer.start();
        init<<<numBlocks, blockSize>>>(dx,N);
        hipDeviceSynchronize();
        timer.stop();

当CPU执行到这个接口时会等待DCU上的现有的任务执行完毕后再继续。对于核函数calculate我们同样加上hipDeviceSynchronize()，用hipcc重新编译并运行，这次的计时是正确的：

Initialization took   35.847 ms.
Calculation    took   34.363 ms.
Validation     took 3321.971 ms.

在介绍了上面的异步问题之后，你可能会担心，比如会不会DCU内存没有申请完就开始执行核函数；核函数init没有执行完，数组dx还有一部分没有初始化就开始执行核函数calculate；calculate还没执行完就开始执行dy到y的拷贝；dy中的内容还没有完全拷贝到y中就开始检查结果。

对于目前的代码，不需要有这些顾虑，因为hipMalloc、hipFree、hipMemcpy这几个运行时编程接口都是同步的，它们会等待DCU上现有的任务全部执行完再执行，并且要等它们自身的操作（分配/释放内存、拷贝）执行完成后才允许CPU继续执行后面的代码。以hipMemcpy为例，即使我们没有改进我们的计时器（没有在calculate后面加hipDeviceSynchronize接口），当程序执行到hipMemcpy时也会等待calculate执行完再执行从DCU上的dy到主存中的y的拷贝，并且会在这个拷贝执行完成后才进入检验y的每个元素是否等于15.0的for循环（因此我们的计时器对检查结果部分的计时基本正确）。此外，DCU会顺序执行核函数，先执行完init之后再执行calculate。

但需要指出，DCU编程有很大的自由度，我们可以让拷贝异步的执行，也可以让多个核函数彼此异步的执行，这些是3.4节中的内容。

思考：请对比核函数calculate后加hipDeviceSynchronize对hipMemcpy耗时的影响并思考原因。

3.1.7 小节

在本节中我们编写了一个包含大量二次函数计算的程序，并使用单个CPU核运行。我们尝试用并行计算对程序进行加速，首先通过OpenMP使用了多个CPU核，之后修改了代码使得程序能使用DCU。在这个过程中，介绍了并行计算的思路，初步尝试了OpenMP编程和DCU编程，并测试了加速效果。

在传统计算机上，一个程序的所有指令都是由单个CPU核心执行的，但在现代计算机上，常常会根据指令进行专门化的分工，将原本由CPU处理的部分工作卸载（offload）给其他硬件或硬件单元。这些专门设计的硬件或硬件单元在处理这部分工作时相比CPU往往有更快的性能，卸载这部分工作之后，还能使CPU继续处理其他工作。对于我们上面的例子而言，通过OpenMP使用多个CPU核其实就是把单个CPU核的两个for循环部分卸载给其他CPU核心，而使用DCU则是直接将CPU的对两个for循环的计算全部卸载给DCU。另一个CPU工作卸载的例子是CPU与DCU之间的数据传输，这类工作通常也不会使用CPU执行，而是卸载给DCU上被称为DMA（Direct Memory Access，直接内存访问）的专用硬件单元。

DCU是专为并行计算设计的硬件。如本节中的DCU版代码所示，要使用DCU，或者说让CPU将程序中适合并行计算的部分卸载给DCU来执行，有两个方面的工作要做，其一是使用以hip为前缀的HIP编程接口，其二就是编写核函数让DCU执行计算。

3.2 DCU编程概述

在上一节中，我们对并行计算做了简要的介绍，并初步尝试了DCU编程，你现在应该对一个DCU程序的代码样式有所了解了。为了能更有效的使用DCU做并行计算，需要对DCU编程做一个系统性的概述。首先介绍DCU编程模型，这里的编程模型是指对计算机系统本身以及计算机如何执行程序的抽象，特别是对DCU硬件以及DCU如何执行核函数；之后简要介绍从编程模型到DCU硬件的映射关系。

3.2.1 异构编程模型

我们的程序主体是用CPU执行的，在程序的某些阶段，CPU会把一些工作交给（卸载）DCU执行，因而我们的程序是用两种处理器执行的，这称为异构计算（heterogeneous computing）。CPU和DCU的关系是主从关系，因此在后文我们会用一个抽象的概念主机（host）来指代CPU，设备（device）来指代DCU。主机使用主存，设备有单独的设备内存。在一台计算机上可以有多个主机与多个设备相互连接（从硬件讲是通过PCIe总线），可以多个主机使用一个设备，也可以一个主机使用多个设备。

主机通过调用运行时编程接口与运行时环境交互，从而操作设备，实现诸如在多个设备时选择使用的设备；为存储在设备内存上的数据分配/释放/清空设备内存；发起主机与设备、单个设备上、设备与设备间的数据拷贝等等。

在前面的例子中已经使用了三个HIP运行时编程接口：hipMalloc、hipMemcpy和hipFree，它们从名称和参数上都尽量做的与C语言标准库函数malloc、memcpy和free相似。这是为了让熟悉C语言的开发者快速上手。后面我们还会用到更多的接口以实现更复杂的功能。

还可以在主机端调用核函数，让设备执行核函数进行并行计算。对于设备端执行的核函数，我们需要单独编写代码，并在核函数前加__global__修饰符以区分，核函数中还可以调用在设备端执行的函数，这些函数前需要加__device__修饰符。主机调用核函数后会继续执行后面的主机端代码，而设备端则会根据配置启动多个线程同时执行核函数。

仍以我们上节中的核函数init自身实现及其调用为例，

// 实现
__global__ void init(float *x, int N)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N) {
        x[i] = 2.0;
    }
}

// 在main中的调用
        int blockSize = 256;
        int numBlocks = (N + blockSize - 1) / blockSize;
        init<<<numBlocks, blockSize>>>(dx,N);

我们在main函数中调用核函数init时，设置了线程块的个数为1M，每个线程块中的线程数为256，它们合起来被称为核函数启动参数（Kernel launch configuration），我们通过将它们置于核函数后的三个尖括号之中，让设备以这样的线程组织结构，启动256M个线程，每个线程都执行一遍init函数，换言之将init函数执行了256M次。

设想一下如果把核函数中的x[i]改成x[0]，则调用核函数会对x[0]赋值256M次，每次都赋值为2.0，这当然是冗余且无意义的；这里的关键是每个设备线程有不同的i，每个线程执行一遍核函数会把数组的不同元素x[i]赋值为2.0，核函数执行256M次就会把长度为256M的数组的所有元素初始化为2.0。即利用对于每个设备线程各不相同的内置变量（blockIdx、blockDim、threadIdx等）作为索引，控制所有线程同时独立地做不同的事。

设备的线程具有层次结构，一组线程可以构成线程块，一组线程块可以构成线程网格，它们由核函数启动参数配置，设备会根据启动参数在核函数启动时给每个线程设置内置变量，包括线程网格包含的线程块数目（gridDim）、线程块在线程网格中的编号/索引（blockIdx）、线程块包含的线程数（blockDim）、线程在线程块中的编号/索引（threadIdx）。这些维度和索引都是三维的，以适应程序中的二维、三维问题。它们有专门的数据类型dim3，可以用1到3个参数来初始化（未指定的维度默认为1），比如我们想设置16×16的二维线程块，可以声明如下：

dim3 dimBlock(16,16);

如果想修改，可以：

dimBlock = dim3(256,1,1);

图3-2 线程层次结构

每个线程块最多包含1024个线程，一个线程块内的所有线程可以同步。不同线程块之间则是彼此独立的。

对应线程层次结构，设备内存也具有层次结构，每个线程有私有的局部内存；同一个线程块中的线程有共享内存（变量前加修饰符__shared__）；所有线程都可以访问同一全局内存。

图3-3 线程层次结构

3.2.2 从模型到硬件

DCU芯片上有64个独立的计算单元，每个计算单元有4个独立的SIMD（Single Instruction，Multiple Data，单指令、多数据）单元，每个SIMD单元又有16个向量ALU（算数逻辑单元），通俗的讲，可以认为一块DCU芯片上有4096个核心，每个向量ALU就对应一个线程。

虽然线程从概念上是独立执行的，但从DCU硬件来说是64个线程一组并行执行的，即一条指令发布后，一个SIMD单元中的16个ALU执行4轮完成，这64个线程被称为线程束（warp）；一个SIMD单元有10个线程束的指令缓冲，所以最多可以同一时间段内并发执行10个线程束（受其他硬件资源的限制）。这样算64个计算单元×4个SIMD单元×10线程束×64线程=163840，这是DCU硬件最大同时支持的线程数。此外，这10个线程束可以来自不同的线程块或者核函数，意味着我们的DCU支持多个核函数并发（在3.4节我们会介绍如何编程实现）。

每个计算单元配有一个标量ALU，负责线程束的流程控制，这里的标量是指ALU对整个线程束处理一个数，而向量ALU则对每个线程处理一个数。举例来说，在流程控制时，每个线程会根据判断条件得到一位（bit）以进行分支或循环，整个线程束是64位，标量ALU负责处理这个64位数。

计算单元上有高速的寄存器文件，每个SIMD单元有64KB的向量寄存器配额供向量ALU使用，负责存储我们在核函数中每个线程独有的局部变量。每个计算单元有64KB的局域数据共享（Local Data Share），负责存储核函数中以__shared__修饰符声明的变量，以实现同一线程块或同一线程束内的数据共享。此外，同现代CPU一样，DCU上还有其他复杂的多级缓存结构以实现对片外存储器的高速访问，但程序开发者是无法直接控制这些缓存的，只能通过修改算法提高访存的连续性以提升缓存命中。上面所介绍的这些存储器件和计算单元一起（当然还包括其他的元件）封装在一块芯片内，构成片上存储。

DCU一般配有16GB或32GB的片外存储，对应我们上面说的设备端全局内存。计算单元请求片外内存上的数据是非常耗时的过程，通常需要几百个时钟周期，因此我们需要充分利用内存层次结构以提高程序访存性能，这也是后文重点要介绍的。

下面我们从硬件角度介绍核函数的执行。编译器在编译时会确定每个核函数的线程束需要多少硬件资源（比如说多少个寄存器，多大的共享内存），在程序运行时，DCU上的调度元件会根据所有计算单元的运行状态将核函数的线程束轮询的发送到空闲的且满足资源需求的计算单元上。这种轮询机制使得即使核函数启动参数中设置的线程数很多（如上节中的256M），仍然能够执行，只是不可以同时并发。同一个线程块中的线程束肯定会发送到同一计算单元上，以利用计算单元上的局域数据共享和专门的同步指令，DCU允许最多16个线程束构成一个线程块，这就是前面我们反复提到单个线程块中的线程最多是1024个的原因。单个线程块中的线程有唯一的ID，按以下公式计算：thread ID = threadIdx.x + threadIdx.y ×blockDim.x + threadIdx.z ×blockDim.x ×blockDim.ythread ID = threadIdx.x + threadIdx.y ×blockDim.x + threadIdx.z ×blockDim.x ×blockDim.y线程按线程ID以64个为一组划分在不同的线程束中，比如对于一维的（128，1，1）线程块，threadIdx.x为0到63的线程在同一线程束，64到127的在另一线程束；对于二维的（16，16，1）线程块，按threadIdx.y从0到3，4到7，8到11，12到15划分在4个线程束上。了解这个划分对访存是有用的，比如说我们知道在这样的二维线程块中，threadIdx.x连续的线程是硬件上连续的，threadIdx.y连续的线程则不然，在后面我们会介绍合并访存的概念，需要连续的线程访问连续的地址；此外还有助于针对线程束的函数的使用。顺便再强调一次，在编程模型中允许有小于线程束的线程块，但在这种情况下DCU依然会启动一个线程束来执行。

3.2.3 小节

在本节中我们简要介绍了DCU编程模型，模型的重点在于主机和设备的主从关系，设备线程的层次结构和设备内存的层次结构。了解DCU编程模型有助于我们构思设计DCU程序，也便于我们后文的介绍。之后我们深入到DCU的硬件，对DCU的架构做了简要的介绍，并说明了硬件单元与编程模型中抽象概念的对应。这方面我们没有过于深入，只需对线程束、线程ID这些源自硬件的概念有所了解即可，在后文中我们会看到这些概念与DCU程序性能的关系。其他硬件方面的信息我们也将在后文中根据需要适时介绍。

3.2 并行规约

在对DCU编程有个概述之后，我们继续使用DCU做并行计算，这次打算计算一个复杂的定积分，以此为例进一步探究DCU编程。我们在网络上搜索积分表，找到如下积分公式：

不要在意上面的公式，我们只想给计算机一些复杂的计算任务。

我们令a=2a=2，b=1b=1，化简一下：

我们取积分限为ππ的整数倍，这样sin都没了，容易判断结果是否正确。我们先计算一个从0到40000ππ的积分，编写代码如下（integral目录下的cpu.cc）：

#include <cstdio>
#include "timer.h"
#include <cmath>
// 定义精度
#ifdef DOUBLE_PREC
typedef double real;
#define SIN(x) sin(x)
#define COS(x) cos(x)
#else
typedef float real;
#define SIN(x) sinf(x)
#define COS(x) cosf(x)
#endif

// 用积分公式计算
double F(double x)
{
    double F = x/4 + sin(2*x)/16 - sin(4*x)/16 -sin(6*x)/48;
    return F;
}

int main()
{
    double lower_bound = 0.0;
    double upper_bound = 40000 * acos(-1); //acos(-1) = pi
    int N = 1e8;
    real dx = (upper_bound - lower_bound) / N;
    double Integral = 0.0;
    my_timer timer;
    timer.start();
    // 数值计算 \int sin^2(2x)*cos^2(2x) dx
    for (int i = 0; i < N; i++) {
        real x = (i + 0.5) * dx;
        // 减少乘法次数
        real sin_2x = SIN(2 * x);
        sin_2x *= sin_2x;
        real cos_x = COS(x);
        cos_x *= cos_x;
        Integral += sin_2x * cos_x * dx;
    }
    timer.stop();
    printf("Numeric result is: %.11f\n", Integral);
    printf("Result form formula is: %.11f\n",
           F(upper_bound) - F(lower_bound));
    printf("Elapsed time: %8.6f ms\n",timer.dt);
    return 0;
}

在上面的代码中，0到40000ππ的区间被分为了10的8次方（一亿）段，取每段的中点x计算sin2(2x)cos2(bx)dxsin2(2x)cos2(bx)dx，最后累加到变量Integral上，并通过一点技巧，使计算sin、cos和乘法的次数减少了一些。我们对积分限和变量Integral采用双精度浮点数存储，利用积分公式计算的F(x)函数使用双精度，其他变量以及sin、cos函数的计算则可以通过条件编译选择不同的精度。

单精度版本编译运行如下：

g++ cpu.cc -o cpu -O3
./cpu
Numeric result is: 31415.92652882586
Result form formula is: 31415.92653589793
Elapsed time: 2773.030029 ms

双精度版本编译运行如下：

g++ cpu.cc -o cpu_dp -O3 -DDOUBLE_PREC
./cpu_dp
Numeric result is: 31415.92653618794
Result form formula is: 31415.92653589793
Elapsed time: 8468.580078 ms

尝试：

（1）使用单精度存储变量Integral并测试结果。

（2）编写DCU版本的代码（先不要看下文）。

3.3.1 原子操作

我们编写DCU版代码（dcu_error.cc）：

#include <cstdio>
#include "timer.h"
#include <cmath>
#include <hip/hip_runtime.h>
// 定义精度
#ifdef DOUBLE_PREC
typedef double real;
#define SIN(x) sin(x)
#define COS(x) cos(x)
#else
typedef float real;
#define SIN(x) sinf(x)
#define COS(x) cosf(x)
#endif

// 用积分公式计算
double F(double x)
{
    double F = x/4 + sin(2*x)/16 - sin(4*x)/16 -sin(6*x)/48;
    return F;
}

// DCU 核函数
__global__
void integral_kernel(double* Integral, double lb, int N, real dx)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int tid = threadIdx.x;
    if (i < N) {
        real x = lb + (i + 0.5) * dx;
        real sin_2x = SIN(2 * x);
        sin_2x *= sin_2x;
        real cos_x = COS(x);
        cos_x *= cos_x;
        *Integral += sin_2x * cos_x * dx;
    }
}

int main()
{
    double lower_bound = 0.0;
    double upper_bound = 40000 * acos(-1); //acos(-1) = pi
    int N = 1e8;
    real dx = (upper_bound - lower_bound) / N;
    double Integral = 0.0;
    double* d_Integral;
    hipMalloc(&d_Integral,sizeof(double));
    // 核函数启动配置
    int blockSize = 256;
    int numBlocks = (N + blockSize - 1) / blockSize;
    // 开始计时
    my_timer timer;
    timer.start();
    hipMemcpy(d_Integral,&Integral,sizeof(double),
              hipMemcpyHostToDevice);
    integral_kernel<<<numBlocks, blockSize>>>
                   (d_Integral,lower_bound,N,dx);
    hipMemcpy(&Integral,d_Integral,sizeof(double),
              hipMemcpyDeviceToHost);
    timer.stop();
    printf("Numeric result is: %.11f\n", Integral);
    printf("Result form formula is: %.11f\n",
           F(upper_bound) - F(lower_bound));
    printf("Elapsed time: %8.6f ms\n",timer.dt);
    hipFree(d_Integral);
    return 0;
}

在上面的代码中，我们在设备上申请内存用来保存积分结果，这块内存通过指针d_Integral访问，在主机端将Integral初始化为0并拷贝到设备上，将循环改为设备端的核函数，以256为线程块大小执行核函数，之后将结果拷贝回Integral中并释放设备内存。

编译运行如下：

hipcc dcu_error.cc -o dcu_error
./dcu_error
Numeric result is: 0.00072969924
Result form formula is: 31415.92653589793
Elapsed time: 4.727000

注意，DCU计算的结果是错误的！仔细检查一下，看起来代码写的并没有什么问题。

造成错误的原因其实是并行计算中的一个常见场景，让多个线程并行做计算，但最后需要把所有结果归并为一个值，这类问题被称为并行规约问题。在上面的代码中，多个线程会同时从d_Integral（所保存的地址）中取值，把本线程的结果加到这个值上再存回d_Integral中。以两个线程为例，假设此时d_Integral的值还是0，线程1算出结果dF1dF1，线程2算出结果dF2dF2，线程1从d_Integral处取值0，0加上dF1dF1还是dF1dF1，此时线程1把dF1dF1写回d_Integral；如果线程1还没完成写回，线程2就从d_Integral处取值，仍然得到0，0加dF2dF2等于dF2dF2，线程将dF2dF2写回d_Integral。这种情况下如果线程1先完成写回，d_Integral最终就是dF1dF1，如果线程2先完成写回，结果就是dF2dF2，但无论如何不是正确的结果dF1+dF2dF1+dF2。

解决方法是使用原子操作，所谓“原子”的意思是上面叙述的一个线程对d_Integral的读-加-写回这三个动作是不可分的，借用了“原子”这个词不可再分的原意。当一个线程执行这个操作时，必须得把读-加-写回这三个操作全部做完，在此期间别的线程不能执行此操作。

修改一下代码，把核函数中的

*Integral += sin_2x * cos_x * dx;

修改为

atomicAdd(Integral,sin_2x * cos_x * dx);

并把分段数N从1e8改为1e6，这是为了使核函数执行不至于太慢。编译并运行：

hipcc dcu_error.cc -o dcu_error.cc
./dcu_error
Numeric result is: 31415.92749216406
Result form formula is: 31415.92653589793
Elapsed time: 9127.658203 ms

从正确性的角度看，虽然精度低了一些，但结果好了很多；从性能的角度看，惨不忍睹，在计算量少了100倍的情况下还用了9.1秒。我们还需继续改进。

3.3.2 利用共享内存优化原子加

在我们的程序中d_Integral保存的是设备端全局内存中的地址。在3.2.2节，我们介绍了线程访问DCU片外存储是非常耗时的过程，此外对于原子操作只有对整型的硬件支持，对浮点数的原子操作则是通过软件实现的，这就是上面让所有线程对d_Integral做原子加性能很差的原因。

在3.2.2节中介绍过，DCU编程的一个重点就是合理利用内存层次结构，在全局内存与局部内存之间，还有一个线程块内所有线程公用的共享内存（对应硬件的LDS），线程对这块内存的访问速度相比全局内存要快很多，下面我们就利用它来改善原子加的性能。

将核函数修改如下：

__global__
void integral_kernel(double* Integral, double lb, int N, real dx)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int tid = threadIdx.x;
    __shared__ double tmp_sh;
    if (tid == 0) tmp_sh = 0.0;
    __syncthreads();
    if (i < N) {
        real x = lb + (i + 0.5) * dx;
        real sin_2x = SIN(2 * x);
        sin_2x *= sin_2x;
        real cos_x = COS(x);
        cos_x *= cos_x;
        atomicAdd(&tmp_sh,sin_2x * cos_x * dx);
    }
    __syncthreads();
    if (tid == 0 ) atomicAdd(Integral,tmp_sh);
}

在这段代码中，首先定义了线程ID的变量tid，这是一种常用的命名，它的值在一维的情况下就等于threadIdx.x。此外多了共享内存中变量tmp_sh的声明（通过__shared__修饰符），我们让线程ID为0的线程（tid == 0）将其初始化为0。后面的__syncthreads()是一个内置函数，作用是令线程块同步，即等待线程ID为0的线程初始化tmp_sh之后，所有线程才能继续往下执行。在最后的原子加阶段，我们让线程块内的所有线程向共享内存做原子加，等完成之后，再由线程ID为0的线程对全局内存中的d_Integral做一次原子加即可。

这里需要多解释几句，前面介绍过，DCU是按64个线程一组的线程束来调度的，当遇到像if (tid == 0)这样的条件判断时，标量ALU会把线程束的所有线程ID与0做比较得出一个64位的掩码，在这里代表线程ID是0的那一位是1，其他都是0，当线程束执行对共享内存中的变量tmp_sh赋值时，通过这个掩码使得只有线程ID为0的线程实际执行。当然在这个过程中，线程束内的其他线程并不能继续向下执行，而需要等线程ID为0的线程执行完。

使用__syncthreads()进行线程块同步对于线程数大于64的线程块是必要的，因为必须将做累加的变量tmp_sh初始化为0，因为线程块中线程束的调度顺序是不确定的，如果不包含线程ID为0的线程束先执行，在不加线程块同步的情况下很可能这些线程束将计算结果加到未初始化的tmp_sh中，之后等包含线程ID为0的线程束执行时，又会把tmp_sh置为0，最后的结果当然是不对的。同理对于最后的原子加阶段，如果线程块内还有线程未完成对共享内存中tmp_sh的原子加，线程0就把tmp_sh的结果加到全局内存d_Integral中去，同样会造成错误。

做了以上修改后，可以把N改回1e8了，编译并运行（shared.cc）：

hipcc shared.cc -o shared
./shared
Numeric result is: 31415.92577126694
Result form formula is: 31415.92653589793
Elapsed time: 185.567001 ms

这样，通过利用共享内存，我们的性能有了大幅的提高。

3.3.3 利用线程束洗牌函数优化原子加

到了这一步，每个线程块对全局内存做一次原子加，如果我们想继续减少全局内存的原子加次数就需要减少线程块的个数，对此，一个自然的想法是增加线程块的大小。将blockSize修改为1024（对于使用大于256的线程块启动的核函数，需要在声明前加__launch_bounds__，这与编译器优化有关，详见下面的核函数代码），编译并运行，得到：

Numeric result is: 31415.92577126688
Result form formula is: 31415.92653589793
Elapsed time: 616.575989 ms

结果并未如我们所愿，耗时反而更长了。我们猜测这可能与对共享内存的原子加增多有关，当blockSize为256时，每个线程块的256个线程向共享内存的同一地址做原子加，blockSize增大为1024时，则是1024个向同一个地址。是否有办法优化对共享内存的原子加呢？

实现的方法有很多种，在此我们通过一组被称为线程束洗牌函数（warp shuffle）的内置函数来实现。它们的作用是让线程获取到线程束内其他线程的局部变量的值，从硬件层面讲仍然是利用共享内存实现的。我们将核函数修改如下：

__global__ __launch_bounds__(1024)
void integral_kernel(double* Integral, double lb, int N, real dx)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int tid = threadIdx.x;
    double tmp = 0.0;
    __shared__ double tmp_sh;
    if (tid == 0) tmp_sh = 0.0;
    __syncthreads();
    if (i < N) {
        real x = lb + (i + 0.5) * dx;
        real sin_2x = SIN(2 * x);
        sin_2x *= sin_2x;
        real cos_x = COS(x);
        cos_x *= cos_x;
        tmp += sin_2x * cos_x * dx;
        int delta = 1;
        for (int j = 0; j < 6; j++) {
            tmp += __shfl_down(tmp,delta,64);
            delta += delta;
        }
        if (tid%64==0) atomicAdd(&tmp_sh,tmp);
        __syncthreads();
        if (tid == 0 ) atomicAdd(Integral,tmp_sh);
    }
}

在这段代码中，我们使用了__shfl_down这个函数，它的作用是使一个线程束内，线程ID为id的线程，获取到线程ID为（id+delta）的线程的局部变量tmp的值，64是一个线程束中线程的个数，对于线程id+delta超过64的线程，还是得到自己的tmp值。举例来说，当delta=1时，对于线程0，这个函数返回线程1的tmp；对于线程2，这个函数返回线程3的tmp；对于线程63，这个函数则返回自己的tmp。当delta=2时，对于线程0，此时它的tmp已是进行洗牌前的线程0和线程1的tmp之和，再用__shfl_down把线程2的tmp（洗牌前的线程2和线程3的tmp之和）传给它，加在一起就是洗牌前线程0-3的tmp之和。当delta=4时，会把线程4的tmp（洗牌前线程4-7的tmp之和）传给线程0，加在一起是洗牌前线程0-7的tmp之和。我们对上述过程示意如下，图中的方块代表进行洗牌前的线程上的tmp值，方块中数字是线程ID：

图3-4 使用__shfl_down函数规约示意图

以此类推，随着j不断循环，delta=8时，线程0的tmp为洗牌前线程0-15的tmp之和，delta=16时，是线程0-31的tmp之和，delta=32时，是线程0-63的tmp之和。

通过调用6次__shfl_down，线程束内的线程0的tmp为所有线程的tmp之和，然后再由每个线程束中的线程0（tid模64为0）将tmp值原子加到共享内存tmp_sh中，这样就把每个线程束中64个线程一起对tmp_sh做原子加改成了6次__shfl_down之后由一个线程对tmp_sh做原子加。我们测试一下效果，编译运行（warpshuffle.cc，注意此时blockSize=1024），结果如下：

hipcc warpshuffle.cc -o warpshuffle
./warpshuffle
Numeric result is: 31415.92577126683
Result form formula is: 31415.92653589793
Elapsed time: 49.224998 ms

现在的效果已经优于使用256的线程块但不使用__shfl_down函数的版本了。

3.3.4 另一种可扩展的核函数写法

现在要问，我们的核函数还有进一步的优化空间吗？答案是肯定的，我们可以进一步的减少线程块的数量，你可能会觉得，每个线程块的线程数已经到1024的上限了，再减少线程块的数量，就没有足够的线程来执行N次核函数了。

解决方法很简单，通过循环让每个线程多执行几次核函数的内容就可以了，将核函数改写如下：

__global__ __launch_bounds__(1024)
void integral_kernel(double* Integral, double lb, int N, real dx)
{
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = gridDim.x * blockDim.x;
    int tid = threadIdx.x;
    double tmp = 0.0;
    __shared__ double tmp_sh;
    if (tid == 0) tmp_sh = 0.0;
    __syncthreads();
    for (int i = index; i < N; i += stride) {
        real x = lb + (i + 0.5) * dx;
        real sin_2x = SIN(2 * x);
        sin_2x *= sin_2x;
        real cos_x = COS(x);
        cos_x *= cos_x;
        tmp += sin_2x * cos_x * dx;
    }
    int delta = 1;
    for (int j = 0; j < 6; j++) {
        tmp += __shfl_down(tmp,delta,64);
        delta += delta;
    }
    if (tid%64==0) atomicAdd(&tmp_sh,tmp);
    __syncthreads();
    if (tid == 0 ) atomicAdd(Integral,tmp_sh);
}

上面的核函数中，CPU版本的for循环又回来了，但这次循环计数i并不是从0开始每次递增1，而是对于每个线程来说从index开始，每次递增stride，index是线程在这个核函数启动的所有线程中的索引，stride是这个核函数启动的总线程数。当核函数启动的总线程数少于循环次数N时，每个线程就从循环计数为index的循环开始，负责循环计数相差stride的所有循环，直到把N次循环全部做完。这是另一种常用的可扩展的核函数写法，区别于我们之前在主机代码中根据工作量动态的确定总线程数，这种写法可以用任意的总线程数启动核函数完成工作。

注意线程束洗牌函数和原子加并不在for循环中，因为我们完全可以让每个线程将它所负责的几次循环结果都累加到tmp之后再进行线程束内的求和，以及向共享内存和全局内存做原子加。

我们不改动主机代码中的启动参数，仅修改核函数并编译执行，效果与之前的完全一样的：

hipcc loop_in_kernel.cc -o loop_in_kernel
./loop_in_kernel
Numeric result is: 31415.92577126695
Result form formula is: 31415.92653589793
Elapsed time: 50.818001 ms

现在可以用更少的线程块来执行核函数了，我们尝试让每个线程做100次循环，为此令：

    int numBlocks = (N / 100 + blockSize - 1) / blockSize;

再次编译运行（loop_in_kernel.cc），得到结果：

Numeric result is: 31415.92577126687
Result form formula is: 31415.92653589793
Elapsed time: 4.646000 ms

对这个核函数的优化暂时告一段落。值得一提的是，这种优化方式实际是一种权衡，即到底是让每个线程多做一些计算，而让原子加的次数减少，还是让每个线程少做一些计算，而让原子加的次数增加。这是优化中常常遇到的一种两者不可兼得的场景，需要结合问题具体分析。

练习：使用hipprof工具查看上面几个版本的核函数运行时间。

3.3.5 用空间换时间

对于我们的程序还有一种可行的思路，既然对全局内存做原子加比较慢，我们可以申请长度为线程总数的数组，这样每个线程可以把计算结果保存在不同的位置，避免多个线程对同一位置的原子加，之后再单独对这个数组中的元素求和。这是一种用空间换时间的方法，同样是一种权衡。

首先编写一个使用CPU串行求和的版本，让blockSize等于线程束的大小，这样numBlocks就是线程束的个数，设置长度为numBlocks的数组res_warp（及其设备端对应d_res_warp）用于存储每个线程束的求和结果，最后在让CPU用一个循环将所有线程束的结果求和。这里我们保留了速度较快的线程束洗牌函数将64个线程的结果做一次求和，这样也将数组长度由总线程数减少到了线程束的个数。核函数改写如下：

__global__
void integral_kernel(double* res_warp, double lb, int N,
                     real dx)
{
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = gridDim.x * blockDim.x;
    int tid = threadIdx.x;
    double tmp = 0.0;
    for (int i = index; i < N; i += stride) {
        real x = lb + (i + 0.5) * dx;
        real sin_2x = SIN(2 * x);
        sin_2x *= sin_2x;
        real cos_x = COS(x);
        cos_x *= cos_x;
        tmp += sin_2x * cos_x * dx;
    }
    int delta = 1;
    for (int j = 0; j < 6; j++) {
        tmp += __shfl_down(tmp,delta,64);
        delta += delta;
    }
    if (tid == 0) res_warp[blockIdx.x] = tmp;
}

将主程序改写如下：

int main()
{
    double lower_bound = 0.0;
    double upper_bound = 40000 * acos(-1); //acos(-1) = pi
    int N = 1e8;
    real dx = (upper_bound - lower_bound) / N;
    double Integral = 0.0;
    // 核函数启动配置
    int blockSize = 64;
    int numBlocks = (N / 100 + blockSize - 1) / blockSize;
    // 用于存储每个warp的结果
    int size = numBlocks * sizeof(double);
    double* res_warp = (double *) malloc(size);
    double* d_res_warp;
    hipMalloc(&d_res_warp,size);
    // 开始计时
    my_timer timer;
    timer.start();
    integral_kernel<<<numBlocks, blockSize>>>
                   (d_res_warp,lower_bound,N,dx);
    hipMemcpy(res_warp,d_res_warp,size,hipMemcpyDeviceToHost);
    for (int i = 0; i < numBlocks; i++) {
        Integral += res_warp[i];
    }
    timer.stop();
    printf("Numeric result is: %.11f\n", Integral);
    printf("Result form formula is: %.11f\n",
           F(upper_bound) - F(lower_bound));
    printf("Elapsed time: %8.6f ms\n",timer.dt);
    hipFree(d_res_warp);
    free(res_warp);
    return 0;
}

编译上面的代码（sum_on_cpu.cc）并运行，结果如下：

hipcc sum_on_cpu.cc -o sum_on_cpu
./sum_on_cpu
Numeric result is: 31415.92577126692
Result form formula is: 31415.92653589793
Elapsed time: 4.539000 ms

可以看到性能表现也非常出色，和前面几小节持续优化单个核函数的效果持平了。

题外话：现在的核函数integral_kernel是以计算为主了，请尝试在编译时通过定义-DDOUBLE_PREC控制计算精度，查看程序使用不同精度时的结果，并通过hipprof工具查看核函数性能。DCU还为单精度计算提供了一组内置快速数学函数，以双下划线__为前缀，如__sinf()、__cos()等（详见开发者社区中的DTK文档），请尝试使用并查看精度与性能。

还可以将CPU求和这步也用DCU来做，为此编写一个核函数sum_kernel，代码如下：

__global__
void sum_kernel(double* res_warp, int numBlocks,
                double* Integral)
{
    int index = blockIdx.x * blockDim.x + threadIdx.x;
    int stride = gridDim.x * blockDim.x;
    int tid = threadIdx.x;
    double tmp = 0.0;
#if 1
    index *= 2;
    stride *= 2;
    // 如果numBlocks是奇数
    // 用0号线程把res_warp[numBlocks-1]单独读取
    if (numBlocks%2==1) {
        if (tid==0) tmp = res_warp[numBlocks-1];
        numBlocks -= 1;
    }
    for (int i = index; i < numBlocks; i += stride) {
        // 从全局内存一条指令读取128位，即两个double型
        tmp += res_warp[i] + res_warp[i+1];
    }
#else
    for (int i = index; i < numBlocks; i += stride) {
        tmp += res_warp[i];
    }
#endif
    int delta = 1;
    for (int j = 0; j < 6; j++) {
        tmp += __shfl_down(tmp,delta,64);
        delta += delta;
    }
    if (tid == 0) *Integral = tmp;
}

在这段核函数的编写中我们额外展示了一个小技巧，即DCU存在可以从全局内存或共享内存一次性读取连续的128位的指令，我们通过代码

        tmp += res_warp[i] + res_warp[i+1];

明确的告诉编译器每个线程会读取两个连续的双精度浮点数，刚好128位，编译器就会采用这条指令，达到更好的访存性能。我们保留了普通的写法，可以通过修改条件编译的#if 0和#if 1在这两种方式间切换，通过hipprof性能分析工具的计时来对比性能的差别。

言归正传，在主机代码中，现在不再需要主机内存中的空间res_warp，仅在设备端为d_res_warp分配内存即可，但这次最终的积分结果d_Integral需要拷贝回主机。至于sum_kernel，我们只使用一个线程束来执行，这样在核函数的最后直接用线程0对d_Integral赋值即可。

篇幅所限仅列出计时起止间的主机代码（其余部分仅删除了对res_warp的malloc()和free()）：

    timer.start();
    integral_kernel<<<numBlocks, blockSize>>>
                   (d_res_warp,lower_bound,N,dx);
    sum_kernel<<<1,blockSize>>>
              (d_res_warp,numBlocks,d_Integral);
    hipMemcpy(&Integral,d_Integral,sizeof(double),
              hipMemcpyDeviceToHost);
    timer.stop();

这次编译运行的结果为（sum_on_dcu.cc）：

Numeric result is: 31415.92577126689
Result form formula is: 31415.92653589793
Elapsed time: 4.341000 ms

3.3.6 小节

本节我们通过一个数值计算定积分的例子介绍了使用DCU做并行计算时的规约问题，即计算虽然可以用很多线程并行，但最终所有线程的计算结果需要归并在一起。为了保证规约的正确，需要使用内置原子函数，即在一个线程做读-改-写这三个为一体的不可分的操作时，其他线程均需等待。大量线程对全局内存中的个别地址同时做浮点数原子操作的性能很慢，但是可以利用DCU的内存层次结构，即利用共享内存、线程束洗牌函数进行优化，效果显著。

之后我们还介绍了另一种编写可扩展的核函数的方法，通过使单个线程执行循环的方式实现核函数对计算量的可扩展：即使核函数的启动参数固定，仍然允许计算量动态扩展。我们利用这种写法通过增加单个线程的计算次数，减少了规约次数，使核函数性能进一步优化。

在本节最后，我们介绍了通过空间换时间优化规约的方法，将一部分规约任务或是直接使用CPU、或是让DCU另外运行一个核函数来执行，对于本节的例子同样都有不错的性能表现。这几种优化规约问题的方法，请读者编程时根据具体问题灵活使用。

3.4 异步并发执行

在本节我们将介绍DCU编程中的流水线技术。仍以上节中的积分计算为例，但在本节中将积分区间扩大8倍，即从0到320000ππ。当然根据sin函数和cos函数的周期性，我们知道结果就是上节中从0到40000ππ的8倍，但在这里我们不利用这个性质，而是让DCU挨个区间的算出来。

我们在3.3.5节中使用CPU做一部分求和的版本（sum_on_cpu.cc）基础上稍加修改，DCU核函数不变，只修改主函数main：

int main()
{
    double lower_bound = 0.0;
    double upper_bound = 40000 * acos(-1); //acos(-1) = pi
    double range = upper_bound - lower_bound;
    int N = 1e8;
    real dx = range / N;
    int numIntervals = 8;
    double Integral = 0.0;
    // 核函数启动配置
    int blockSize = 64;
    int numBlocks = (N  + blockSize - 1) / blockSize;
    // 用于存储每个warp的结果
    int size = numBlocks * sizeof(double);
    double* res_warp = (double *) malloc(size);
    double* d_res_warp;
    hipMalloc(&d_res_warp,size);
    // 开始计时
    my_timer timer;
    timer.start();
    for (int j = 0; j < numIntervals; j++) {
        integral_kernel<<<numBlocks, blockSize>>>
                       (d_res_warp,lower_bound,N,dx);
        hipMemcpy(res_warp,d_res_warp,size,
                  hipMemcpyDeviceToHost);
        for (int i = 0; i < numBlocks; i++) {
            Integral += res_warp[i];
        }
        lower_bound += range;
    }
    timer.stop();
    printf("Numeric result is: %.11f\n", Integral);
    printf("Result form formula is: %.11f\n",
           F(range*numIntervals) - F(0));
    printf("Elapsed time: %8.6f ms\n",timer.dt);
    hipFree(d_res_warp);
    free(res_warp);
    return 0;
}

在这段代码中定义了积分区间数numIntervals，积分区间的长度range，此外将调用核函数integral_kernel、设备到主机的拷贝以及主机的求和部分循环8次，每次循环之后将积分区间的下界递增range。此外，为了方便展示本节中要讨论的主题，将线程块个数numBlocks改回到(N + blockSize - 1) / blockSize，这样数组res_warp会更长，因而增加拷贝时间。

3.4.1 用核函数掩盖主机端的求和

在前面我们使用hipprof工具进行性能分析时，只看了统计数据，这次要深入一些，利用它来可视化地观察程序的时间线。编译（integral_async路径下的sum_on_cpu.cc）并使用性能分析工具运行：

hipcc sum_on_cpu.cc -o sum_on_cpu
hipprof ./sum_on_cpu

可以看到路径下生成了几个文件，可能你之前使用时就注意到，其中hip-prof-进程号.db是hipprof记录的原始数据文件，result进程号.hiptrace.csv和result进程号.kernel.csv分别是CSV格式的API和核函数的统计，result_进程号.json则是JSON格式的程序详细记录（你可能会发现文件名中的进程号与可视化界面中的主机线程号Thread并不相等，因为文件名中的进程号是hipprof进程，而我们程序是作为其子进程启动的）。

通过Chrome浏览器（地址栏输入chrome://tracing）或Edge浏览器（edge://tracing）打开这个JSON文件（界面左上角的Load按钮）：

图3-5 使用浏览器打开hipprof生成的JSON文件的界面

可以看到界面分上下两部分，上半部分包括时间标尺和三行时间线，每行时间线都是hipprof记录的一系列事件，用一个个色块表示，色块的长度就是其执行的时间。这三行中的第一行是主机API，第二行是拷贝，第三行则是核函数。图中中间偏右的纵向小面板用于控制鼠标的行为，箭头是选择模式，在这种模式下，单击鼠标可以选中一个色块，界面下半部分就会出现对应事件的具体细节；十字箭头为滚动模式，用于上下左右滚动时间线；前后箭头为缩放模式，用于缩放时间尺度；最下方的是计时模式，用于在时间线上测量时间。

从时间线可以看出，主机在调用核函数之后就开始执行设备到主机的拷贝，但因为拷贝是同步的，因此实际的拷贝是在核函数执行完之后才开始，并且要等拷贝结束后，hipMemcpy才算执行完，这与3.1.6节中对同步异步行为的叙述是一致的。在hipMemcpy之后会有一段空隙（用红框标出），根据代码可知这段是主机在执行res_warp到Integral的累加。在主机执行这段工作时，DCU上也是空闲的，直到下一个循环调用核函数之后才会再次工作起来。

图3-6 程序sum_on_cpu的时间线

思考一下，DCU等待这段时间是没有必要的，它可以立刻继续执行核函数integral_kernel，这只会改变设备上d_res_warp的内容而不会改变主机上res_warp的内容；在DCU执行核函数的同时主机来做对Integral的累加，完成后再将DCU新计算的d_res_warp的结果拷贝过来，这样才更有效率。

当前的代码无法实现这种想法，因为主机是串行执行，只有在当前循环的累加执行完毕后，才能做下一次循环中的核函数调用。为此，一种思路是使用主机的多线程并行，一个线程负责累加，另一个线程调用核函数。但其实使用一个主机也可以完成，修改主机代码如下（只需修改timer.start()和timer.stop()之间的代码）：

    timer.start();
    // 对第一个区间计算并拷贝后，不立刻在主机端求和
    // 而是先启动下一个区间计算的核函数
    integral_kernel<<<numBlocks, blockSize>>>
                   (d_res_warp,lower_bound,N,dx);
    hipMemcpy(res_warp,d_res_warp,size,
              hipMemcpyDeviceToHost);
    lower_bound += range;
    for (int j = 1; j < numIntervals; j++) {
        // 启动核函数计算当前区间
        integral_kernel<<<numBlocks, blockSize>>>
                       (d_res_warp,lower_bound,N,dx);
        // 主机端同时对上一个区间求和
        for (int i = 0; i < numBlocks; i++) {
            Integral += res_warp[i];
        }
        // 求和完毕后，拷贝当前区间的结果
        hipMemcpy(res_warp,d_res_warp,size,
                  hipMemcpyDeviceToHost);
        lower_bound += range;
    }
    // 最后一个区间的求和
    for (int i = 0; i < numBlocks; i++) {
        Integral += res_warp[i];
    }
    timer.stop();

改动很简单，在对j的循环中将拷贝和求和调换一下顺序即可。这样每个循环开始时直接调用核函数在DCU上对当前区间做计算，调用之后主机开始对上一个区间求和（此时res_warp还是上一个区间的结果），求和完毕后再调用hipMemcpy将当前区间的结果拷贝到主机端，这样就实现了DCU上核函数与主机端求和在时间上的重叠。第一个区间和最后一个区间比较特殊，因为对于第一个区间来说没有上一个区间，因此在调用核函数之后就进行拷贝；而最后一个区间的求和没有包含在j的循环中，因此在对j的循环结束后需要补上。

编译运行（hide_cpu.cc）并使用hipprof生成JSON文件：

hipcc hide_cpu.cc -o hide_cpu
hipprof ./hide_copy

打开后时间线显示如下：

图3-7 程序hide_cpu的时间线

仍以红框标出主机上做求和的时间，可以看到在这些时间段内DCU不再空闲，而是在执行integral_kernel核函数了。从整段时间来看，DCU始终在不间断的执行核函数与拷贝，优化效果体现为总耗时从87 ms缩短到了73 ms。

3.4.2 流与异步拷贝

在上节中，我们利用核函数执行与主机执行的异步性来优化程序，使得当主机执行上一个区间求和的同时，DCU也在执行当前区间的计算，使得相邻区间的计算在时间上有重叠，因此缩短了总的计算时间。思考一下，是否能让相邻区间的计算在时间上有更多的重叠呢？比如说让当前区间的核函数与上一个区间的拷贝重叠，甚至与上一个区间的核函数重叠。

我们在3.1.6节的最后提到过，在DCU编程中主机与设备间的拷贝可以是异步的，多个核函数也可以异步执行，在本节我们介绍如何实现这些异步行为，并利用它们进一步优化程序。

在DCU编程中，可以创建被称为“流（stream）”的设备端队列，并将拷贝以及核函数发布在流上。同一个流上的任务是按入列的顺序串行执行的，不同流上的任务则是异步并发的，当然如果一个流的任务已占据了所需的软硬件资源，另一个流的任务则需要等待。

现在看如何将拷贝和核函数发布在指定流上。此前我们使用的拷贝API是hipMemcpy，它有一个异步版本hipMemcpyAsync，相比同步版本多一个流参数。对于核函数，在前面的例子中三个尖括号“<<<>>>”中的核函数启动配置仅包含线程网格和线程块的配置，但启动配置还包括动态共享内存的大小以及流，分别是三个尖括号中的第三和第四个参数。这两个参数的默认值都是0，所以之前没有配置也不影响核函数调用；动态共享内存将在后文介绍，而流参数的默认值代表所使用的是默认流。此前我们的核函数总是发布在默认流上，如上面时间线第三行所标的“Stream0”。

回到我们的程序，简单思考一下你会发现问题，如果将当前区间的核函数与上一个区间的拷贝重叠，会导致上一个区间从d_res_warp向res_warp拷贝的过程中当前区间的核函数integral_kernel同时在修改d_res_warp的内容，这样拷贝的结果很可能出问题；如果将当前区间的核函数与上一个区间的核函数也重叠，则会导致两个核函数同时修改d_res_warp，这同样不行。

有了上面的分析，解决方法自然是为每个区间分配不同d_res_warp，这里我们准备使用双缓冲（double buffering）的方法。定义缓冲结构体Buffer如下：

struct Buffer
{
    double* res_warp;
    double* d_res_warp;
    hipStream_t stream;
    void alloc(int size) {
        hipHostMalloc(&res_warp,size);
        hipMalloc(&d_res_warp,size);
        hipStreamCreate(&stream);
    }
    void dealloc() {
        hipHostFree(res_warp);
        hipFree(d_res_warp);
        hipStreamDestroy(stream);
    }
};

我们将res_warp和d_res_warp都包含在Buffer中，此外还有流变量stream。在Buffer中还有两个成员函数alloc和dealloc，用于申请/释放空间以及创建/销毁流。对于使用hipMalloc和hipFree申请和释放设备内存无需多言；对于主机内存，这里使用了hipHostMalloc和hipHostFree，而非之前的malloc和free，原因是对于涉及主机端内存的拷贝，如果想要异步，则主机端内存必须是页锁定的（page-locked），否则即使使用异步版本的拷贝hipMemcpyAsync，拷贝仍然是同步的。目前用到的hipHostMalloc和hipHostFree是专门用来申请与释放页锁定主机内存的API，关于页锁定主机内存会在下节专门介绍。剩下的hipStreamCreate和hipStreamDestory两个API则是关于流的，分别用于创建与销毁。

定义结构体Buffer后，我们将代码中的res_warp和d_res_warp的创建与内存分配替换为：

    Buffer buffer[2];
    buffer[0].alloc(size);
    buffer[1].alloc(size);

将释放res_warp和d_res_warp替换为：

    buffer[0].dealloc();
    buffer[1].dealloc();

并且将计时器起止间的代码如下修改：

timer.start();
integral_kernel<<<numBlocks, blockSize,0,buffer[0].stream>>>
               (buffer[0].d_res_warp,lower_bound,N,dx);
hipMemcpyAsync(buffer[0].res_warp,buffer[0].d_res_warp,size,
               hipMemcpyDeviceToHost,buffer[0].stream);
lower_bound += range;
int b1, b2; // buffer索引
for (int j = 1; j < numIntervals; j++) {
    b1 = j%2; // 当前区间使用的buffer
    b2 = (j-1)%2; // 上一个区间使用的buffer
    integral_kernel<<<numBlocks, blockSize,0,buffer[b1].stream>>>
                       (buffer[b1].d_res_warp,lower_bound,N,dx);
    hipMemcpyAsync(buffer[b1].res_warp,buffer[b1].d_res_warp,size,
                   hipMemcpyDeviceToHost,buffer[b1].stream);
    // 等待上个区间的结果拷贝完毕
    hipStreamSynchronize(buffer[b2].stream);
    for (int i = 0; i < numBlocks; i++) {
        Integral += buffer[b2].res_warp[i];
    }
    lower_bound += range;
}
// 最后一个区间的求和
hipStreamSynchronize(buffer[b1].stream);
for (int i = 0; i < numBlocks; i++) {
    Integral += buffer[b1].res_warp[i];
}
timer.stop();

在这段代码中，我们根据循环计数j确定使用的Buffer，在每次循环中主机先将当前积分区间的核函数和拷贝发布在该Buffer中的流上，之后调用流同步接口hipStreamSynchronize让主机等待上一个区间的核函数和拷贝结束（使用另一个Buffer中的流和内存空间）。这里的流同步是必要的，因为此时的拷贝是异步的；如果不加同步，主机会在d_res_warp中的结果还未拷贝到res_warp前就执行res_warp到Integral的累加。hipStreamSynchronize只会让主机等待指定的流，只要指定流上的任务完成后主机就可以继续往下执行，而与其他流的执行状态无关。

对于第一个区间和最后一个区间同样需要单独处理，不再赘述。

尝试：将对上一个区间的流同步与当前区间的核函数和拷贝调换顺序，并使用hipprof观察时间线。

编译（double_buffer.cc），并使用hipprof运行，因为我们使用了两个流，我们通过选项—group-stream将时间线按流来分组，即将拷贝和核函数按流显示在同一行：

hipcc double_buffer.cc -o double_buffer
hipprof --group-stream ./double_buffer

这样生成的JSON文件命名为result_进程号_stream.json，使用浏览器打开（在图中，为了清晰区分，我们修改了拷贝的颜色，原本与核函数是相同的配色）：

图3-8 程序double_buffer的时间线

这次的耗时进一步缩短了，而且的确实现了当前区间的核函数在时间上与上一个区间的核函数和拷贝重叠。附带提醒一下主机端的空当是每次循环中使用主机求和的部分。

你可能觉得耗时从73 ms缩短到69 ms这样的优化效果并不明显，其中有一部分因素在于当使用多流执行时，会有额外的软硬件调度开销，需要一定的“预热”时间；表现为图中Stream1上第一个核函数从主机端调用到设备端开始执行之间大约8 ms的空当，以及Stream1上第一个核函数与第一次拷贝之间大约3.6 ms的空当。在预热之后，两个流就维持高效执行了。由于这里的例子总的执行时间较短，使得预热占比较大，确实对优化效果产生了一定影响，对于执行时间较长的场景（比如增加积分区间数），优化效果会更加显著。

3.4.3 主机-设备间拷贝与页锁定内存

现在我们来仔细分析上节使用两个流异步并发的具体优化效果，本节先从设备到主机的拷贝开始。

建议：将hide_cpu.cc和double_buffer.cc编译并使用hipprof工具运行，在浏览器中打开生成的JSON文件，按照本节实际查看一下。

使用可视化界面中的选择模式，选择DeviceToHost的色块查看具体耗时，对于我这次的执行，8次拷贝的耗时基本上都是0.885 ms；打开hide_cpu的JSON文件，同样查看拷贝的时间，得到第一次拷贝大约4.5 ms，后续7次拷贝大约是2.5 ms，对比下来我们使用异步拷贝要快不少。但这里需要说明，导致拷贝变快的原因并不在于同步或异步，而在于我们使用了页锁定主机内存。

要说明这个问题，我们需要详细介绍主机与设备间的内存拷贝。在3.1.7节我们提到过主机与内存间的拷贝会被CPU卸载给DMA（直接内存访问）硬件，DMA在获取到hipMemcpy指定的源地址和目的地址后，不需要CPU参与就可以通过PCIe传递数据。这与将核函数卸载给DCU类似，都是相对CPU异步执行的。但这里的问题是现代操作系统大都使用虚拟内存技术进行内存管理，虚拟内存被分成称为“页”的固定大小的块，操作系统会根据需要（比如缺页时）将页面从内存换出到磁盘，或将磁盘中的页面换入到内存。DMA是无法禁止操作系统进行页面调度的，因此在DCU执行主机与设备之间的拷贝时，都会将拷贝的内容通过主机上一块页面锁定的缓存区进行中转，以保证传输过程中操作系统不会移动内容所在的页面。换句话说，在DCU内存与用malloc申请的、可以被页面调度的主机内存之间进行内存拷贝，存在主机端将数据在可以被页面调度的主机内存与页面锁定的主机内存间搬运的开销；因为主机上的搬运工作需要CPU来做，导致拷贝是同步的，也使得整个拷贝过程耗时更长。具体到我们的代码，如果没有给res_warp分配页锁定内存，当发布从d_res_warp到res_warp的拷贝后，主机肯定得等待这段拷贝执行完毕，因为DMA只是完成了d_res_warp到页锁定缓存区的拷贝，主机还需要将数据从缓冲区拷贝到res_warp中。

清楚了主机与设备间的拷贝机制后，使用页锁定内存拷贝快的原因就很清楚了。当我们使用hipHostMalloc给res_warp分配内存后，这些内存本来就是页锁定的，因此DMA可以直接将数据从d_res_warp拷贝到res_warp中，省去了主机从缓存区的搬运；此外，由于不需要主机的参与，拷贝也可以是异步的了，这正是上节中强调涉及主机端的异步拷贝，主机端内存必须是页锁定的。顺带提一下，使用页锁定内存还有一种方式，可以通过hipHostRegister接口将使用malloc分配的内存页锁定，使用hipHostUnregister解锁。

现在你可能觉得，既然页锁定内存有这些优势，为什么不把主机端的内存都申请为页锁定内存呢？我们不建议这么做，原因有二，首先是申请页锁定内存本身有一定的时间开销，查看主机端API的耗时，对于我这次运行中两次hipHostMalloc的时间分别是2.198 ms和2.501 ms；另一个原因是如果在主机端分配大量页锁定内存会导致主机端没有足够的可调度的页面，从而影响主机的性能。我们的建议是，仅为需要经常在主机与设备端拷贝的变量分配页锁定内存，分配也尽量在程序的开头与结尾，避免在程序运行中反复分配与释放。

此外，还有两点需要说明，第一是DCU有两个DMA引擎，可以支持同时从DCU传入传出数据，这意味着不同流上的HostToDevice和DeviceToHost拷贝是可以并发的，但如果是同向的拷贝则不能并发；第二，因为拷贝是通过DMA，不需要计算单元参与，因此异步并发拷贝和核函数不会对彼此的性能造成影响，这也正是我们让核函数与拷贝在时间上重叠的意义所在。为了确认这一点，可以观察我们程序的时间线，前7次拷贝的同时都有核函数执行，最后一次拷贝时则没有核函数执行了，但这8次拷贝的时间几乎相同；至于核函数耗时增加（仅使用默认流版本中大约5.985 ms，而使用两个流的版本中核函数最长耗时大约为9 ms）则是由于多个核函数并发所导致的。

3.4.4 事件

在上节的介绍之后，你可能想验证一下，如果仅是核函数与拷贝并发时，核函数是否确实与仅使用默认流时的耗时一致。但目前看想控制核函数仅与另一个流上的拷贝并发也并不容易，因为不同流上的任务在发布之后，执行顺序是不一定的。比如说从程序逻辑来看，主机端肯定先发布当前区间的核函数与拷贝，再发布下个区间的核函数与拷贝；但从时间线来看，下一个区间的核函数却比上一个区间的拷贝先执行，导致与上一个区间的核函数有部分重叠。

对于这个问题，我们可以通过在流中插入事件（event）的方式来解决，它相当于流执行的一个节点。我们不做过多抽象的解释，还是结合代码来说明。首先，我们在Buffer结构体中添加事件成员event，并在alloc和dealloc中加入对它的创建与销毁：

struct Buffer
{
    double* res_warp;
    double* d_res_warp;
    hipStream_t stream;
    hipEvent_t event;
    void alloc(int size) {
        hipHostMalloc(&res_warp,size);
        hipMalloc(&d_res_warp,size);
        hipStreamCreate(&stream);
        hipEventCreate(&event);
    }
    void dealloc() {
        hipHostFree(res_warp);
        hipFree(d_res_warp);
        hipStreamDestroy(stream);
        hipEventDestroy(event);
    }
};

之后，我们对在第一个区间的核函数调用之后插入事件：

    integral_kernel<<<numBlocks, blockSize,0,buffer[0].stream>>>
                   (buffer[0].d_res_warp,lower_bound,N,dx);
    hipEventRecord(buffer[0].event,buffer[0].stream);
    hipMemcpyAsync(buffer[0].res_warp,buffer[0].d_res_warp,size,
                   hipMemcpyDeviceToHost,buffer[0].stream);

在j循环中调用核函数之前，我们插入对事件的等待，并在核函数调用之后也插入事件：

for (int j = 1; j < numIntervals; j++) {
    b1 = j%2; // 当前区间使用的buffer
    b2 = (j-1)%2; // 上一个区间使用的buffer
    //等待上一个区间的核函数执行完毕
    hipStreamWaitEvent(buffer[b1].stream,buffer[b2].event);
    // 或使用
    // hipEventSynchronize(buffer[b2].event);
    integral_kernel<<<numBlocks, blockSize,0,buffer[b1].stream>>>
                   (buffer[b1].d_res_warp,lower_bound,N,dx);
    hipEventRecord(buffer[b1].event,buffer[b1].stream);
    hipMemcpyAsync(buffer[b1].res_warp,buffer[b1].d_res_warp,size,
                   hipMemcpyDeviceToHost,buffer[b1].stream);
    // 等待上个区间的结果拷贝完毕
    hipStreamSynchronize(buffer[b2].stream);
    for (int i = 0; i < numBlocks; i++) {
        Integral += buffer[b2].res_warp[i];
    }
    lower_bound += range;
}

先看一下效果，编译（hide_copy.cc）并使用hipprof运行：

hipcc hide_copy.cc -o hide_copy
hipprof --group-stream ./hide_copy

打开时间线（清晰起见，仍然修改了拷贝的颜色）：

图3-9 程序hide_copy的时间线

可以看到时间线符合我们的初衷，抛开前面的预热，后面的核函数只与前一个区间的拷贝在时间上重叠，而与前一个区间的核函数错开了。选中核函数可以看到执行时间基本都是5.985 ms，和仅使用默认流版本（hide_cpu）的执行时间一致，这也就验证了核函数与拷贝重叠确实不会影响彼此的性能。

现在来介绍上面用到的事件相关API就比较容易了，hipEventRecord的作用是在指定流上记录一个事件，用于监视主机在调用hipEventRecord之前发布到该流上任务的状态，在上面的代码中，我们调用hipEventRecord的作用是监视它前面发布到同一流上的integral_kernel是否完成。hipStreamWaitEvent则是保证主机在调用hipStreamWaitEvent之后发布到该流上的任务必须等待事件前的任务完成后再执行，在上面的代码中，当前区间（使用流buffer[b1].stream）的核函数integral_kernel调用在hipStreamWaitEvent之后，因此它必须等待上一个区间buffer[b2].event记录前发布在buffer[b2].stream上的任务，也就是上一个区间的integral_kernel执行完毕，因此相邻区间的integral_kernel在执行时间上就总是错开的；而主机发布上一个区间的拷贝是在调用hipEventRecord记录事件之后，因此当前区间的核函数不必等待上一个区间的拷贝结束就可以执行，两者在时间上是重叠的。

当然还可以使用时间同步接口hipEventSynchronize，它的作用是让主机等待事件记录前相应流上的任务完成，之后再去执行后续的发布，对于我们的代码，这相当于让主机等上一个区间的核函数执行完再发布当前区间的核函数，同样可以避免两者的重叠。

归结起来，在我们的代码中实际上是利用事件来处理了两个流之间的依赖问题，即我们要求一个流上的核函数要等另一个流上的核函数执行完才能开始。除了这方面的功能外，事件的另一个功能是计时，有一个专门用来计时的API：

hipError_t hipEventElapsedTime(float* ms, hipEvent_t start,
                               hipEvent_t stop)

它将事件start和stop记录的时间间隔以毫秒为单位写入ms所指的地址中，注意在调用这个API之前需要先使用hipEventSynchronize等待事件stop完成。

练习：使用hipEventElapsedTime计时，并与hipprof的计时做比较。

3.4.5 核函数并发执行

现在我们将核函数与拷贝并发搞清楚了，下面将专注于核函数与核函数并发。在3.2.2节中我们提到计算单元中的每个SIMD单元支持并发10个线程束，这些线程束可以来自不同的核函数，这正是多个核函数并发的硬件基础。前面已经对比了核函数不重叠（hide_cpu和hide_copy）中的核函数性能（5.985 ms）与核函数部分重叠时（double_buffer）中的核函数性能（6.590 ~ 9.035 ms）。我们暂且不考虑预热期间的第一个核函数，用浏览器可视化界面中的计时模式（小面板中4个按钮最下面的一个）来测量从第2个核函数开始到第8个核函数结束的耗时，对于hide_copy这段耗时是42.119 ms，而double_buffer中则是41.96 ms。这样看来，我们让核函数与核函数部分重叠获得的收益并不高，出现这种情况的原因是单个核函数已经将计算单元占的比较满了（单看线程数需要1e8个，这就已经大于支持的最大并发线程数163840），这时再并发更多的核函数也很难得到更高的性能，因为硬件资源需要分配给不同核函数。

我们再尝试一个例子，记得再3.3.6节中，我们还编写了使用核函数sum_kernel求和每个线程束的计算结果的版本（sum_on_dcu.cc），我们也将它的积分区间扩大8倍。

练习：先不看下面的代码，仿照对sum_on_cpu.cc的修改，完成对sum_on_dcu.cc的修改，再与下面提供的代码对照。

核函数同样不需要改动，主机代码修改如下：

int main()
{
    double lower_bound = 0.0;
    double upper_bound = 40000 * acos(-1); //acos(-1) = pi
    double range = upper_bound - lower_bound;
    int N = 1e8;
    real dx = range / N;
    int numIntervals = 8;
    double Integral[numIntervals];
    double* d_Integral;
    hipMalloc(&d_Integral,numIntervals*sizeof(double));
    // 核函数启动配置
    int blockSize = 64;
    int numBlocks = (N + blockSize - 1) / blockSize;
    // 用于存储每个warp的结果
    int size = numBlocks * sizeof(double);
    double* d_res_warp;
    hipMalloc(&d_res_warp,size);
    // 开始计时
    my_timer timer;
    timer.start();
    for (int j = 0; j < numIntervals; j++) {
        integral_kernel<<<numBlocks, blockSize>>>
                       (d_res_warp,lower_bound,N,dx);
        sum_kernel<<<1,blockSize>>>
                  (d_res_warp,numBlocks,&d_Integral[j]);
        lower_bound += range;
    }
    hipMemcpy(&Integral,d_Integral,numIntervals*sizeof(double),
              hipMemcpyDeviceToHost);
    // 对所有区间的积分求和
    for (int i = 1; i < numIntervals; i++) {
        Integral[0] += Integral[i];
    }
    timer.stop();
    printf("Numeric result is: %.11f\n", Integral[0]);
    printf("Result form formula is: %.11f\n",
           F(range*numIntervals) - F(0));
    printf("Elapsed time: %8.6f ms\n",timer.dt);
    hipFree(d_Integral);
    hipFree(d_res_warp);
    return 0;
}

在上面的代码中，我们为d_Integral分配了8个双精度浮点数的空间，原因在于核函数sum_kernel中最后是对d_Integral直接赋值而非累加，因此我们对每个积分区间的结果分配一个双精度浮点数。在所有区间计算结束后，使用hipMemcpy一次性将这些结果都拷贝到主机端，最后由主机将结果累加在一起。你完全可以修改一下sum_kernel使其最后对d_Integral累加，这样为d_Integral分配一个双精度浮点数的空间就够了（但别忘了通过hipMemcpy或hipMemset初始化d_Integral为0）。我们采用上面的写法，是想展示在申请设备内存与拷贝时应尽量合并。对于上面的d_Integral，高效的做法是一次性申请64字节的空间，并一次性拷贝这64字节；而非分8次，每次申请/拷贝8个字节，因为在3.4.3节已经介绍过，每次拷贝都会有一系列开销。

言归正传，这次仍然让integral_kernel的总线程数是1e8，编译运行（integral_async路径下的sum_on_dcu.cc）：

hipcc sum_on_dcu.cc -o sum_on_dcu
hipprof ./sum_on_dcu

打开时间线如下（我们修改了sum_kernel色块的颜色）：

图3-10 程序sum_on_dcu的时间线

统计数据中sum_kernel平均耗时是6.104 ms，integral_kernel是5.864 ms。

下面我们同样使用双缓冲的方法将程序改为使用两个流。这次的结构体Buffer只需要d_res_warp和stream两个成员即可：

struct Buffer
{
    double* d_res_warp;
    hipStream_t stream;
    void alloc(int size) {
        hipMalloc(&d_res_warp,size);
        hipStreamCreate(&stream);
    }
    void dealloc() {
        hipFree(d_res_warp);
        hipStreamDestroy(stream);
    }
};

对于主机端代码的修改，相比于double_buffer.cc要容易许多：

struct Buffer
{
    double* d_res_warp;
    hipStream_t stream;
    void alloc(int size) {
        hipMalloc(&d_res_warp,size);
        hipStreamCreate(&stream);
    }
    void dealloc() {
        hipFree(d_res_warp);
        hipStreamDestroy(stream);
    }
};
对于主机端代码的修改，相比于double_buffer.cc要容易许多：
// 用于存储每个warp的结果
int size = numBlocks * sizeof(double);
Buffer buffer[2];
buffer[0].alloc(size);
buffer[1].alloc(size);
// 开始计时
my_timer timer;
timer.start();
int b; // buffer索引
for (int j = 0; j < numIntervals; j++) {
    b = j%2; // 当前区间使用的buffer
    integral_kernel<<<numBlocks, blockSize,0,buffer[b].stream>>>
                   (buffer[b].d_res_warp,lower_bound,N,dx);
    sum_kernel<<<1,blockSize,0,buffer[b].stream>>>
              (buffer[b].d_res_warp,numBlocks,&d_Integral[j]);
    lower_bound += range;
}
hipMemcpy(&Integral,d_Integral,numIntervals*sizeof(double),
          hipMemcpyDeviceToHost);
// 对所有区间的积分求和
for (int i = 1; i < numIntervals; i++) {
    Integral[0] += Integral[i];
}
timer.stop();
printf("Numeric result is: %.11f\n", Integral[0]);
printf("Result form formula is: %.11f\n",
       F(range*numIntervals) - F(0));
printf("Elapsed time: %8.6f ms\n",timer.dt);
buffer[0].dealloc();
buffer[1].dealloc();
return 0;

上面的代码不需要像double_buffer.cc中那样包含复杂的逻辑，在j循环中既发布当前区间的任务又等待上一个区间的任务，这是因为这里的求和是在设备端而非主机端，不需要在循环中添加同步，因而不必担心主机因等待设备端的执行导致无法继续发布任务。上面代码的j循环中，只需要让主机一前一后的调用计算当前区间的integral_kernel和sum_kernel，并通过索引b让相邻区间的计算发布在不同流上即可。

编译并运行（kernel_overlap.cc）：

hipcc kernel_overlap.cc -o kernel_overlap
hipprof ./kernel_overlap

打开时间线：

图3-11 程序kernel_overlap的时间线

这次使用两个流的收益就大多了，如果抛开预热的因素，只算第一次sum_kernel开始到最后一次sum_kernel的耗时，这里是48.636 ms，而sum_on_dcu中则是90.685 ms，几乎快了一倍。从统计数据来看，sum_kernel平均耗时是6.035 ms，integral_kernel是6.047 ms；integral_kernel的平均耗时相比于sum_on_dcu中的只增加了约0.2 ms，却换来与sum_kernel在执行时间上完美的重叠。这里的原因在于sum_kernel使用硬件资源很少（只用了一个线程束），因此与integral_kernel并发时对integral_kernel的性能影响也很小（相当于DCU上少了一个线程束）。

尝试：通过修改numBlocks使sum_on_dcu和kernel_overlap中integral_kernel的总线程数是1e5（小于最大并发线程数163840），编译运行并查看并发效果并分析原因。

最后，对于多个核函数并发，我们做一些总结说明。当核函数通过多流在DCU上并发执行时，硬件会根据核函数的资源需求和硬件当前使用情况进行调度，当硬件资源已占满时，硬件会以在并发的核函数间来回切换的方式执行，这样核函数的执行时间相比单独执行时要长，具体长多少取决于并发的几个核函数各自占用的硬件资源，在这种情况下多个核函数并发并不一定会带来更高的性能。顺便提一句，CPU在有大量进程同时运行时也会使用这样的方式，出现多个进程仍在同时在运行，但进程都变得卡顿的现象，不过相比于CPU，DCU在核函数间切换的开支要小得多。但当一个核函数没有占满硬件资源时，多个核函数并发通常能带来性能的大幅提升。

当我们在编写DCU程序时，很多时候未必清楚一个核函数是否能占满硬件资源，而且这种情况还常常随执行时问题规模的动态变化而改变，这时将多个核函数并发执行，将问题扔给硬件调度通常是不错的选择，如果硬件资源已满，结果通常也就是提速不明显，如果硬件资源没用满，就会有很好的收益。

不过仍需提醒，这里依旧存在着权衡，比如在本节中见到的，多流调度起始阶段有一定的预热时间，此外使用两个流时需要双缓冲区，如果使用三个流则需要三缓冲区，等等，仍需仔细权衡。

3.4.6 默认流

在本节的最后，我们来看将核函数发布到默认流会怎样。我们将kernel_overlap.cc中sum_kernel的核函数执行配置改为：

sum_kernel<<<1,blockSize,0,0>>>

编译运行并打开时间线：

图3-12 将sum_kernel发布到默认流时kernel_overlap的时间线

可以看到，在默认流执行核函数并不会导致主机端同步，即主机线程发布核函数后会继续执行并发布其他任务，而不会等待默认流的核函数执行完毕后再继续；但默认流会让设备端的执行同步，一旦核函数开始在默认流执行，其他流上的核函数必须等待其执行完毕后才能开始。当然你可以利用默认流的性质，使得个别核函数总是单独执行，比如说以此保证它在执行时数据不会被其他核函数搞乱，但是在异步并发的场景下使用默认流会让设备上的执行流程变得更为复杂。

3.4.7 小节

在本节中，我们介绍了如何利用异步并发执行来优化程序。首先，利用了核函数执行与主机执行的异步性，避免了在主机执行计算时DCU处于不必要的空闲状态，以此提升了程序性能。之后，介绍了如何将核函数和拷贝发布在流上以实现异步并发，我们介绍了对于拷贝和核函数由于使用的是不同的硬件单元，因此在并发执行时两者可以在时间上重叠而不会影响彼此的性能，达到很好的优化效果；而对于核函数与核函数的并发，则相当于我们将多个核函数一块交给硬件，由硬件来调度，对于使用硬件资源较多的核函数，效果不一定很好，但如果并发的核函数都只用部分的硬件资源则可以得到很好的收益。

本节中为了能异步执行主机与设备间的拷贝，我们还介绍了页锁定主机内存，并以此为引介绍了主机与设备间拷贝的详细机制；为了处理不同流间的依赖，我们介绍了事件的用法；最后我们还展示了默认流上执行的核函数的行为，这些都是比较重要的。

最后我们建议如果是初学者在程序中使用多流并发时，尽量少的包含同步（除hipStreamSynchronize、hipDeviceSynchronize外，还包括hipMalloc、hipMemcpy、hipMemset、hipFree这些同步的接口），它们会导致主机无法持续向设备发布任务，使并发状态被破坏，还使得代码本身的逻辑以及实际执行的逻辑变得复杂。对于初学者来说，本节中kernel_overlap程序是比较好的范例。

3.5 矩阵乘法

本节我们要介绍如何使用DCU计算矩阵乘，很多领域的计算都会涉及到矩阵乘。对于一个m×n 的矩阵A，其矩阵元记为Aij，一个n×p的矩阵B，其矩阵元为Bij，A与B相乘，得到m×p的矩阵C，其矩阵元为：

即每计算CC的一个矩阵元，就需要做nn次乘法和n−1次加法，所以整个矩阵乘需要mnp次乘法和mp(n−1)次加法。对于AA、BB、CC都是n×n方阵的情况，这是一个O(n3)复杂度的问题，计算量会随nn快速增长。观察上面的公式，我们发现每个矩阵元Cij的计算是相互独立的，因此非常适合用并行计算来加速，但想要最大化的利用DCU的性能，在编写DCU程序时又有很多值得关注的点。像前几节一样，我们会在本节中编写DCU程序计算矩阵乘并不断改进，以此进一步介绍DCU编程。

们先简单编写一版CPU代码（matrixMul目录下的cpu.cc），使用my_timer计时看看效果，简单起见，让矩阵都是N维方阵：

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "timer.h"

void matrix_multiply(int N, double *A, double *B, double *C)
{
    int i, j, k;
    for (j = 0; j < N; j++) {
        for (i = 0; i < N; i++) {
            double Cij = 0.0;
            for (k = 0; k < N; k++) {
                // C_{ij} = \sum_k A_{ik} * B_{kj}
                Cij += A[i + k*N] * B[k + j*N];
            }
            C[i + j*N] = Cij;
        }
    }
}

void print_matrix(int N, double* M)
{
    int i, j;
    int max_size = 5; // 打印的最大维度
    for (i = 0; i < N && i < max_size; i++) {
        for (j = 0; j < N && j < max_size; j++) {
            printf("%13.6f ", M[i + j*N]);
        }
        if (N > max_size) { // 超过了最大维度，打印省略号
            printf("... ");
        }
        printf("\n");
    }
    if (N > max_size) { // 超过了最大维度，打印省略号
        printf(".\n.\n.\n");
    }
}

void generate_matrix(int N, double* M)
{
    int i, j;
    for (j = 0; j < N; j++) {
        for (i = 0; i < N; i++) {
            M[i + j*N] = (double) rand() / (double) RAND_MAX;
        }
    }
}

int main(int argc, char **argv)
{
    if (argc != 2) {
        fprintf(stderr, "Usage: %s N\n", argv[0]);
        return 1;
    }
    int N = atoi(argv[1]);
    if (N < 2) {
        fprintf(stderr, "N must larger than 2.\n");
        return 1;
    }

    // 初始化随机数发生器
    srand(time(NULL));

    // 计时器
    my_timer timer;

    // 分配矩阵内存空间
    double *A = (double *) malloc(N * N * sizeof(double));
    double *B = (double *) malloc(N * N * sizeof(double));
    double *C = (double *) malloc(N * N * sizeof(double));

    // 随机生成矩阵A和B
    timer.start();
    generate_matrix(N, A);
    generate_matrix(N, B);
    timer.stop();
    printf("Generating Matrices took %8.3f ms.\n", timer.dt);

    // 打印矩阵A,B
    printf("Matrix A:\n");
    print_matrix(N, A);
    printf("\n");
    printf("Matrix B:\n");
    print_matrix(N, B);
    printf("\n");

    // 矩阵乘法，将结果存储在矩阵C中
    timer.start();
    matrix_multiply(N, A, B, C);
    timer.stop();
    printf("Matrix multiplication took %8.3f ms.\n", timer.dt);

    // 打印矩阵C
    printf("Matrix C:\n");
    print_matrix(N, C);

    // 释放内存
    free(A);
    free(B);
    free(C);

    return 0;
}

在这段代码里，我们没有使用C语言中的二维数组来存储矩阵，而是视为一维数组，并且采用了列主序的存储方式（注意函数print_matrix和matrix_multiply中引用矩阵元的方式），即同一列的元素在内存中是连续存储的。这是有历史原因的，因为诸如BLAS、LAPACK这些线性代数库都是使用Fortran语言编写的，而Fortran语言是列主序的。我们使用列主序可以方便在后续代码中使用按照BLAS接口标准的库。

矩阵A和B是使用函数generate_matrix随机生成的，利用了C语言标准库中的rand()函数（包含在头文件stdlib.h中，作用是产生0到RAND_MAX间的伪随机数），使得所有矩阵元都是0~1之间的双精度浮点数。在生成矩阵后，使用函数print_matrix将矩阵输出到屏幕上（对于大矩阵，只打印前5行5列的子矩阵）。之后使用matrix_multiply函数将A和B相乘并将结果保存在矩阵C中，之后打印结果。在matrix_multiply中没有采用任何优化手段，仅是按照矩阵乘法的公式编写了最简单的实现。

编译并使用我们的程序计算两个2×22×2矩阵相乘：

g++ cpu.cc -o cpu -O3
./cpu 2
Generating Matrices took    0.001 ms.
Matrix A:
     0.041022      0.575630 
     0.975030      0.516855 

Matrix B:
     0.170363      0.635518 
     0.882355      0.999212 

Matrix multiplication took    0.001 ms.
Matrix C:
     0.514898      0.601246 
     0.622158      1.136096

可以使用计算器验证一下，我们的程序是没有问题的。再试试两个1000×10001000×1000的矩阵相乘：

./cpu 1000
Generating Matrices took   18.755 ms.
Matrix A:
     0.404158      0.718073      0.138790      0.999525      0.904998 ... 
     0.200418      0.410106      0.639114      0.602230      0.319843 ... 
     0.208472      0.098193      0.212182      0.299061      0.239892 ... 
     0.331144      0.516972      0.145787      0.030868      0.595335 ... 
     0.812551      0.437160      0.547514      0.278571      0.824977 ... 
.
.
.

Matrix B:
     0.647890      0.089632      0.320305      0.104900      0.003846 ... 
     0.721201      0.027903      0.387046      0.265118      0.312013 ... 
     0.775124      0.951991      0.262372      0.277300      0.874698 ... 
     0.663479      0.147379      0.351739      0.047627      0.450676 ... 
     0.369385      0.029823      0.998284      0.239315      0.682138 ... 
.
.
.

Matrix multiplication took 1257.347 ms.
Matrix C:
   256.739709    251.385397    271.003417    258.734088    251.304090 ... 
   248.442402    237.363746    258.863536    246.663164    244.148036 ... 
   235.748062    238.223707    248.997626    243.717928    242.310350 ... 
   249.627830    246.672617    260.516637    255.078815    246.545102 ... 
   240.161690    238.072973    249.207898    249.527630    243.100153 ... 
.
.
.

尝试：使用我们的程序计算更多方阵的乘积（建议不要过大），并尝试不同的编译器优化级别，观察耗时与矩阵规模N的关系，思考当矩阵规模越来越大时程序有哪些瓶颈。

3.5.1 使用DTK中的库

下面我们开始编写DCU程序，在2.1.2节中，我们提到过DTK软件栈中有很多库，这里我们就先使用DTK软件栈中的hipRAND库来随机生成矩阵，再使用hipBLAS库来计算矩阵乘。这样，一则简单的展示一下怎样使用DTK软件栈中的库，二来也看看我们在DCU上能达到的上限，并作为我们自己编写的DCU程序的优化目标。

编写DCU程序如下（use_libs.cc）：

#include <stdio.h>
#include <hip/hip_runtime.h>
#include <hiprand/hiprand.h>
#include <hipblas.h>
#include <time.h>

void print_matrix(int N, double* M)
{
    int i, j;
    int max_size = 5; // 打印的最大维度
    for (i = 0; i < N && i < max_size; i++) {
        for (j = 0; j < N && j < max_size; j++) {
            printf("%13.6f ", M[i + j*N]);
        }
        if (N > max_size) { // 超过了最大维度，打印省略号
            printf("... ");
        }
        printf("\n");
    }
    if (N > max_size) { // 超过了最大维度，打印省略号
        printf(".\n.\n.\n");
    }
}

int main(int argc, char **argv)
{
    if (argc != 2) {
        fprintf(stderr, "Usage: %s N\n", argv[0]);
        return 1;
    }
    int N = atoi(argv[1]);
    if (N < 2) {
        fprintf(stderr, "N must larger than 2.\n");
        return 1;
    }

    // 初始化hipRAND随机数发生器
    hiprandGenerator_t gen;
    hiprandCreateGenerator(&gen, HIPRAND_RNG_PSEUDO_DEFAULT);
    hiprandSetPseudoRandomGeneratorSeed(gen, time(NULL));

    // 计时器
    hipEvent_t start, stop;
    hipEventCreate(&start);
    hipEventCreate(&stop);
    float dt;

    // 分配矩阵内存空间
    double *A = (double *) malloc(N * N * sizeof(double));
    double *B = (double *) malloc(N * N * sizeof(double));
    double *C = (double *) malloc(N * N * sizeof(double));
    double *dA, *dB, *dC;
    hipMalloc(&dA, N * N * sizeof(double));
    hipMalloc(&dB, N * N * sizeof(double));
    hipMalloc(&dC, N * N * sizeof(double));

    // 随机生成矩阵A和B
    hipEventRecord(start);
    hiprandGenerateUniformDouble(gen, dA, N * N);
    hiprandGenerateUniformDouble(gen, dB, N * N);
    hipEventRecord(stop);
    hipEventSynchronize(stop);
    hipEventElapsedTime(&dt, start, stop);
    printf("Generating Matrices took %8.3f ms.\n", dt);

    // 打印矩阵A,B
    hipMemcpy(A, dA, N * N * sizeof(double),
              hipMemcpyDeviceToHost);
    hipMemcpy(B, dB, N * N * sizeof(double),
              hipMemcpyDeviceToHost);
    printf("Matrix A:\n");
    print_matrix(N, A);
    printf("\n");
    printf("Matrix B:\n");
    print_matrix(N, B);
    printf("\n");

    // 创建hipBLAS句柄
    hipblasHandle_t handle;
    hipblasCreate(&handle);

    // 矩阵乘法，将结果存储在矩阵C中
    const double alpha = 1.0;
    const double beta = 0.0;
    // 预热
    hipblasDgemm(handle, HIPBLAS_OP_N, HIPBLAS_OP_N, 2, 2, 2,
                 &alpha, dA, N, dB, N, &beta, dC, N);
    // 正式计算
    hipEventRecord(start);
    hipblasDgemm(handle, HIPBLAS_OP_N, HIPBLAS_OP_N, N, N, N,
                 &alpha, dA, N, dB, N, &beta, dC, N);
    hipEventRecord(stop);
    hipEventSynchronize(stop);
    hipEventElapsedTime(&dt, start, stop);
    printf("Matrix multiplication took %8.3f ms.\n", dt);

    // 打印矩阵C
    hipMemcpy(C, dC, N * N * sizeof(double),
              hipMemcpyDeviceToHost);
    printf("Matrix C:\n");
    print_matrix(N, C);

    // 清理
    hipFree(dA);
    hipFree(dB);
    hipFree(dC);
    hipEventDestroy(start);
    hipEventDestroy(stop);
    hiprandDestroyGenerator(gen);
    hipblasDestroy(handle);
    free(A);
    free(B);
    free(C);

    return 0;
}

在上面的代码中，我们为矩阵在DCU上申请了内存dA、dB、dC。对于库的使用方面，代码中首先include了用到的hipRAND库和hipBLAS库的头文件，创建了hipRAND的随机数生成器gen，并为其设置了伪随机数种子，之后使用了hiprandGenerateUniformDouble函数，其作用是产生0~1之间均匀分布的双精度伪随机数，我们用它生成了dA和dB的矩阵元；对于矩阵乘，我们先创建了hipBLAS的句柄handle，之后使用hipBLAS的库函数hipblasDgemm计算dA与dB的乘积。hipblasDgemm的功能更具普遍性，可以对m×n的矩阵Op(A)、n×k的矩阵Op(B)，和m×k的矩阵C计算形式为C:=αOp(A)Op(B)+βC的问题，其中Op()可以是转置或转置复共轭（Herimtian共轭），也可以不做操作，这里我们只是用它计算了很简单的情况。

我们没有使用之前的my_timer计时器，而是利用了3.4.4事件一节中提到的hipEventElapsedTime接口来计时。注意到，我们在调用参数m、n、k都为N的hipblasDgemm之前，先令m、n、k都为2调用了一次，这是为了计时准确而做的“预热”。因为矩阵乘法非常重要，而其性能又同矩阵规模密切相关，因此hipBLAS中包含了大量的核函数，用于在程序运行时，根据hipblas?gemm（"?"代表不同精度，目前支持半精度H、单精度S、双精度D、单精度复数C和双精度复数Z）的规模选择性能最优的核函数来执行，而当前版本的hipBLAS库会在第一次调用hipblas?gemm时将这些核函数全部加载到内存里，在我的计算机里这需要大约3秒的时间，我们需要将这段时间排除。

使用hipcc编译，这次需要链接hipBLAS和hipRAND的动态库（此外低于23.04的DTK版本需要额外指定头文件路径-I$ROCM_PATH/rocrand/include）：

hipcc use_libs.cc -o use_libs -I$ROCM_PATH/rocrand/include -L$ROCM_PATH/lib -lhipblas -lhiprand

运行程序use_lib，令N等于4096：

Generating Matrices took    5.450 ms.
Matrix A:
     0.798041      0.043402      0.178880      0.356740      0.575366 ... 
     0.836317      0.664488      0.877761      0.791172      0.515096 ... 
     0.499405      0.304049      0.046138      0.740579      0.088408 ... 
     0.500394      0.360916      0.580063      0.067178      0.363289 ... 
     0.010420      0.664991      0.032631      0.443933      0.858081 ... 
.
.
.

Matrix B:
     0.708751      0.067173      0.965586      0.871611      0.974055 ... 
     0.552622      0.625481      0.081746      0.607098      0.587923 ... 
     0.742625      0.765052      0.133508      0.116854      0.348009 ... 
     0.044137      0.759527      0.606604      0.640436      0.216237 ... 
     0.386603      0.200040      0.233368      0.187415      0.674187 ... 
.
.
.

Matrix multiplication took   28.288 ms.
Matrix C:
  1015.087370   1021.600700   1007.313120   1033.348788   1026.125069 ... 
  1048.421677   1032.788823   1011.338992   1035.188871   1021.101648 ... 
  1010.991994   1013.646259   1018.475126   1015.856618   1020.496487 ... 
  1004.955157   1020.468667   1015.666630   1010.192140   1015.513084 ... 
  1012.398788   1012.408923   1029.628456    997.582603   1015.931811 ... 
.
.
.

我们来计算一下，按DCU在28.288 ms内计算了4096340963次乘加运算（乘加算做了两次浮点运算），DCU达到了4.858 TFLOPS（floating point operations per second）的性能。

令N等于16384，hipblasDgemm耗时：

Matrix multiplication took 1736.475 ms.

这次DCU达到了5.065 TFLOPS的性能。

尝试使用BLAS、OpenBLAS、MKL等使用CPU计算的线性代数库中的Dgemm函数，与hipBLAS做一些性能对比。

3.5.2 初步编写核函数

我们编写一个简单的核函数，看看与hipblasDgemm的性能差距。

__global__
void matrix_multiply(int N, double *A, double *B, double *C)
{
    int tidx = blockIdx.x * blockDim.x + threadIdx.x;
    int stridex = gridDim.x * blockDim.x;
    int tidy = blockIdx.y * blockDim.y + threadIdx.y;
    int stridey = gridDim.y * blockDim.y;
    int i, j, k;
    for (j = tidy; j < N; j += stridey) {
        for (i = tidx; i < N; i += stridex) {
            double Cij = 0.0;
            for (k = 0; k < N; k++) {
                // C_{ij} = \sum_k A_{ik} * B_{kj}
                Cij += A[i + k*N] * B[k + j*N];
            }
            C[i + j*N] = Cij;
        }
    }
}

删去打印矩阵的函数，转而添加一个主机端函数来比较核函数与hipblasDgemm的计算结果：

void compare_matrices(int N, double* C1, double* C2)
{
    int i, j;
    double max_error = 0.0;
    for (j = 0; j < N; j++) {
        for (i = 0; i < N; i++) {
            double error = fabs(C1[i+j*N] - C2[i+j*N]);
            if (error > max_error) {
                max_error = error;
            }
        }
    }
    printf("Max error: %.4e\n", max_error);
}

在main函数中，申请C_kernel和C_hipblas分别存储核函数与hipBLAS的计算结果：

    double *C_kernel = (double *) malloc(N * N * sizeof(double));
    double *C_hipblas= (double *) malloc(N * N * sizeof(double));

计算部分改动如下：

    // 正式计算
    hipEventRecord(start);
    hipblasDgemm(handle, HIPBLAS_OP_N, HIPBLAS_OP_N, N, N, N,
                 &alpha, dA, N, dB, N, &beta, dC, N);
    hipEventRecord(stop);
    hipEventSynchronize(stop);
    hipEventElapsedTime(&dt, start, stop);
    printf("Matrix multiplication (hipblasDgemm)"
           "took %8.3f ms.\n", dt);

    // 将hipBLAS计算的结果拷贝回主机
    hipMemcpy(C_hipblas, dC, N * N * sizeof(double),
              hipMemcpyDeviceToHost);

    // 将dC清零
    hipMemset(dC, 0, N * N * sizeof(double));

    // 使用核函数计算
    dim3 blockSize = dim3(Nsub, Nsub, 1);
    int Nblk = (N + Nsub - 1) / Nsub;
    dim3 numBlocks = dim3(Nblk, Nblk, 1);
    hipEventRecord(start);
    matrix_multiply<<<numBlocks, blockSize>>>(N, dA, dB, dC);
    hipEventRecord(stop);
    hipEventSynchronize(stop);
    hipEventElapsedTime(&dt, start, stop);
    printf("Matrix multiplication (kernel)"
           "took %8.3f ms.\n", dt);

    // 将核函数计算的结果拷贝回主机并与hipBLAS的结果做比较
    hipMemcpy(C_kernel, dC, N * N * sizeof(double),
              hipMemcpyDeviceToHost);
    compare_matrices(N, C_kernel, C_hipblas);

先使用hipblasDgemm计算dA与dB相乘，保存在dC中，计算完成后将dC中的结果拷贝到C_hipblas；接着使用核函数matrix_multiply再算一遍，计算完成后将结果拷贝到C_kernel中，最后用compare_matrices逐个元素的比较两者的计算结果。

图3-13示意了核函数matrix_multiply的线程分工，我们使用2维的线程块和线程网格，线程索引的x分量对应矩阵的行，y分量对应矩阵的列，每个线程负责计算C的一个矩阵元Cij，为此需要读取矩阵A的i行和矩阵B的j列，并做乘累加。目前我们先使用（8，8，1）的线程块（一个线程束）来启动核函数，为此在文件一开始定义Nsub为8：

#define Nsub 8

我们启动足够多的线程块，让每个线程计算一个矩阵元即可。从分块矩阵乘的角度，每个线程块计算了矩阵A的8×N的子矩阵与矩阵B的N×8的子矩阵的乘积，得到矩阵C的8×8的子矩阵（因此将线程块单个维度命名为Nsub）。

图3-13 核函数matrix_multiply

编译运行如下（dcu.cc）：

hipcc dcu.cc -o dcu -I$ROCM_PATH/rocrand/include -L$ROCM_PATH/lib -lhipblas -lhiprand
./dcu 4096
Generating Matrices took    7.533 ms.
Matrix multiplication (hipblasDgemm) took   28.389 ms.
Matrix multiplication (kernel) took  478.200 ms.
Max error: 0.0000e+00

我们编写的核函数的计算结果与hipBLAS库的完全一致，但是性能相去甚远，需要设法改进。观察图3-13，我们很快会发现多个线程会使用矩阵A和B的同一部分，即同一行的N个线程（tidx相同，tidy不同）都需要矩阵A的i行，因此A的每一行都需要从全局内存读取N次，同一列的N个线程（tidy相同，tidx不同）都需要矩阵B的j列，因此B的每一列也都需要从全局内存读取N次。这样算下来，整个核函数执行下来，矩阵A和B的每个矩阵元都需要从全局内存读取到寄存器N次，当然在这些读取中DCU的多级缓存结构会起到作用，这些矩阵元并非每次都得从片外内存获取，但这一过程完全由硬件决定，我们是无法控制的。如果我们想要自己将反复用到的矩阵元缓存下来，需要使用共享内存。

3.5.3 使用共享内存优化

在上节中，我们看到按目前的核函数写法和执行配置（包含一个线程束的8×8的线程块），每个线程块是在做矩阵A的8行构成的子矩阵（8×N）和矩阵B的8列构成的子矩阵（N×8）的矩阵乘，得到矩阵C的一个8×8的子矩阵。在计算过程中A的行和B的列会被反复使用，每次都需要从全局内存中读取。本节我们尝试使用共享内存利用这种数据复用的特点对程序进行优化，先将子矩阵从全局内存读取到共享内存，之后多个线程用到同一矩阵元都从共享内存中读取，以利用共享内存相比全局内存更快的访存速度。

对于我们目前的程序，一个简单的想法是将上述这两个子矩阵全部保存到共享内存中，首先计算一下所需的内存大小，对于矩阵元为双精度浮点数的两个子矩阵需要2×8×N×8B=128NB的空间。在3.2.2节中我们提到每个计算单元有64KB的共享内存，所以只有当N小于等于512时才有可能一次性将这两个子矩阵读入共享内存，对于N更大的情况，需要将子矩阵进一步分块，因此我们不妨直接从分块矩阵乘的角度来考虑。

对于分块矩阵乘法，公式

然成立，但此时Cij、Aik、Bkj不是一个矩阵元，而是子矩阵。这些子矩阵不能随意划分，要求Aik的列数与Bkj的行数相等，这样两者才能相乘，得到的Cij的行数与Aik的行数相等，列数与Bkj的列数相等。我们按这个公式重新编写核函数，将N×N的矩阵A、B、C划分为Nblock×Nblock个子矩阵，每个子矩阵都是Nsub×Nsub的方阵。

我们先编写一个设备端函数用于将矩阵MM的元素从全局内存读到共享内存中的子矩阵Mij中，每个线程负责将矩阵M中Row行Col列的元素读取到Mij中的row行col列，编写这个函数的关键是搞清楚Row与row、Col与col的关系：

__device__
void read_matrix(double* M, int N, int i, int j,
                 double* M_ij, int LD)
{
    // 矩阵元在子矩阵M_ij中是row行col列
    int row = threadIdx.x;
    int col = threadIdx.y;
    // 矩阵元在M_ij中的偏移
    int offset = row + col * LD;
    // 矩阵元在M中是Row行Col列
    int Row = i * Nsub + row;
    int Col = j * Nsub + col;
    // 矩阵元在M中的偏移
    int Offset = Row + Col * N;
    // 没有越界则读取，否则置零
    if ((Row < N) && (Col < N)) {
        M_ij[offset] = M[Offset];
    } else {
        M_ij[offset] = 0.0;
    }
}

在函数read_matrix中，我们使用了LD（Leading Dimension）这个变量，这同样是按照BLAS中命名方式，在列主序中，这个量是同一行相邻两列元素在内存中的间隔。

注意到，在read_matrix中我们让线程索引的x分量threadIdx.x对应矩阵的行row，这是因为在3.2.2中提到的，让线程ID连续的线程访问内存中连续的地址时，DCU可以将这些访存请求合并，以此提高访存效率。如果反过来，让threadIdx.x对应矩阵的列col，则线程ID连续的线程访问的内存地址间隔可能很远（当N增大时），DCU就无法将这些访存请求合并了。此外，我们考虑了NsubNsub不能整除NN的情况，将MijMij中超过MM边界的元素都置为0。现在可以将核函数改写如下：

__global__
void matrix_multiply(int N, int Nblk, double *A, double *B,
                     double *C)
{
    // 声明A_sub,B_sub用于保存A,B的子矩阵
    __shared__ double A_sub[Nsub * Nsub];
    __shared__ double B_sub[Nsub * Nsub];
    // 每个线程块负责计算子矩阵C_ij
    for (int i = blockIdx.x; i < Nblk; i += gridDim.x) {
    for (int j = blockIdx.y; j < Nblk; j += gridDim.y) {
        // 每个线程负责计算子矩阵C_ij中row行col列的矩阵元
        int row = threadIdx.x;
        int col = threadIdx.y;
        double C_sub = 0.0;
        for (int k = 0; k < Nblk; k++) {
            // 将子矩阵A_ik从全局内存读到共享内存
            read_matrix(A, N, i, k, A_sub, Nsub);
            // 将子矩阵B_kj从全局内存读到共享内存
            read_matrix(B, N, k, j, B_sub, Nsub);
            // 计算矩阵A_ik和B_kj的乘
            for (int l = 0; l < Nsub; l++) {
                C_sub += A_sub[row + l * Nsub]
                         * B_sub[l + col * Nsub];
            }
        } // k
        // 将矩阵元写到C中
        int Row = i * Nsub + row;
        int Col = j * Nsub + col;
        if ((Row < N) && (Col < N)) {
            C[Row + Col * N] = C_sub;
        }
    } // j
    } // i
}

在这个版本的核函数中，每个线程块负责计算子矩阵Cij，每个线程负责计算子矩阵中的一个矩阵元（变量C_sub,位于子矩阵Cij的row行col列）。计算流程是先将子矩阵Aik和子矩阵Bkj从全局内存读取到共享内存，将这两个子矩阵相乘后累加到C_sub上，将k从0至Nblk循环最终完成计算，将结果写回矩阵CC中。

这样，对全局内存的读取变成了线程网格中同一行的Nblk个线程块（blockIdx.x相同，blockIdx.y不同）都需要矩阵A的i行子矩阵，因此A的每一行子矩阵都需要从全局内存读取Nblk次，同一列的Nblk个线程块（blockIdx.y相同，blockIdx.x不同）都需要B的j列子矩阵，因此B的每一列子矩阵也都需要从全局内存读取Nblk次。整个核函数执行下来，矩阵A和矩阵B的每个子矩阵都需要从全局内存读取到共享内存Nblk次，也就是每个矩阵元都需要从全局内存读取到共享内存Nblk次，对比上节最后的分析，可知对全局内存的访问由N次减少为Nblk次，变为之前的1/Nsub1/Nsub，而线程块内多个线程对同一矩阵元的重复读取都是从共享内存中读的。

注意核函数的参数中增加了矩阵在单个维度的分块数Nblk，因此修改main函数中对核函数的调用为：

    matrix_multiply<<<numBlocks, blockSize>>>
                    (N, Nblk, dA, dB, dC);

编译运行（shared.cc）如下：

hipcc shared.cc -o shared -I$ROCM_PATH/rocrand/include -L$ROCM_PATH/lib -lhipblas -lhiprand
./shard 4096
Generating Matrices took    5.643 ms.
Matrix multiplication (hipblasDgemm) took   28.706 ms.
Matrix multiplication (kernel) took  465.965 ms.
Max error: 0.0000e+00

效果不明显，我们试试增大线程块的大小Nsub，进一步减少对全局内存的读取，定义Nsub为16。但注意此时我们需要在核函数中加上同步了，因为16×16的线程块包含4个线程束，我们需要让所有线程束完成从全局内存到共享内存的读取之后（这里的关键是一个线程束会用到由其他线程束读取到共享内存的数据），再开始利用共享内存中的数据进行计算，还需要让当前k循环计算结束后才能进入下一个k循环做下一组全局内存到共享内存的读取。为此，加入两个同步：

            __syncthreads();
            // 计算矩阵A_ik和B_kj的乘
            for (int l = 0; l < Nsub; l++) {
                C_sub += A_sub[row + l * Nsub]
                         * B_sub[l + col * Nsub];
            }
            __syncthreads();

再次编译运行，结果如下：

./shared 4096
Generating Matrices took    5.777 ms.
Matrix multiplication (hipblasDgemm) took   28.971 ms.
Matrix multiplication (kernel) took  174.049 ms.
Max error: 0.0000e+000

可以看到速度提升了两倍多，再将Nsub定义为32（记得为核函数加上__launch_bounds__(1024)），结果如下：

./shared 4096
Generating Matrices took    5.597 ms.
Matrix multiplication (hipblasDgemm) took   28.669 ms.
Matrix multiplication (kernel) took  130.636 ms.
Max error: 0.0000e+00

相比shared.cc中Nsub的大小对核函数性能的显著影响，Nsub的大小对于dcu.cc中核函数性能的影响则小很多，可以尝试修改并测试。

3.5.4 让每个线程计算更多的矩阵元

当你想继续增加Nsub时，可能会想到32×32=1024已经达到单个线程块中线程数的极限了。但有了之前章节的经验，你应该知道解决的办法，即每个线程不必只计算子矩阵Cij的一个矩阵元，而是可以多算几个，这样子矩阵就可以比线程块更大了。比如说每个线程负责将子矩阵Aik和Bkj的两行两列的4个元素从全局内存读到共享内存，计算时负责从共享内存中的Aik中取两行，Bkj中取两列，算出子矩阵Cij的4个元素，这样就可用32×32的线程块处理64×64的子矩阵了。

不过64×6464×64的子矩阵又达到了共享内存的极限，因为此时A_sub和B_sub加起来所需的空间为2×64×64×8B=64K刚好为单个计算单元中共享内存的配额。此时如果仍想继续增加Nsub，需要将矩阵AA与矩阵BB重新分块，令Aik为Nsub×Nsub′，Bkj为Nsub′×Nsub的矩阵，这样矩阵C的分块仍保持不变（线程网格也保持不变），通过Nsub′可以调整所需的共享内存空间。

有了这些编程思路，我们可以进一步优化程序了，但为了让程序不显得太复杂，突出我们的优化重点，从本节开始我们会省去边界检查，只考虑Nsub和Nsub′都能整除N的情况。

作为练习，你可以尝试编写包含边界检查的版本，以处理NN任意取值的情况。类似我们在shared.cc的read_matrix中所写的，如果行列没有超过全局内存中矩阵M的边界则将元素读入共享内存中的子矩阵的适当位置，否则将子矩阵该位置的元素置零。但我们之前的写法是有优化空间的，因为只有当子矩阵包含矩阵的最后一行或列时才需要做边界检查，而不需要对所有子矩阵做边界检查。

现在开始优化，首先通过宏定义一些常数，包括Nsub、Nsub′（记为Nsub_）、线程块的维度BLOCKDIM_X和BLOCKDIM_Y以及每个线程需要处理的A的子矩阵的行数Nx、B的子矩阵的列数Ny：

#define Nsub  64
#define Nsub_ 64
#define BLOCKDIM_X 32
#define BLOCKDIM_Y 32
#define Nx Nsub / BLOCKDIM_X
#define Ny Nsub / BLOCKDIM_Y

比如上面的定义就表明我们想用32×32的线程块来处理矩阵A、B、C都划分成一系列64×64的子矩阵的情况，每个线程需要处理子矩阵Aik中的两行（Nx等于2）与子矩阵Bkj中的两列（Ny也等于2）的内积，算出子矩阵Cij的4个（Nx乘Ny）元素。

接下来是改写read_matrix函数，因为现在Aik和Bkj的行列数可以是不同的，此外着眼于性能，我们希望在编译时就把更多的值确定下来，因此使用函数模板：

template <int nrows, int ncols>
__device__
void read_matrix(double* M, int N, int i, int j,
                 double* M_ij, int LD)
{
    int row, col, offset, Row, Col, Offset;
    // 子矩阵行列数与线程块行列数的比值
    constexpr int Rrow = nrows / BLOCKDIM_X;
    constexpr int Rcol = ncols / BLOCKDIM_Y;
    for (int J = 0; J < Rcol; J++) {
        for (int I = 0; I < Rrow; I++) {
            // 矩阵元在子矩阵M_ij中的行列
            row = threadIdx.x + I * BLOCKDIM_X;
            col = threadIdx.y + J * BLOCKDIM_Y;
            offset = row + col * LD;
            // 矩阵元在矩阵M中的行列
            Row = i * nrows + row;
            Col = j * ncols + col;
            Offset = Row + Col * N;
            // 将矩阵元从全局内存读入共享内存
            M_ij[offset] = M[Offset];
        } // I
    } // J
}

在这一版的read_matrix中，每个线程总共需要读取的线程数由子矩阵行列数与线程块行列数的比值Rrow和Rcol来确定（R代表ratio），constexpr表示我们要求这两个值在编译时就确定下来。举例来说，当子矩阵的行数nrows和列数ncols都为64，线程块的行列数BLOCKDIM_X和BLOCKDIM_Y都为32时，Rrow和Rcol都等于2，因此对I和J的两重循环总共循环四次，所以每个线程读取了4个矩阵元。

将子矩阵读取到共享内存之后，下一步可以计算子矩阵的乘积了，由于现在每个线程需要计算多个矩阵元，相比之前的代码要复杂一些，因此我们也专门编写一个设备端函数multiply_submatrices：

__device__
void multiply_submatrices(double* A_sub, double* B_sub,
                          double* c)
{
    for (int k = 0; k < Nsub_; k++) {
        for (int J = 0; J < Ny; J ++) {
            for (int I = 0; I < Nx; I ++) {
                // 每个线程计算Nx * Ny个c
                int i = threadIdx.x + I * BLOCKDIM_X;
                int j = threadIdx.y + J * BLOCKDIM_Y;
                c[I + J * Nx] += A_sub[i + k * Nsub]
                                *B_sub[k + j * Nsub_];
            } // I
        } // J
    } // k
}

在这个函数中，我们让每个线程负责Nx行，Ny列，行与行之间间隔BLOCKDIM_X，列与列间隔BLOCKDIM_Y，计算出Nx乘Ny个C的元素，保存在局部内存中的数组c中。

类似的，再编写一个用于将矩阵元写回矩阵C中的函数write_matrix：

__device__
void write_matrix(double* C, int N, int i, int j, double* c)
{
    int row, col, Row, Col;
    for (int J = 0; J < Ny; J++) {
        for (int I = 0; I < Nx; I++) {
            row = threadIdx.x + I * BLOCKDIM_X;
            col = threadIdx.y + J * BLOCKDIM_Y;
            Row = i * Nsub + row;
            Col = j * Nsub + col;
            if ((Row < N) && (Col < N)) {
                C[Row + Col * N] = c[I + J * Nx];
            }
        }
    }
}

此时的核函数matrix_multiply简单明了：

__global__ __launch_bounds__(BLOCKDIM_X * BLOCKDIM_Y)
void matrix_multiply(int N, int Nblk, double *A, double *B,
                     double *C)
{
    // 共享内存中的A_sub,B_sub用于保存A,B的子矩阵
    __shared__ double A_sub[Nsub * Nsub_];
    __shared__ double B_sub[Nsub_ * Nsub];
    // 每个线程块负责计算子矩阵C_ij
    for (int i = blockIdx.x; i < Nblk; i += gridDim.x) {
    for (int j = blockIdx.y; j < Nblk; j += gridDim.y) {
        // 每个线程负责计算子矩阵C_ij中Nx*Ny个矩阵元
        double c[Nx * Ny];
        for (int k = 0; k < Nx * Ny; k++) {
            c[k] = 0.0;
        }
        // 假设Nsub_能整除N
        for (int k = 0; k < N / Nsub_; k++) {
            // 将子矩阵A_ik从全局内存读到共享内存
            read_matrix<Nsub, Nsub_>(A, N, i, k, A_sub, Nsub);
            // 将子矩阵B_kj从全局内存读到共享内存
            read_matrix<Nsub_, Nsub>(B, N, k, j, B_sub, Nsub_);
            __syncthreads();
            // 计算矩阵A_ik和B_kj的乘
            multiply_submatrices(A_sub, B_sub, c);
            __syncthreads();
        } // k
        // 将矩阵元写到C中
        write_matrix(C, N, i, j, c);
    } // j
    } // i
}

流程与之前版本是一致的，区别主要有三点：A_sub与B_sub的大小不同；需要一个局部内存中的数组c用于存储每个线程负责的多个矩阵元；k的循环范围。注意将主机代码中的线程块配置改为：

    dim3 blockSize = dim3(BLOCKDIM_X, BLOCKDIM_Y, 1);

编译运行如下（more_elements.cc，Nsub和Nsub_都为64，BLOCKDIM_X和BLOCKDIM_Y都为32）：

hipcc more_elements.cc -o more_elements -I$ROCM_PATH/rocrand/include -L$ROCM_PATH/lib -lhipblas -lhiprand
./more_elements 4096
Generating Matrices took    5.634 ms.
Matrix multiplication (hipblasDgemm) took   29.711 ms.
Matrix multiplication (kernel) took   72.035 ms.
Max error: 0.0000e+00
./more_elements 16384
Generating Matrices took   11.246 ms.
Matrix multiplication (hipblasDgemm) took 1734.654 ms.
Matrix multiplication (kernel) took 4598.284 ms.
Max error: 4.7294e-11

调整一下这些规模，当Nsub为128，Nsub_为16，BLOCKDIM_X和BLOCKDIM_Y都为16时：

./more_elements 4096
Generating Matrices took    5.592 ms.
Matrix multiplication (hipblasDgemm) took   28.318 ms.
Matrix multiplication (kernel) took   23.790 ms.
Max error: 0.0000e+00
./more_elements 16384
Generating Matrices took   11.211 ms.
Matrix multiplication (hipblasDgemm) took 1734.478 ms.
Matrix multiplication (kernel) took 1515.077 ms.
Max error: 5.1841e-11

你可能会觉得我们编写的核函数比hipBLAS库的性能更好了。只能说在目前版本（DTK-22.10）中，针对本节中的简单情况（又未做边界检查），是这样的。

你可能想，我们继续增加Nsub岂不更好，Nsub等于256，Nsub_等于16对于共享内存容量来说是允许的，但这样修改之后（线程块保持16×16）运行，会出现这样的结果：

Generating Matrices took    5.635 ms.
Matrix multiplication (hipblasDgemm) took   29.254 ms.
Matrix multiplication (kernel) took 2057.917 ms.
Max error: 0.0000e+00

这里性能骤降的原因在于我们碰到了又一条限制，即寄存器的数量是有限的，在3.2.2节中提到过，每个SIMD单元有64KB的向量寄存器配额，一个计算单元总计256KB，每个寄存器是32位的（4B），即所有线程总共可以使用64K个向量寄存器，对于大小为1024的线程块，这意味着每个线程最多只能使用64个向量寄存器；对于我们这里大小为256的线程块（也包括更小的线程块），每个线程最多只能使用256个向量寄存器。但注意到此时局部内存中的数组c的长度是256了（Nx和Ny都是16），考虑到数组元素是双精度浮点数来说，仅这个数组就需要512个寄存器，这是不允许的，编译器会将这个数组分配到片外内存中去，结果函数multiply_submatrices变成了反复对全局内存而非局部内存做累加，如此性能下降就是意料之中的了。

3.5.5 局部内存中的数组

本节我们详细分析一下上节中的程序more_elements.cc，取性能最优的配置（Nsub为128，Nsub_为16，BLOCKDIM_X和BLOCKDIM_Y都为16）。

首先，你可能会对函数read_matrix和multiply_submatrices中循环的写法有疑问，按照之前章节中的写法，这两个函数中的循环应该是这样的：

template <int nrows, int ncols>
__device__
void read_matrix(double* M, int N, int i, int j,
                 double* M_ij, int LD)
{
    int row, col, offset, Row, Col, Offset;
    for (col = threadIdx.y; col < ncols; col += blockDim.y) {
    for (row = threadIdx.x; row < nrows; row += blockDim.x) {
            offset = row + col * LD;
            Row = i * nrows + row;
            Col = j * ncols + col;
            Offset = Row + Col * N;
            M_ij[offset] = M[Offset];
        } // row
    } // col
}

以及这样的：

__device__
void multiply_submatrices(double* A_sub, double* B_sub,
                          double* c)
{
    for (int k = 0; k < Nsub_; k++) {
        for (int j = threadIdx.y; j < Nsub; j+= blockDim.y) {
        for (int i = threadIdx.x; i < Nsub; i+= blockDim.x) {
            int I = i / blockDim.x;
            int J = j / blockDim.y;
            c[I + J * Nx] += A_sub[i + k * Nsub]
                            *B_sub[k + j * Nsub_];
        } // i
        } // j
    } // k
}

我们按上面的写法来修改read_matrix并编译运行：

./more_elements 4096
Generating Matrices took    5.624 ms.
Matrix multiplication (hipblasDgemm) took   28.977 ms.
Matrix multiplication (kernel) took   35.649 ms.
Max error: 0.0000e+00

结果没问题，但性能下降了。再按上面的写法修改multiply_submatrices，运行结果如下：

./more_elements 4096
Generating Matrices took    5.680 ms.
Matrix multiplication (hipblasDgemm) took   28.496 ms.
Matrix multiplication (kernel) took 2097.088 ms.
Max error: 0.0000e+00

这次导致了很严重的性能下降，原因是什么呢？

答案是编译器不知道该如何为局部内存中的数组c的元素分配寄存器，所以只能将它放到片外了。仔细思考一下，因为我们是使用C语言编程，数组实际上是一块连续的内存，通过指针来访问数组中的元素，但寄存器是不同于其他内存的，不存在“寄存器的地址”这个概念，如果你熟悉X86架构的CPU，你可能知道每个寄存器都有特殊的名称；对于DCU的向量寄存器，每个都有特定的编号，但这种编号只可以理解为代号，不能理解为索引，不存在某个寄存器保存着其他寄存器的编号这种概念。

从这个角度来讲，本质上是不可以把数组放到寄存器上的。上一节中的写法之所以能让编译器把数组元素分配到一个个寄存器上，是因为编译器在编译时就能解析语句的索引，知道c[I + J * Nx]具体是c的哪一个元素。更通俗的讲，编译器只能将“假”的数组放到寄存器上，数组中每一个元素就是一个的独立变量，不具备根据索引访问元素的功能。

具体来说，按上节中的写法，编译器在编译时可以做循环展开，可以将I从0到7，J从0到7直接代入，等价于我们写类似这样的代码（将其他的宏也替换掉）：

    // I = 0, J = 0
    int i = threadIdx.x;
    int j = threadIdx.y;
    c[0] += A_sub[i + k * 128] * B_sub[k + j * 16];
    // I = 1, J = 0
    int i = threadIdx.x + 16;
    int j = threadIdx.y;
    c[1] += A_sub[i + k * 128] * B_sub[k + j * 16];
    // I = 2, J = 0
    int i = threadIdx.x + 32;
    int j = threadIdx.y;
    c[2] += A_sub[i + k * 128] * B_sub[k + j * 16];
    // ...

此时，我们用三个变量c0、c1、c2替换上面的c[0]、c[1]和c[2]是没有任何问题的。编译器可以按给c0、c1、c2分配寄存器的方式，给c[0]、c[1]、c[2]分配寄存器。但如果按本节的写法，虽然逻辑上看没有任何问题，可是在编译时，编译器是无法做循环展开的。i和j的初始值分别是threadIdx.x 和threadIdx.y，这是在运行时才能确定的，而且对每个线程各不相同；i和j在编译时无法确定，I和J是多少就更不知道了。对于

c[I + J * 8] += A_sub[i + k * 128] * B_sub[k + j * 16];

这样的代码，应该用c?来替换掉c[I + J * 8]，这在编译时是不知道的，因此数组c的元素没法分配到寄存器上。

使用局部内存中的数组是有好处的，可以让我们的代码简洁而灵活，设想一下，对于我们目前的配置，每个线程需要计算64个矩阵元，如果一句句的写：

    c0  += ...;
    c1  += ...;
    c2  += ...;
    ...
    c63 += ...;

代码会非常繁杂，而且随着Nx和Ny的改变，变量的数目还会发生变化。但按照上节的写法，使用局部内存中的数组加上循环展开，简单的几行代码就可以达到效果。

最后我们顺带提一下循环展开，在hipcc默认-O3的优化级别下，大部分情况下编译器会自动将可以展开的循环进行适当的展开，如果需要手动指定，则在循环之上添加编译器指导语句#pragma unroll，或用#pragma nounroll要求编译器不做展开。循环展开除了能让我们简洁的使用局部内存的数组之外，还有很多优点，比如说不必反复判断循环条件，根据条件进行跳转；展开后的代码也更有利于编译器进行指令级的优化。我们在本书中不想涉及汇编层次的机器级代码，因此对这部分也不做更详细的介绍了，值得一提的是循环展开对算法上有一定要求，比如对于本节中的矩阵乘法问题，由于算法的仔细设计，循环次数在编译时就完全给定了；对于一些一般性的算法，循环指标的初始化、循环次数等在运行时才能确定，就难以利用循环展开进行优化了。

3.5.6 共享内存的bank冲突

在本节的最后部分，我们来讨论共享内存的访存问题。在3.2.2节，我们提到过共享内存在硬件上对应每个计算单元上容量为64KB的局域数据共享存储（方便起见，下文直接用共享内存来称呼它），它有32个bank，每个bank有512个32位的条目。每个bank在一个时钟周期内可以做一次32位的读或写，32个bank可以同时进行，因此共享内存有很高的带宽，此外共享内存上还有专门处理32位整型原子加的硬件单元。

DCU的一个线程束有64个线程，因此如果线程束访问共享内存中连续的64个4字节数据，则每半个线程束中的32个线程刚好与共享内存中的32个bank一一对应，在两个时钟周期就可以完成访存请求，这种情况下对共享内存的访存是最高效的。但如果64个线程都访问一个bank中的不同条目时，则是最低效的情况，这种情况下会产生严重的bank冲突，访存只能串行进行，需要64个时钟周期才能完成。因此在使用共享内存时，需要注意数据在共享内存中的布局和线程的访存方式，避免多个线程在一个时钟周期内访问同一bank中的不同条目（多个线程访问同一bank中的同一条目不会带来bank冲突，因为共享内存支持广播方式的访存），这样才能充分利用共享内存的带宽。

下面仍以我们的函数multiply_submatrices为例，分析一下子矩阵大小以及线程块大小对于共享内存访存效率的影响：

__device__
void multiply_submatrices(double* A_sub, double* B_sub,
                          double* c)
{
    for (int k = 0; k < Nsub_; k++) {
        for (int J = 0; J < Ny; J ++) {
            for (int I = 0; I < Nx; I ++) {
                // 每个线程计算Nx * Ny个c
                int i = threadIdx.x + I * BLOCKDIM_X;
                int j = threadIdx.y + J * BLOCKDIM_Y;
                c[I + J * Nx] += A_sub[i + k * Nsub]
                                *B_sub[k + j * Nsub_];
            } // I
        } // J
    } // k
}

子矩阵A_sub和B_sub是共享内存中的两个数组，由于是双精度浮点数（64位），因此每一个元素占用连续的两个bank。对于16×16的线程块，刚好在一个周期内threadIdx.x从0到15的16个线程可以从32个bank中读取A_sub的16个元素，不存在bank冲突，最大化的利用了共享内存的带宽；对于B_sub，也是16个线程访问同一个元素，属于广播方式的访存，同样不存在bank冲突。

再分析一种情况，使用8×88×8的线程块，Nsub、Nsub_也都等于8时，这与我们使用共享内存优化的最初版本的shared.cc中Nsub等于8时的程序基本一致。可以修改Nsub、Nsub_、BLOCKDIM_X和BLOCKDIM_Y后重新编译，测试一下这种情况下的性能：

./more_elements 4096
Generating Matrices took    5.654 ms.
Matrix multiplication (hipblasDgemm) took   28.996 ms.
Matrix multiplication (kernel) took  428.407 ms.
Max error: 0.0000e+00

在这种情况下，对于读取A_sub的8个元素来说也不存在bank冲突，因为这8个元素保存在连续的16个bank中，同行不同列（threadIdx.x相同、threadIdx.y不同）的线程通过广播方式读取同一元素；对于B_sub，同一列（threadIdx.y相同）的8个线程读取相同的元素，不存在bank冲突，但对于不同列读取的元素则存在bank冲突。考虑threadIdx.y等于0到3的四列32个线程，分别读取B_sub[k]、B_sub[k+8]、B_sub[k+16]和B_sub[k+24]四个元素，因为每16个元素就会占用全部的32个bank，所以B_sub[k]与B_sub[k+16]是同一bank中的不同元素，B_sub[k+8]和B_sub[k+24]也是同一bank中的不同元素。如何解决这种bank冲突呢？

一种简单的处理方式是将B_sub的空间多申请一行，即为B_sub申请9行8列的空间，但只使用8行8列：

__shared__ double B_sub[(Nsub_ + 1) * Nsub];

修改read_matrix的参数LD：

// 将子矩阵B_kj从全局内存读到共享内存
read_matrix<Nsub_, Nsub>(B, N, k, j, B_sub, Nsub_ + 1);

并修改multiply_submatrices中对于B_sub的引用：

c[I + J * Nx] += A_sub[i + k * Nsub]
                *B_sub[k + j * (Nsub_ + 1)];

分析一下，此时threadIdx.y等于0到3的四列32个线程，分别读取B_sub[k]、B_sub[k+9]、B_sub[k+18]和B_sub[k+27]四个元素，均在不同的bank中，就不会冲突了。编译运行，结果如下：

./more_elements 4096
Generating Matrices took    5.628 ms.
Matrix multiplication (hipblasDgemm) took   28.953 ms.
Matrix multiplication (kernel) took  355.662 ms.

练习：分析其他线程块规模与子矩阵规模下的共享内存bank冲突情况。

最后还要补充一点，一个线程在一个周期内最多产生两个4字节（64位）的访存请求，如果访存请求不满4字节（比如只读取一个半精度浮点数，FP16），仍然会占用一个bank的32位带宽；如果是128位的访存请求，比如读取两个双精度浮点数、双精度复数等，也只能用到共享内存一半的带宽。举例来说，假设线程ID为0的线程需要读取128位的数据在0~3号bank中，线程ID为1的线程需要的数据在4~7号bank中，这些访存请求需要拆到两个时钟周期，第一个周期线程ID为0的线程读取到0~1号bank的数据，线程ID为1的线程读取到4~5号bank的数据，此时2~3号、6~7号bank没有访存请求，因此有一半的带宽没有用上，第二个周期同理，只有2~3号、6~7号bank在处理访存请求，0~1号、4~5号bank没有用上。

3.5.7 小节

在本节中，我们以双精度矩阵乘法为例进一步介绍了DCU编程与优化方法。我们首先使用hipRAND库随机生成了两个矩阵，利用hipBLAS库计算了矩阵乘，以此作为在DCU程序中使用库的简单范例，并对DCU在处理这些问题时的性能有了初步的印象。之后我们自己编写了矩阵乘法的程序，通过分析，发现矩阵乘法中存在大量的数据复用。抓住这一点，我们尽可能地利用共享内存进行优化，在维持代码相对简洁可读的情况下，达到了不错的性能。

在优化的过程中，我们也看到共享内存的大小，共享内存的访存方式，局部内存的有效利用，还有程序的编写方法对编译器的影响都会对程序的性能产生影响，根据算法的不同，这些因素的影响或多或少，有时还会彼此关联、彼此制约，导致问题很复杂。本节也仅仅是抛砖引玉，希望读者参考本节的优化案例，在编程解决自己的问题时，根据算法特点，首先抓住少量影响性能的主要因素进行优化，再处理其他因素。

思考：在本节的程序中，有哪些硬件资源的使用情况有关联。

举例来说：当核函数使用64KB的共享内存时，每个计算单元最多并发多少个线程束。

补充一点，本节中申请共享内存的大小在编译时就确定了（Nsub、Nsub_都通过宏定义替换为具体的数），这种方式称为静态共享内存。有时共享内存的大小可能难以在编译时确定，此时还可以动态申请共享内存，在调用核函数时通过三个尖括号中的第三个参数来指定，例如：

    int shmemSize = 2 * Nsub * Nsub_ * sizeof(double);
    matrix_multiply<<<numBlocks, blockSize, shmemSize>>>
                    (N, Nblk, dA, dB, dC);

此时Nsub和Nsub_可以是变量。但在核函数中使用共享内存会麻烦一些：

    // 获取动态共享内存的指针
    extern __shared__ float shmem[];
    // 声明A_sub,B_sub用于保存A,B的子矩阵
    double* A_sub = (double*)shmem;
    double* B_sub = (double*)&shmem[2 * Nsub * Nsub_];

需要通过指针手动实现。