wjin

CUDA教程10 -- 多GPU编程

 7519  16 分钟

§多GPU编程模型

在超级计算机上，一个系统可能安装有多个GPU，如何利用这些GPU协同工作需要我们在代码上作出一定处理。GPU属于计算设备，本身与操作系统内程序的执行流程互不相干。无论程序以单线程运行、多线程运行，亦或以多进程运行，甚至跨操作系统协作，GPU对于程序而言都只是一个外设，CUDA提供的接口使我们可以方便地操纵GPU而已。

典型的多GPU编程模型有单线程对应多GPU、多线程对应多GPU、多进程对应多GPU，这几种模型之间的区别我们将一一讲述。在实践中，选择哪种模型要看计算机的类型，如果是多机系统就必须选择多进程模型，如果是单机系统应当看调用者使用哪种模型。

§cuda流

CUDA在执行GPU代码时是以任务队列的形式异步执行的，这个任务队列在CUDA中被称为cuda流。每当我们通过启动核函数的语法在GPU上执行一个函数时，仅仅是将任务提交给了任务队列，函数不一定是立即执行的，程序的执行流程也会立即返回，继续执行后续代码，不会等待核函数执行完成。

在CUDA中，任务队列的数据结构是cudaStream，这个类型的细节无须关心，我们只使用其指针类型cudaStream_t来操作任务队列。和cudaStream_t相关的常用接口有：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

// 创建一个cuda流
cudaError_t cudaStreamCreate(cudaStream_t *pStream);

// 销毁一个cuda流
cudaError_t cudaStreamDestroy(cudaStream_t stream);

// 等待一个cuda流中的任务全部完成
cudaError_t cudaStreamSynchronize(cudaStream_t stream);

// 等待cuda流中的某个事件的发生，flags必须为0
cudaError_t cudaStreamWaitEvent(cudaStream_t stream, cudaEvent_t event,
                                unsigned int  flags);

一个cuda流就是一个任务队列，因此通过创建多个cuda流就可以将GPU的执行流程分成互不干扰的多个流程，使不同的GPU同时执行任务。

§单线程和多线程模型

单线程操作多个GPU通过设置cuda流，然后在循环中依次操作不同的GPU来实现。例如一个向量加法的实现：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

// 查询GPU的数量
int num_devices;
CHECK(cudaGetDeviceCount(&num_devices));

// 创建cuda流
std::vector<cudaStream_t> streams(num_devices);
for (int i = 0; i < num_devices; ++i) {
    CHECK(cudaSetDevice(i));   // 设置当前使用的GPU
    CHECK(cudaStreamCreate(&streams[i]));
}

// 数据
constexpr int N;
std::vector<double> A(N);
std::vector<double> B(N);

// 传输到GPU上
// 每个GPU上的数据个数
std::vector<int> length(num_devices);
for (int i = 0; i < N % num_devices; ++i)
    length[i] = N / num_devices + 1;
for (int i = N % num_devices; i < num_devices; ++i)
    length[i] = N / num_devices;

// 每个GPU上第一个数据的偏移
std::vector<int> offset(num_devices + 1);
for (int i = 0, s = 0; i < num_devices; ++i) {
    offset[i] = s;
    s += length[i];
}
offset.back() = N;

// GPU上的向量只有一部分
std::vector<double *> dev_A(num_devices, nullptr);
std::vector<double *> dev_B(num_devices, nullptr);
for (int i = 0; i < num_devices; ++i) {
    CHECK(cudaSetDevice(i));      // 设置GPU
    CHECK(cudaMalloc(&dev_B[i], length[i] * sizeof(double)));  // 分配显存
    CHECK(cudaMalloc(&dev_B[i], length[i] * sizeof(double)));  // 分配显存
    // CHECK(cudaMemcpy(&dev_A[i], &A[offset[i]],
    //                 length[i] * sizeof(double),
    //                 cudaMemcpyHostToDevice));   // 复制数据
    // CHECK(cudaMemcpy(&dev_B[i], &B[offset[i]],
    //                 length[i] * sizeof(double),
    //                 cudaMemcpyHostToDevice));   // 复制数据
}

// 复制数据不需要设置需要使用的GPU，因为不同GPU的编址空间不相交，与内存也不相交，
// CUDA是可以从指针本身判断出它指针的显存位于哪个GPU的
for (int i = 0; i < num_devices; ++i) {
    CHECK(cudaMemcpy(&dev_A[i], &A[offset[i]],
                    length[i] * sizeof(double),
                    cudaMemcpyHostToDevice));   // 复制数据
    CHECK(cudaMemcpy(&dev_B[i], &B[offset[i]],
                    length[i] * sizeof(double),
                    cudaMemcpyHostToDevice));   // 复制数据
}

// 计算C = A + B
std::vector<double *> dev_C(num_devices, nullptr);
for (int i = 0; i < num_devices; ++i) {
    CHECK(cudaSetDevice(i));      // 设置GPU
    CHECK(cudaMalloc(&dev_C[i], length[i] * sizeof(double)));  // 分配显存
}

for (int i = 0; i < num_devices; ++i) {
    CHECK(cudaSetDevice(i));
    AddKernel<<<grid, block, 0, streams[i]>>>(   // 不同GPU在不同的cuda流上执行
        dev_C[i], dev_A[i], dev_B[i], length[i]);
}

// 销毁cuda流
for (int i = 0; i < num_devices; ++i) {
    CHECK(cudaStreamDestroy(streams[i]));
    CHECK(cudaFree(dev_A[i]));
    CHECK(cudaFree(dev_B[i]));
    CHECK(cudaFree(dev_C[i]));
}

其中的AddKernel和单GPU的完全一样：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

__global__ void AddKernel(
        double *result,
        const double *v1,
        const double *v2,
        const int N
) {
    unsigned int i = threadIdx.x + blockIdx.x * blockDim.x;
    unsigned int stride = gridDim.x * blockDim.x;
    for (; i < N; i += stride)
        result[i] = v1[i] + v2[i];
}

多线程的方法和单线程类似，在实践中，一般需要执行很多个核函数，因此最好采用线程池或者openmp、tbb等技术来实现，以降低创建线程的开销。多线程模型中也需要创建多个cuda流来实现GPU间的并行，向量加法的多线程多GPU实现留作思考问题。

§异步传输

前面说过，核函数的执行与CPU的执行流是异步的，具体表现在<<<>>>的调用将立即返回，任务进入cuda流排队完成。那么分配显存和数据传输是不是异步的呢？答案是否定的，并且在进行数据传输时，为了保证被传输的数据已经正确无误的生成在了指针指向的内存里面，cuda将对设备进行同步，其效果等价于调用cudaDeviceSynchronize，这显然妨碍了我们用计算掩盖数据传输的手段。为此，cuda引入了异步传输接口cudaMemcpyAsync，其原型如下：

1
2
3

cudaError_t cudaMemcpyAsync(void *dst, const void *src,
                            size_t count, cudaMemcpyKind kind,
                            cudaStream_t stream = 0);

需要注意，不是所有内存都能用于异步传输的，因为现代计算机的内存都使用请求调页系统，程序中的内存地址都是虚拟地址，运行时通过MMU翻译成物理地址，如果物理地址对应的页面不在内存中，则发出缺页异常，由操作系统负责将其调入内存再继续运行程序，而长久不用的内存则会被替换到磁盘上的swap文件或者swap分区。这就导致一个问题，如何保证异步传输时指针所指的内存页面始终有效呢？操作系统一般都提供了分配锁页内存的接口，以保证分配的内存始终在物理内存中不被换出。cuda为我们提供了cudaMallocHost和cudaFreeHost为我们屏蔽了操作系统间的差异，其原型为：

1
2
3
4
5

// 分配锁页内存
cudaError_t cudaMallocHost(void **ptr, size_t size);

// 释放锁页内存
cudaError_t cudaFreeHost(void *ptr);

使用示例：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

// 创建cuda流

// 数据
double *A = nullptr;
double *B = nullptr;
double *C = nullptr;

// 分配内存
CHECK(cudaMallocHost(&A, N * sizeof(double)));
CHECK(cudaMallocHost(&B, N * sizeof(double)));
CHECK(cudaMallocHost(&C, N * sizeof(double)));

// 生成数据
// 计算length、offset，分配显存

// 传输并计算
for (int i = 0; i < num_devices; ++i) {
    CHECK(cudaSetDevice(i));
    
    CHECK(cudaMemcpyAsync(&dev_A[i], &A[offset[i]],
                          length[i] * sizeof(double),
                          cudaMemcpyHostToDevice, streams[i]));   // 异步复制数据
    CHECK(cudaMemcpyAsync(&dev_B[i], &B[offset[i]],
                          length[i] * sizeof(double),
                          cudaMemcpyHostToDevice, streams[i]));   // 异步复制数据
    
    AddKernel<<<grid, block, 0, streams[i]>>>(   // 不同GPU在不同的cuda流上执行
        dev_C[i], dev_A[i], dev_B[i], length[i]);
    
    CHECK(cudaMemcpyAsync(&C[offset[i]], &dev_C[i],
                          length[i] * sizeof(double),
                          cudaMemcpyHostToDevice, streams[i]));   // 异步传输回来
}

// 释放内存
CHECK(cudaFreeHost(A));
CHECK(cudaFreeHost(B));
CHECK(cudaFreeHost(C));

// 销毁cuda流

值得注意是的显存有虚拟内存但没有请求调页机制，因此在显存间传输数据可以直接使用异步接口。

§多进程模型和CUDA-aware MPI

多进程模型就简单多了，由于进程间的数据是完全隔离的，因此不需要创建cuda流，和单GPU编程基本没有区别，但有个问题是如何在进程间交换数据。普通的MPI接口并不能直接用于GPU，因为数据都在显存上保存着，无法通过memcpy这样的经典函数进行复制，因此必须设置一个内存中的缓冲区，用于交换数据，但这样做的问题在于数据传输时总是需要在内存中转存一遍，非常影响性能。

为了解决这个问题，产生了CUDA-aware MPI技术。CUDA-aware MPI可以将MPI的接口直接用于显存，尽量避免了在内存中的转存，并且针对可能的硬件作了优化，如infiniband、nvlink等。本文中，我们选用openmpi，编译时需要加上参数--with-cuda=/path/to/cuda，如果更换了MPI(如将bashrc中的WM_MPLIB更改成SYSTEMOPENMPI)导致OpenFOAM报错可以重新编译Pstream，切换到目录OpenFOAM/src/Pstream，加载bashrc，然后执行wmake即可。

接下来我们将采用CUDA-aware MPI为OpenFOAM编写一个多GPU的预处理共轭梯度法。

§选择GPU

和单/多线程模型类似，我们还是需要使用cudaSetDevice来为后续的计算代码设置一个合适的GPU。前面已经说过，系统中GPU的编号为0,1,...,n-1，其中n是GPU的数量。但这仅限于单机系统，单/多线程模型也只能用于单机系统，不可能有两个GPU编号相同，而对于多机系统就不一样了，每个节点机上的GPU编号都是0到n-1，如节点1上有4个GPU，那么节点1上的GPU编号就是0-3，节点2上有6个GPU，则节点2上的GPU编号为0-5，那么如何为进程选择合适的GPU呢？

首先，各节点机上的进程数量不是随机安排的，我们应当看节点的硬件如何，此处当然就是看节点上的GPU数量，如节点1和节点2分别有4和6个GPU，那么我们应当分别给节点1和节点4安排4和6个进程，算法的思路非常简单，我们将MPI进程号模GPU数量传给cudaSetDevice即可。但要注意一个问题，如果进程号在节点上不是连续分布的，那么分配的结果不太好了，如节点1的4个进程编号为0、1、4、8，则节点1分配的GPU为0、1、0、0，负载严重不均。我们应当先计算进程在节点机内的编号，然后再计算应当使用哪个GPU。

这里，我们将进程按照是否可以共享内存进行分组，一般而言，组内进程就会位于同一个节点机上，然后我们使用组内的进程号选择GPU就可以避免前面说的负载不均的问题了。代码如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46

//========================= GPUSelector.h =========================
class GPUSelector {
public:
    static void Select();

private:
    int _gpu_no;   // 本进程使用哪个GPU

    GPUSelector();

    static GPUSelector *_gpu_selector;
};

//========================= GPUSelector.cpp =========================
GPUSelector *GPUSelector::_gpu_selector = nullptr;

GPUSelector::GPUSelector() {
    // 按照进程间是否可以共享内存来拆分通信器，如果可以共享内存则是同一个节点机
    MPI_Comm local_comm;
    MPI_Comm_split_type(
            MPI_COMM_WORLD,
            MPI_COMM_TYPE_SHARED,
            0, MPI_INFO_NULL,
            &local_comm
    );


    int local_rank;
    MPI_Comm_rank(local_comm, &local_rank);

    MPI_Comm_free(&local_comm);  // 只需要_local_rank，不需要通信器


    // 获取本机GPU数量
    int num_devices = 0;
    CHECK(cudaGetDeviceCount(&num_devices));

    _gpu_no = local_rank % num_devices;
}

void GPUSelector::Select() {
    if (_gpu_selector == nullptr)
        _gpu_selector = new GPUSelector{};

    CHECK(cudaSetDevice(_gpu_selector->_gpu_no));
}

§OpenFOAM Processor Boundary

在OpenFOAM中，有许多各类的边界，我们在此关注一类称为Processor的边界。Processor边界不是由人设置的，而是在执行decomposePar时由OpenFOAM生成的，它是那些原本是内部的边经过OpenFOAM拆分计算区域后变成外部的边，保存在processor0、processor1等文件夹下的constant/polyMesh中。例如windAroudBuildings的constant/polyMesh/boundary内容为：

// ...
5
(
    inlet
    {
        type            patch;
        nFaces          494;
        startFace       558513;
    }
    outlet
    {
        type            patch;
        nFaces          200;
        startFace       559007;
    }
    ground
    {
        type            wall;
        inGroups        1(wall);
        nFaces          6310;
        startFace       559207;
    }
    frontAndBack
    {
        type            symmetry;
        inGroups        1(symmetry);
        nFaces          1000;
        startFace       565517;
    }
    buildings
    {
        type            wall;
        inGroups        1(wall);
        nFaces          22721;
        startFace       566517;
    }
)

而processor1/constant/polyMesh/boundary的内容为：

// ...
9
(
    inlet
    {
        type            patch;
        nFaces          311;
        startFace       91498;
    }
    outlet
    {
        type            patch;
        nFaces          0;
        startFace       91809;
    }
    ground
    {
        type            wall;
        inGroups        1(wall);
        nFaces          1060;
        startFace       91809;
    }
    frontAndBack
    {
        type            symmetry;
        inGroups        1(symmetry);
        nFaces          215;
        startFace       92869;
    }
    buildings
    {
        type            wall;
        inGroups        1(wall);
        nFaces          2311;
        startFace       93084;
    }
    procBoundary1to0
    {
        type            processor;
        inGroups        1(processor);
        nFaces          702;
        startFace       95395;
        matchTolerance  0.0001;
        transform       unknown;
        myProcNo        1;
        neighbProcNo    0;
    }
    procBoundary1to2
    {
        type            processor;
        inGroups        1(processor);
        nFaces          1469;
        startFace       96097;
        matchTolerance  0.0001;
        transform       unknown;
        myProcNo        1;
        neighbProcNo    2;
    }
    procBoundary1to3
    {
        type            processor;
        inGroups        1(processor);
        nFaces          2;
        startFace       97566;
        matchTolerance  0.0001;
        transform       unknown;
        myProcNo        1;
        neighbProcNo    3;
    }
    procBoundary1to4
    {
        type            processor;
        inGroups        1(processor);
        nFaces          400;
        startFace       97568;
        matchTolerance  0.0001;
        transform       unknown;
        myProcNo        1;
        neighbProcNo    4;
    }
)

Processor Boundary对应的数据结构为Foam::processorFvPatchField，其详细信息请参考https://www.openfoam.com/documentation/guides/latest/doc/guide-bcs-constraint-processor.html，此处不再赘述。我们关注一下怎样保存和处理Processor边界信息。

对于一个线性方程而言，如果将其未知数拆分为多个，然后按未知数的个数将矩阵拆分成多个块，则非对角块就对应Processor Boundary，如下图所示：

couple_interface

非对角块的特点是非常稀疏且某些行没有非零元素，如windAroundBuildings中，进程1对应的几个块的非零元素个数仅为702、1469、2、400、0，而进程1对应的未知数有37086个，对角块的非零元素个数为97969。这就启示我们不必将它们当成普通的矩阵存储和计算，且在进行矩阵与向量乘法时，不必将它们对应的向量元素全部传输过来。

这里设计了一个类型Interface用于保存这些边界，Interface不仅保存边界，同时还负责进程间的数据传输，以及进行矩阵向量乘法时的部分运算。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

namespace MultiGPU {
class Interfaces {
public:
    // 本进程ID
    int my_proc_no;

    // 有多少个进程与此进程的区域有共同边界
    int n_par_interfaces;

    // 与此进程有共同边界的进程id
    std::vector<int> neighbour_proc_no;

    // 边界上各个面对应的单元行号，在对面进程中的行号是本进程的列号
    std::vector<SingleGPU::VectorInt> row_indices;

    // 边界上各个面的非零元素值
    std::vector<SingleGPU::Vector> values;

    // 向各个进程发送数据时的发送缓冲区，根据行号进行压缩
    std::vector<SingleGPU::Vector> send_buffers;

    // 从各个进程接收数据的接收缓冲区
    std::vector<SingleGPU::Vector> recv_buffers;

    // 使用非阻塞的方式发送和接收数据，记录发送状态和接收状态
    std::vector<MPI_Request> send_requests;

    std::vector<MPI_Request> recv_requests;

    // 是否是第一次发送，如果是则发送缓冲区空闲，不需要待前面的发送操作完成
    bool is_first_send = true;

    // 根据共同边界的个数设置各种缓冲区的数量
    void resize(int num_interfaces);

    // 压缩并发送数据
    void initMatrixInterfaces(
            const SingleGPU::Vector &v
    );

    // 作用到矩阵上，矩阵与向量乘法的add传入false，残量传入true
    void updateMatrixInterfaces(
            SingleGPU::Vector &Av,
            bool add = false
    );
};
} // namespace MultiGPU

Interfaces有个前置条件：边界上的面在其两侧进程中的编号顺序完全相同，即左侧进程中的row_indices[i]是非零元素的行号，则右侧进程中的row_indices[i]是对应元素的列号，反之亦然。因此进程可以根据自已的row_indices[i]来决定对面的进程在做矩阵与向量乘法时需要向量的哪些元素，我们按照顺序将这些元素放入send_buffers[i]中，然后使用MPI的接口将它发送出去。在此处，我们采用非阻塞的MPI接口，以尽量保证传输过程被计算掩盖。

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181

namespace MultiGPU {
__global__ void CompressVectorKernel(
        ValueType *__restrict__ compressed,
        const ValueType *__restrict__ vec,
        const IndexType *__restrict__ row_indices,
        const IndexType num_entries
) {
    IndexType i = threadIdx.x + blockIdx.x * blockDim.x;
    const int stride = blockDim.x * gridDim.x;

    for (; i < num_entries; i += stride)
        compressed[i] = vec[row_indices[i]];
}

void Interfaces::initMatrixInterfaces(const SingleGPU::Vector &v) {
    if (is_first_send) {
        for (int i = 0; i < n_par_interfaces; ++i) {
            const int num_entries = row_indices[i].size;
            send_buffers[i].resize(num_entries);
            recv_buffers[i].resize(num_entries);

            // 压缩数据
            LaunchKernelWithNumThreads(
                    CompressVectorKernel, num_entries, 0, 0,
                    send_buffers[i].data(),
                    v.data(), row_indices[i].data(), num_entries
            );
        }

        is_first_send = false;
    } else {
        int count = 0;

        std::vector<int> finished(n_par_interfaces);
        while (count < n_par_interfaces) {
            int out_count = 0;
            MPI_Waitsome(
                    (int)send_requests.size(),
                    send_requests.data(),
                    &out_count,
                    finished.data(),
                    MPI_STATUSES_IGNORE
            );

            count += out_count;

            for (int i = 0; i < out_count; ++i) {
                int proc = finished[i];

                const int num_entries = row_indices[proc].size;
                send_buffers[proc].resize(num_entries);

                // 压缩数据
                LaunchKernelWithNumThreads(
                        CompressVectorKernel, num_entries, 0, 0,
                        send_buffers[proc].data(),
                        v.data(), row_indices[proc].data(), num_entries
                );
            }
        }
    }

    // 使用cuda with mpi，则无需将数据传输到Host，可以直接在GPU之间交换数据，
    // 由MPI决定传输细节

    // 使用非阻塞接口通信必须保证数据已经计算到缓冲内了，
    // MPI_Send会调用cuDeviceSynchronize进行同步，所以非阻塞接口不需要
    CHECK(cudaStreamSynchronize(0));

    for (int i = 0; i < n_par_interfaces; ++i) {
        const int num_entries = row_indices[i].size;

        recv_buffers[i].resize(num_entries);

        if (my_proc_no > neighbour_proc_no[i]) {
            // 将压缩后的数据发送给对方
            MPI_Isend(
                    send_buffers[i].data(), num_entries, MPI_DOUBLE,
                    neighbour_proc_no[i], 0, MPI_COMM_WORLD, &send_requests[i]
            );

            // 从对方接收数据
            MPI_Irecv(
                    recv_buffers[i].data(), num_entries, MPI_DOUBLE,
                    neighbour_proc_no[i], 0, MPI_COMM_WORLD, &recv_requests[i]
            );
        } else {
            // 从对方接收数据
            MPI_Irecv(
                    recv_buffers[i].data(), num_entries, MPI_DOUBLE,
                    neighbour_proc_no[i], 0, MPI_COMM_WORLD, &recv_requests[i]
            );

            // 将压缩后的数据发送给对方
            MPI_Isend(
                    send_buffers[i].data(), num_entries, MPI_DOUBLE,
                    neighbour_proc_no[i], 0, MPI_COMM_WORLD, &send_requests[i]
            );
        }
    }
}

__global__ void UpdateInterfacesAddKernel(
        ValueType *__restrict__ Av,
        const ValueType *__restrict__ values,
        const IndexType *__restrict__ row_indices,
        const ValueType *__restrict__ vec,
        const IndexType num_entries
) {
    IndexType i = threadIdx.x + blockIdx.x * blockDim.x;
    const int stride = blockDim.x * gridDim.x;

    for (; i < num_entries; i += stride) {
        ValueType t = values[i] * vec[i];

        atomicAdd(&Av[row_indices[i]], t);
    }
}

__global__ void UpdateInterfacesMinusKernel(
        ValueType *__restrict__ Av,
        const ValueType *__restrict__ values,
        const IndexType *__restrict__ row_indices,
        const ValueType *__restrict__ vec,
        const IndexType num_entries
) {
    IndexType i = threadIdx.x + blockIdx.x * blockDim.x;
    const int stride = blockDim.x * gridDim.x;

    for (; i < num_entries; i += stride) {
        ValueType t = -values[i] * vec[i];

        atomicAdd(&Av[row_indices[i]], t);
    }
}

void Interfaces::updateMatrixInterfaces(
        SingleGPU::Vector &Av,
        bool add
) {
    int count = 0;

    std::vector<int> finished(n_par_interfaces);
    while (count < n_par_interfaces) {
        int out_count = 0;

        MPI_Waitsome(
                (int)recv_requests.size(),
                recv_requests.data(),
                &out_count,
                finished.data(),
                MPI_STATUSES_IGNORE
        );

        count += out_count;

        for (int i = 0; i < out_count; ++i) {
            int proc = finished[i];

            const IndexType num_entries = recv_buffers[proc].size;

            // 矩阵与向量乘法，用减法
            if (!add) {
                LaunchKernelWithNumThreads(
                        UpdateInterfacesMinusKernel, num_entries, 0, 0,
                        Av.data(),
                        values[proc].data(), row_indices[proc].data(),
                        recv_buffers[proc].data(), num_entries
                );
            } else {  // 残量用加法
                LaunchKernelWithNumThreads(
                        UpdateInterfacesAddKernel, num_entries, 0, 0,
                        Av.data(),
                        values[proc].data(), row_indices[proc].data(),
                        recv_buffers[proc].data(), num_entries
                );
            }
        }
    }
}
} // namespace MultiGPU

§矩阵和向量基本操作

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58

namespace MultiGPU {
// 向量加法不需要进程间通信，和单GPU完全一样

// return v1^T * v2
ValueType Dot(
        const SingleGPU::Vector &v1,
        const SingleGPU::Vector &v2
) {
    ValueType local_dot = SingleGPU::Dot(v1, v2);

    ValueType global_dot = 0;

    // 支持CUDA-aware MPI的MPI接口也可以用于内存
    MPI_Allreduce(
            &local_dot, &global_dot, 1, MPI_DOUBLE,
            MPI_SUM, MPI_COMM_WORLD
    );

    return global_dot;
}

ValueType Norm2(const SingleGPU::Vector &v) {
    return sqrt(MultiGPU::Dot(v, v));
}

// result = A * v
void SpMv(
        SingleGPU::Vector *result,
        const SingleGPU::CsrMatrix &A,
        const SingleGPU::Vector &v,
        Interfaces &interfaces
) {
    if (A.num_cols != v.size)
        throw std::runtime_error{"MultiGPU::SpMv: A.num_cols != v.size"};

    interfaces.initMatrixInterfaces(v);
    SingleGPU::SpMv(result, A, v);
    interfaces.updateMatrixInterfaces(*result);
}

// r = b - Ax
void Residual(
        SingleGPU::Vector *r,
        const SingleGPU::CsrMatrix &A,
        const SingleGPU::Vector &b,
        const SingleGPU::Vector &x,
        Interfaces &interfaces
) {
    if (A.num_cols != x.size)
        throw std::runtime_error{"MultiGPU::Residual: A.num_cols != x.size!"};
    if (A.num_rows != b.size)
        throw std::runtime_error{"MultiGPU::Residual: A.num_rows != b.size!"};

    interfaces.initMatrixInterfaces(x);
    SingleGPU::Residual(r, A, b, x);
    interfaces.updateMatrixInterfaces(*r, true);  // 注意这里传true
}
} // namespace MultiGPU

§预处理共轭梯度法

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

namespace MultiGPU {
void PCG(
        SingleGPU::Vector *x,
        const SingleGPU::CsrMatrix &A,
        MultiGPU::Interfaces &interfaces,
        const SingleGPU::Vector &b,
        const SingleGPU::Preconditioner &P,
        ValueType error = 1e-5,
        IndexType max_steps = 0,
) {
    if (A.num_cols != b.size)
        throw std::runtime_error{"PCGFoam: A.num_cols != b.size!"};

    const IndexType N = A.num_rows;

    // r_0 = b - Ax_0
    SingleGPU::Vector r(N);
    MultiGPU::Residual(&r, A, b, *x, interfaces);

    ValueType norm_r = MultiGPU::Norm2(r);

    // 初始解精度已经足够高
    if (norm_r < error)
        return;

    // z_0 = P^-1 * r_0
    SingleGPU::Vector z(N);
    P.precondition(&z, r);

    // p_1 = z_0
    SingleGPU::Vector p = z;

    // r_k^T * z_k and r_{k-1}^T * z_{k-1}
    ValueType rho0 = MultiGPU::Dot(r, z);
    ValueType rho1 = 0;

    // q_k
    SingleGPU::Vector q(N);
    
    if (max_steps <= 0)
        max_steps = std::numeric_limits<IndexType>::max();

    for (IndexType k = 1; k <= max_steps; ++k) {
        // q_k = A * p_k
        MultiGPU::SpMv(&q, A, p, interfaces);

        // xi_k = (r_{k-1}^T * z_{k-1}) / (p_k^T * q_k)
        ValueType xi = rho0 / MultiGPU::Dot(p, q);

        // x_k = x_{k-1} + xi_k * p_k
        SingleGPU::Add(x, *x, p, xi);

        // r_k = r{k-1} - xi_k * q_k
        SingleGPU::Add(&r, r, q, -xi);

        norm_r = MultiGPU::Norm2(r);
        
        if (norm_r < error)
            break;

        // z_k = P^-1 * r_k
        P.precondition(&z, r);

        rho1 = rho0;
        rho0 = MultiGPU::Dot(r, z);

        // p_{k+1} = z_k + u_k * p_k
        SingleGPU::Add(&p, z, p, rho0 / rho1);
    }
}
} // namespace MultiGPU

值得注意的是，在OpenFOAM的多进程实现中，预处理子只用了对角块进行构造，没有使用边界信息，这是因为边界信息本就影响甚微而且考虑边界信息构造预处理子非常麻烦，通信频繁，并不划算，这里也遵循OpenFOAM的实现。

§集成

集成到OpenFOAM可以通过继承类型Foam::lduMatrix::solver来实现：

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

//========================= cuFoamPCG.H ==========================//
#ifndef CUFOAM_PCG_H
#define CUFOAM_PCG_H

#include "lduMatrix.H"

namespace Foam {

class cuFoamPCG : public lduMatrix::solver {
    cuFoamPCG(const cuFoamPCG &) = delete;

    cuFoamPCG &operator=(const cuFoamPCG &) = delete;

public:
    TypeName("cuFoamPCG");

    cuFoamPCG(
            const word &fieldName,
            const lduMatrix &matrix,
            const FieldField<Field, scalar> &interfaceBouCoeffs,
            const FieldField<Field, scalar> &interfaceIntCoeffs,
            const lduInterfaceFieldPtrsList &interfaces,
            const dictionary &solverControls
    );

    virtual ~cuFoamPCG() {}

    virtual solverPerformance solve(
            scalarField &psi,
            const scalarField &source,
            const direction cmpt = 0
    ) const;
};

}  // namespace Foam

#endif // ~CUFOAM_CG_H


//========================= cuFoamPCG.C ==========================//
namespace Foam {

defineTypeNameAndDebug(cuFoamPCG, 0);

// 只能用于对称的(symmetric)方程
lduMatrix::solver::addsymMatrixConstructorToTable<cuFoamPCG>
        addCGSolverSymMatrixConstructorToTable_;

} // ~namespace Foam

Foam::cuFoamPCG::cuFoamPCG(
        const word &fieldName,
        const lduMatrix &matrix,
        const FieldField<Field, scalar> &interfaceBouCoeffs,
        const FieldField<Field, scalar> &interfaceIntCoeffs,
        const lduInterfaceFieldPtrsList &interfaces,
        const dictionary &solverControls
) :
        lduMatrix::solver(
                fieldName,
                matrix,
                interfaceBouCoeffs,
                interfaceIntCoeffs,
                interfaces,
                solverControls
        ) {}


Foam::solverPerformance Foam::cuFoamPCG::solve(
        scalarField &psi,
        const scalarField &source,
        const direction cmpt
) const {
    IndexType N = matrix().diag().size();

    // 将ldu格式的矩阵转换成csr格式
    Host::CsrMatrix A = getSparseCSRMat(matrix());

    // 设置使用的GPU，每个进程对应一个GPU
    if (Pstream::parRun())
        GPUSelector::Select();


    // 将数据传输到GPU上
    SingleGPU::CsrMatrix d_A{A};

    SingleGPU::Vector d_b(N);
    SingleGPU::Vector d_x(N);

    uploadVector(d_b, source);
    uploadVector(d_x, psi);

    // 预处理子
    word precond_name = controlDict_.lookupOrDefault<word>("preconditioner", "none");

    SingleGPU::Preconditioner *P = nullptr;

    if (precond_name == "DILU")
        P = new SingleGPU::DILUPreconditioner(d_A);
    else if (precond_name == "DIC")
        P = new SingleGPU::DICPreconditioner(d_A);
    else if (precond_name == "Diagonal")
        P = new SingleGPU::DiagonalPreconditioner(d_A);
    else if (precond_name == "none")
        P = new SingleGPU::NoPreconditioner(d_A);
    else {
        FatalErrorInFunction
                << "Preconditioner must be one of "
                   "DILU, DIC, Diagonal, and none!"
                << nl << exit(FatalError);
    }

    IndexType max_steps = maxIter_;  // 最大迭代步数
    IndexType min_steps = minIter_;  // 最小迭代步数
    ValueType error = tolerance_;    // 绝对误差
    ValueType rel_error = relTol_;   // 相对误差

    std::vector<ValueType> res(max_steps);

    if (Pstream::parRun()) {
        MultiGPU::Interfaces gpu_interface = getInterfaces(
                matrix(), interfaceBouCoeffs(), interfaces()
        );

        // 这个函数请参考源代码文件，上节仅给出了核心功能，
        // 未实现OpenFOAM需要的许多控制功能
        MultiGPU::PCGFoam(
                &d_x, d_A, gpu_interface, d_b, *P,
                error, rel_error, max_steps, min_steps, &res
        );
    } else {
        SingleGPU::PCGFoam(
                &d_x, d_A, d_b, *P,
                error, rel_error, max_steps, min_steps, &res
        );
    }

    solverPerformance solverPerf("cuFoam" + precond_name + "PCG", fieldName());
    solverPerf.initialResidual() = res.front();
    solverPerf.finalResidual() = res.back();
    solverPerf.nIterations() = static_cast<int>(res.size()) - 1;

    delete P;

    return solverPerf;
}

这段代码向OpenFOAM注册了一个名为cuFoamPCG的求解器，在fvSolution的solver字段中指定为此名字就可以调用这个求解器。需要实现的接口主要为Foam::lduMatrix::solve，这里给出了大致流程，其中的getSparseCsrMat和getInterfaces没有给出，它们位于openfoam/common下，读者请自行查看。需要注意OpenFOAM中的Processor边界对应的类型为Foam::processorFvPatchField，其继承关系非常复杂，阅读OpenFOAM源码时需要仔细梳理。

相关文章：

加载评论