新的运算符 + OpenMP 动态调度子句

Question

我一直在使用 C++ 代码来执行量子化学、原子和分子任务，其中隐含着大量使用数组（1D、2D、3D 等）的工作。我有一整堂课叫array 来处理这个问题。当然，从一开始，最基本的成员函数就是那些为这些数组动态分配内存、调整大小或删除它们的函数。

data = new double **[row]();
#pragma omp for schedule(static) nowait
for(unsigned int i = 0; i < row; ++i)
{
    data[i] = new double *[column]();
}

现在我正在做的是使用 OpenMP 加速这些例程。对于大部分例程，我一直在使用schedule(static) nowait 子句将我的循环分成step/threads 的块，因为这些块花费几乎相同的时间由它们的线程处理。

但是对于像上面这样的循环，通过多次调用new 运算符，我有一种（不好的）感觉，这些循环的块在它们的线程中执行的时间并不相同，从某种意义上说，我应该考虑改用schedule(dynamic, chunk_size)。

你们同意吗？动态分配不是那么简单的任务，而且可能很昂贵，因此动态分配块的执行时间可能会有所不同。

实际上，我不确定我是否在堆栈碎片或类似的事情上犯了任何错误。欢迎提出建议。

PS.：我正在使用nowait 子句来尽量减少隐式障碍的瓶颈。

Answer 1

如果您使用默认的new 运算符，您的特定循环可能不会为并行提供太多机会，因为堆是单个资源并且需要通过互斥体对其进行序列化。但是，假设您有其他希望使用 OpenMP 的循环，以下内容应该会有所帮助。

来自 OpenMP 3.1 specification：

static 当 schedule(static, chunk_size) 被指定时，迭代被分成大小为 chunk_size，chunk按照线程号的顺序循环分配给团队中的线程。

当没有指定chunk_size时，迭代空间被分成大小大致相等的chunk，每个线程最多分配一个chunk。请注意，在这种情况下，块的大小是未指定的。

dynamic 当 schedule(dynamic, chunk_size) 指定时，迭代次数为 当线程请求它们时，以块的形式分发给团队中的线程。每个线程执行一个迭代块，然后请求另一个块，直到没有块需要分发。

每个块都包含 chunk_size 次迭代，但要分发的最后一个块除外，它的迭代次数可能更少。

当没有指定chunk_size时，默认为1。

在您的情况下，您没有指定 chunk_size，因此未指定每个任务的迭代次数。

一般来说，我更喜欢对线程数和每个任务正在执行的迭代次数进行一些控制。我发现（在 Windows 上，使用 mingw-w64 编译），开始新的工作块的任务有很大的开销，所以给它们尽可能大的块是有益的。我倾向于使用动态（尽管我可以将静态用于固定执行时间的任务），并将 chunk_size 设置为循环计数除以线程数。在您的情况下，如果您怀疑任务执行时间不均匀，您可以将其除以 2 或 4。

// At the top of a C++ file:
static int NUM_THREADS = omp_get_num_procs();

// Then for your loop construct (I'm using a combined parallel for here):
#pragma omp parallel for num_threads(NUM_THREADS) \
   schedule(dynamic, row / NUM_THREADS / 2)
for(unsigned int i = 0; i < row; ++i)
{
   data[i] = new double *[column]();
}

另请注意，如果您不设置 num_threads，则默认值为 nthreads-var，由 omp_get_max_threads 确定。

关于 nowait 子句，显然要确保您没有在 loop 构造之外使用 data。我在上面使用了组合并行循环结构，这意味着不能指定nowait。

Answer 2

有一种更简洁的方法可以完成所有这些操作：

std::vector<double> actual_data(omp_get_num_procs() * column);

std::vector<double *> data(omp_get_num_procs());
for (unsigned i = 0; i < row; i++) {
    data[i] = &(actual_data[i * column]);
}

您现在有一个一次性分配的数组，其中包含指向该数组的指针数组。在并行算法中，您可以使用data[i][j] 来获取您想要的成员，几乎为零开销（开销发生在编译时）。

唯一的潜在风险是false sharing，因为矩阵的行可能在它们的端点共享一个缓存线。

Memory management is automatic;不需要free 任何东西。

Answer 3

如果您打算进行大量基本的线性代数运算 (BLAS)，那么我建议您不要将数组数组用于多维数组。不幸的是，多维静态和动态数组的 C 语法破坏了对称性。这导致 Java 和 C# 使用相同的语法，因此人们经常认为，当他们分配 2D 数组时，他们得到的东西与静态数组相同，但来自堆而不是堆栈。你的二维数组实际上是Jagged arrays。

根据我在 BLAS 操作（以及图像处理）方面的经验，您需要为 2D 数组提供连续的内存块。因此，对于 nxn 矩阵，您应该像 double[n*n] 一样分配并以 data[i*n+j] 的形式访问它，或者创建一个 C++ 矩阵类，您可以在其中访问像 matrix.get(i,j) 这样的矩阵。

不要并行分配内存。无论如何，大多数 BLAS 操作都将受内存限制，因此 OpenMP 仅对 3 级 BLAS 操作有用，例如矩阵乘法、LU 分解、cholesky 分解，如O(n^3)。