OpenMP/Tasks

如果你跟着上一章并做了练习，你就会发现我们开发的求和函数很不稳定。（公平地说，我们设置了数据来展示这种不稳定性。）实际上，简单地求和的 $O(n)$ 舍入误差：误差与我们要累加的元素数量成正比。我们可以通过更改求和算法来解决这个问题。有三个候选的稳定求和算法

将数字从小到大排序，然后像以前一样累加它们。问题：这将问题的复杂度从线性变为 $O(n\log n)$ .
使用Kahan 算法。虽然该算法需要线性时间并且非常稳定，但在实践中它相当缓慢，并且比我们的第三种选择更难并行化，第三种选择是
分治递归。这种算法的舍入误差是 $O(\log n)$ ，即与元素数量的对数成正比。

基本的分治求和在 C 中很容易表达

float sum(const float *a, size_t n)
{
    // base cases
    if (n == 0) {
        return 0;
    }
    else if (n == 1) {
        return *a;
    }

    // recursive case
    size_t half = n / 2;
    return sum(a, half) + sum(a + half, n - half);
}

如果你在上一章为该程序开发的程序中使用这个 sum 的定义，你会发现它产生了完全预期的结果。但是这个算法没有循环，那么我们如何使用 OpenMP 制作并行版本呢？

我们将使用 OpenMP 中的任务构造，将问题视为任务并行而不是数据并行。在第一个版本中，sum 的任务递归版本如下所示

float sum(const float *a, size_t n)
{
    // base cases
    if (n == 0) {
        return 0;
    }
    else if (n == 1) {
        return 1;
    }

    // recursive case
    size_t half = n / 2;
    float x, y;

    #pragma omp parallel
    #pragma omp single nowait
    {
        #pragma omp task shared(x)
        x = sum(a, half);
        #pragma omp task shared(y)
        y = sum(a + half, n - half);
        #pragma omp taskwait
        x += y;
    }
    return x;
}

我们引入了两个任务，每个任务都设置一个被声明为 shared 的变量，该变量与另一个任务共享。如果我们没有声明共享变量，每个任务都会设置自己的局部变量，然后丢弃结果。然后我们使用 #pragma omp taskwait 等待任务完成，并结合递归结果。

你可能会对紧跟着 #pragma omp single nowait 的 #pragma omp parallel 感到惊讶。问题是第一个pragma会导致池中的所有线程执行下一段代码。single 指令会导致除了一个线程（通常是第一个遇到该代码块的线程）以外的所有线程不执行它，而 nowait 会关闭 single 上的屏障；在封闭的 parallel 区域上已经有一个屏障，其他线程会蜂拥而至。

不幸的是，如果你真的尝试运行这段代码，你会发现它仍然不是特别快。原因是任务太细粒度了：在递归树的底部附近， ${\tfrac {n}{2}}$ 次调用将两个元素的数组分成每个处理一个元素的子任务。我们可以通过引入除了基本情况和递归情况之外的“中间情况”来解决这个问题，该情况是递归的，但不涉及设置并行任务：如果递归达到预先指定的截止值，它将不再尝试为 OpenMP 线程池设置任务，而是只进行递归求和。

练习：在递归中引入额外的情况，并衡量程序的速度。不要提前查看下一个程序，因为它包含了这个练习的答案。

现在我们实际上有两个递归合二为一：一个带并行递归情况，另一个是串行递归。我们可以通过在每个级别进行更少的检查来解开这两个递归，以获得更好的性能。我们还将并行设置代码分离到驱动函数中。

#include <stddef.h>

#define CUTOFF 100  // arbitrary

static float parallel_sum(const float *, size_t);
static float serial_sum(const float *, size_t);

float sum(const float *a, size_t n)
{
    float r;

    #pragma omp parallel
    #pragma omp single nowait
    r = parallel_sum(a, n);
    return r;
}

static float parallel_sum(const float *a, size_t n)
{
    // base case
    if (n <= CUTOFF) {
        return serial_sum(a, n);
    }

    // recursive case
    float x, y;
    size_t half = n / 2;

    #pragma omp task shared(x)
    x = parallel_sum(a, half);
    #pragma omp task shared(y)
    y = parallel_sum(a + half, n - half);
    #pragma omp taskwait
    x += y;

    return x;
}

static float serial_sum(const float *a, size_t n)
{
    // base cases
    if (n == 0) {
        return 0.;
    }
    else if (n == 1) {
        return a[0];
    }

    // recursive case
    size_t half = n / 2;
    return serial_sum(a, half) + serial_sum(a + half, n - half);
}

当并行任务内的代码花费更多时间进行计算，而花费更少时间进行内存访问时，这种技术效果更好，因为这些内存访问可能需要在处理器之间进行同步。