CUDA C编程（三十四）多GPU上的有限差分

接下来通过使用有限差分的方法求解二维波动方程，将会学习到如何跨设备重叠计算和通信。

二 次 波 动 方 程 的 模 板 计 算

二维波的传播由以下波动方程来决定：

其中,u(x,y,t)是波场，v(x,y)是介质的速度。这是一个二阶偏微分方程。求解这种偏微分方程的典型方法是使用规则的笛卡尔网格上的有限差分法。更简单地说，有限差分法近似于使用一个模板求导以计算规则网格中单一点的倒数，具体方法是围绕该点的多个局部点应用一个函数。下图展示了17点的模板。如果要求解中心点的导数，则需要使用16个离中心点最近的局部点。

按计算方式，偏导数可以由一个泰勒展开式来表示，其在一个维度的实现用以下伪代码来表示。这个伪代码使用一维数据der_u来从当前元素u[i]前面的4个元素(u[i+d])和4个后面的元素(u[i-d])中积累贡献，c数组存储导数稀疏，u[i]是计算的中心点，der_u[i]是得到的中心点的导数。

der_u[i] = c[0] * u[i];
for(int d = 1; d <= 4; d++)
   der_u[i] += c[d] * (u[i-d] + u[i+d]);

多 GPU 程 序 的 典 型 模 式
为了准确地模拟通过不同戒指的波传播，需要大量的数据。但单GPU的全局内存没有足够的空间存储模拟过程的状态。这就需要跨多个GPU的数据域分解。假设在二维数组中z轴是最内层的维度，那么可以沿y轴分割数据使其分布在多个GPU上。因为对一个给定点的计算需要它两侧最近的4个点，所以需要为存储在每个GPU上的数据添加填充区域，如下图所示。此填充区域或halo区域，使得在波传播计算的每次循环中，相邻设备之间进行数据交换。

下图所示的域分解，使用了多GPU来求解波动方程，在模拟的每一步时都会使用以下的通用模式：

1. 在一个流中使用相邻的GPU计算halo区域和交换halo数据；
2. 在不同的流中计算内部区域；
3. 在进行下一次循环之前，在所有设备上进行同步计算。

如果使用两个不同的流，一个halo计算和通信，另一个用于内部区域的计算，步骤1可以用步骤2重叠。如果内部计算所需的计算时间比halo操作所需的时间长，可以通过使用多个GPU隐藏halo通信的性能影响来实现线程加速。在两个GPU上进行模板计算的伪代码如下：

for(int istep = 0; istep < nsteps; istep++){
   for(int i = 0; i < 2; i++){
      cudaSetDevice(i);
      2dfd_kernel<<<..., stream_halo[i]>>>(...);
   }
    
   cudaMemcpyAsync(..., cudaMemcpyDeviceToDevice, stream_halo[0]);
   cudaMemcpyAsync(..., cudaMemcpyDeviceToDevice, stream_halo[1]);
   
   for(int i = 0; i < 2; i++){
      cudaSetDevice(i);
      2dfd_kernel<<<..., stream_internal[i]>>>(...);
   }
    
   for(int i = 0; i < 2; i++){
      cudaSetDevice(i);
      cudaDeviceSynchronize();
   }
}

不需要为复制操作设置当前设备，但是在启动内核之前，必须指定当前设备。

多 GPU 上 的 二 维 模 板 计 算
在二维模板计算中使用两个设备数组，一个保存当前的波场，另一个保存更新的波场。如果将x定义为最内层的数组维度，并且将y作为最外层的维度，那么可以沿着y轴跨设备均匀地分配计算。

因为更新一个点需要访问9个最近点，所以很多点都将共享输入数据。因此，使用共享内存可以减少全局内存的访问。共享内存地使用量等同于保存相邻数据的线程块大小，该线程块被8个点填充，如下图所示，代码如下：__shared__ float line[4 + BDIMX + 4];用于存储y轴模板值的9个浮点值被声明为一个核函数的本地数组，并因此存储在寄存器中。当沿y轴在当前元素的前后加载元素时，用到的寄存器很像用来减少冗余访问的共享内存。

下图说明了单线程中沿x轴存储模板值的共享内存和沿y轴存储模板值的9个寄存器：

一旦输入数据被分配并初始化，在每个GPU线程上实现有限差分(FD)的模板计算可以写成如下代码：

float tmp = coef[0] * line[stx] * 2.0f;

for(int d = 1; d <= 4; d++){
   tmp += coef[d] * (line[stx - d] + line[stx + d]);
}

for(int d = 1; d <= 4; d++){
   tmp += coef[d] * (yval[stx - d] + yval[stx + d]);
}

g_u1[idx] = 2.0f * current - g_u1[idx] + alpha * tmp;

二维模板计算的完整内核代码如下所示：

__global__ void kernel_2dfd(float *g_u1, float *g_u2, const int nx,
                            const int iStart, const int iEnd)
{
    // global to line index
    unsigned int ix  = blockIdx.x * blockDim.x + threadIdx.x;

    // smem idx for current point
    unsigned int stx = threadIdx.x + NPAD;
    unsigned int idx  = ix + iStart * nx;

    // shared memory for x dimension
    __shared__ float line[BDIMX + NPAD2];

    // a coefficient related to physical properties
    const float alpha = 0.12f;

    // register for y value
    float yval[9];

    for (int i = 0; i < 8; i++) yval[i] = g_u2[idx + (i - 4) * nx];

    // skip for the bottom most y value
    int iskip = NPAD * nx;

#pragma unroll 9
    for (int iy = iStart; iy < iEnd; iy++)
    {
        // get yval[8] here
        yval[8] = g_u2[idx + iskip];

        // read halo part
        if(threadIdx.x < NPAD)
        {
            line[threadIdx.x]  = g_u2[idx - NPAD];
            line[stx + BDIMX]    = g_u2[idx + BDIMX];
        }

        line[stx] = yval[4];
        __syncthreads();

        // 8rd fd operator
        if ( (ix >= NPAD) && (ix < nx - NPAD) )
        {
            // center point
            float tmp = coef[0] * line[stx] * 2.0f;

#pragma unroll
            for(int d = 1; d <= 4; d++)
            {
                tmp += coef[d] * ( line[stx - d] + line[stx + d]);
            }

#pragma unroll
            for(int d = 1; d <= 4; d++)
            {
                tmp += coef[d] * (yval[4 - d] + yval[4 + d]);
            }

            // time dimension
            g_u1[idx] = yval[4] + yval[4] - g_u1[idx] + alpha * tmp;
        }

#pragma unroll 8
        for (int i = 0; i < 8 ; i++)
        {
            yval[i] = yval[i + 1];
        }

        // advancd on global idx
        idx  += nx;
        __syncthreads();
    }
}

重 叠 计 算 与 通 信

此二维模板的执行配置使用一个具有一维线程块的一维网格，在主机上的声明如下所示：

dim3 block(BDIMX);
dim3 grid(nx / block.x);

在主机上波的传播时间通过使用nsteps次迭代的时间循环来控制。在第一次时间步骤中，核函数kernel_add_wavelet引入GPU0介质中一个扰动。随着时间的变化进一步的迭代传递干扰。因为halo区域计算和数据交换被安排在每个设备的stream_halo流中，内部区域的计算被安排在每个设备的stream_internal流中，所以在此二维模板上计算和通信可以重叠。

// add wavelet only onto gpu0
        if (istep == 0)
        {
            CHECK(cudaSetDevice(0));
            kernel_add_wavelet<<<grid, block>>>(d_u2[0], 20.0, nx, iny, ngpus);
        }

        // halo part
        for (int i = 0; i < ngpus; i++)
        {
            CHECK(cudaSetDevice(i));

            // compute halo
            kernel_2dfd<<<grid, block, 0, stream_halo[i]>>>(d_u1[i], d_u2[i],
                    nx, haloStart[i], haloEnd[i]);

            // compute internal
            kernel_2dfd<<<grid, block, 0, stream_body[i]>>>(d_u1[i], d_u2[i],
                    nx, bodyStart[i], bodyEnd[i]);
        }

        // exchange halo
        if (ngpus > 1)
        {
            CHECK(cudaMemcpyAsync(d_u1[1] + dst_skip[0], d_u1[0] + src_skip[0],
                        iexchange, cudaMemcpyDefault, stream_halo[0]));
            CHECK(cudaMemcpyAsync(d_u1[0] + dst_skip[1], d_u1[1] + src_skip[1],
                        iexchange, cudaMemcpyDefault, stream_halo[1]));
        }

        for (int i = 0; i < ngpus; i++)
        {
            CHECK(cudaSetDevice(i));
            CHECK(cudaDeviceSynchronize());

            float *tmpu0 = d_u1[i];
            d_u1[i] = d_u2[i];
            d_u2[i] = tmpu0;
        }
    }

编 译 和 执 行

从Wrox.com中可以下载包含完整示例文件simple2DFD.cu。通过以下代码进行编译：

$ nvcc -arch=sm_20 -O3 -use_fast_math simple2DFD.cu simple2DFD

以下时系统输出的例子，该系统有两个M2090设备：

使用的性能用Mcells/s表示。对于二维情况，其定义如下：

只用一个设备运行时，产生以下结果：

simple2DFD显示了从一个设备移动到两个设备的近似线形扩展(有效率为96%)，由此，可以得出结论，转移halo区域增加的通信开销可以采用CUDA流在多个GPU上有效地隐藏。

在simple2DFD中检查内核地资源使用情况也是很有趣的，在nvcc编译器中添加标志可以报告内核资源使用情况：

$ nvcc -arch=sm_20 -Xptxas -v simple2DFD.cu -o simple2DFD

nvcc编译器输出一下信息：

上面代码输出的最后一行表明，核函数kernel_2dfd使用了26个寄存器，160个字节的共享内存和每个线程上的一些常量内存(用来存储核函数的一些参数)。