Intel TBB 开发指南 3 parallel_reduce

原文

循环可以进行归约，如以下求和：

float SerialSumFoo( float a[], size_t n ) {
    float sum = 0;
    for( size_t i=0; i!=n; ++i )
        sum += Foo(a[i]);
    return sum;
}

如果迭代是独立的，你可以使用模板类 parallel_reduce 并行化此循环，如下所示：

float ParallelSumFoo( const float a[], size_t n ) {
    SumFoo sf(a);
    parallel_reduce( blocked_range<size_t>(0,n), sf );
    return sf.my_sum;
}

SumFoo 类指定了归约的细节，例如如何累加和组合它们。 这是类 SumFoo 的定义：

class SumFoo {
    float* my_a;
public:
    float my_sum;
    void operator()( const blocked_range<size_t>& r ) {
        float *a = my_a;
        float sum = my_sum;
        size_t end = r.end();
        for( size_t i=r.begin(); i!=end; ++i )
            sum += Foo(a[i]);
        my_sum = sum;
    }


    SumFoo( SumFoo& x, split ) : my_a(x.my_a), my_sum(0) {}


    void join( const SumFoo& y ) {my_sum+=y.my_sum;}


    SumFoo(float a[] ) :
        my_a(a), my_sum(0)
    {}
};

请注意与来自 parallel_for 的类 ApplyFoo 的区别。 首先，operator() 不是 const。 这是因为它必须更新 SumFoo::my_sum。 其次，SumFoo 有一个拆分构造函数和一个方法 join，必须存在才能使 parallel_reduce 工作。 拆分构造函数将原始对象的引用和类型为 split 的伪参数作为参数，该参数由库定义。 虚拟参数将拆分构造函数与复制构造函数区分开来。

在示例中，operator() 的定义将局部临时变量（a、sum、end）用于循环内访问的标量值。 这种技术可以通过让编译器清楚地知道这些值可以保存在寄存器而不是内存中来提高性能。 如果这些值太大而无法放入寄存器中，或者它们的地址以编译器无法跟踪的方式获取，则该技术可能无济于事。 对于典型的优化编译器，仅对写入的变量（例如示例中的 sum）使用局部临时变量就足够了，因为这样编译器就可以推断出循环不会写入任何其他位置，并将其他读取提升到外部 循环。

当工作线程可用时，由任务调度程序决定，parallel_reduce 调用拆分构造函数为工作线程创建子任务。 当子任务完成时，parallel_reduce 使用方法 join 来累积子任务的结果。 下图顶部的图形显示了当工作程序可用时发生的拆分连接序列：

上图中的箭头表示时间顺序。 拆分构造函数可能会在对象 x 用于缩减的前半部分时并发运行。 因此，创建 y 的拆分构造函数的所有操作都必须相对于 x 成为线程安全的。 因此，如果拆分构造函数需要递增与其他对象共享的引用计数，则应使用原子递增。
如果 worker 不可用，则使用减少前半部分的相同主体对象减少迭代的后半部分。 那就是下半场的减持开始，上半场的减持结束。

由于如果 worker 不可用，则不使用 split/join，parallel_reduce 不一定会进行递归拆分。
由于同一个主体可能用于累积多个子范围，因此 operator() 不丢弃较早的累积是至关重要的。 下面的代码显示了 SumFoo::operator() 的错误定义。

class SumFoo {
    ...
public:
    float my_sum;
    void operator()( const blocked_range<size_t>& r ) {
        ...
        float sum = 0;  // WRONG – should be 'sum = my_sum".
        ...
        for( ... )
            sum += Foo(a[i]);
        my_sum = sum;
    }
    ...
};

由于错误，主体返回最后一个子范围的部分和，而不是 parallel_reduce 应用它的所有子范围。
parallel_reduce 的分区器和粒度规则与 parallel_for 相同。
parallel_reduce 推广到任何关联操作。 通常，拆分构造函数会做两件事：

• 复制运行循环体所需的只读信息。
• 将归约变量初始化为操作的标识元素。
  join 方法应该进行相应的合并。 你可以同时进行多次缩减：可以使用单个 parallel_reduce 收集最小值和最大值。

归约运算可以是非交换的。 如果浮点加法被字符串连接替换，该示例仍然有效。

Advanced Example

更高级的关联操作的一个例子是找到 Foo(i) 最小化的索引。 串行版本可能如下所示：

long SerialMinIndexFoo( const float a[], size_t n ) {
    float value_of_min = FLT_MAX;        // FLT_MAX from <climits>
    long index_of_min = -1;
    for( size_t i=0; i<n; ++i ) {
        float value = Foo(a[i]);
        if( value<value_of_min ) {
            value_of_min = value;
            index_of_min = i;
        }
    }
    return index_of_min;
}

该循环记录迄今为止找到的最小值以及该值的索引。 这是循环迭代之间携带的唯一信息。 要将循环转换为使用 parallel_reduce，函数对象必须跟踪携带的信息，以及当迭代分布在多个线程上时如何合并这些信息。 此外，函数对象必须记录指向 a 的指针以提供上下文。
以下代码显示了完整的函数对象。

class MinIndexFoo {
    const float *const my_a;
public:
    float value_of_min;
    long index_of_min;
    void operator()( const blocked_range<size_t>& r ) {
        const float *a = my_a;
        for( size_t i=r.begin(); i!=r.end(); ++i ) {
           float value = Foo(a[i]);
           if( value<value_of_min ) {
               value_of_min = value;
               index_of_min = i;
           }
        }
    }


    MinIndexFoo( MinIndexFoo& x, split ) :
        my_a(x.my_a),
        value_of_min(FLT_MAX),    // FLT_MAX from <climits>
        index_of_min(-1)
   {}


    void join( const SumFoo& y ) {
        if( y.value_of_min<value_of_min ) {
            value_of_min = y.value_of_min;
            index_of_min = y.index_of_min;
        }
    }


    MinIndexFoo( const float a[] ) :
        my_a(a),
        value_of_min(FLT_MAX),    // FLT_MAX from <climits>
        index_of_min(-1),
    {}
};

现在 SerialMinIndex 可以使用 parallel_reduce 重写，如下所示：

long ParallelMinIndexFoo( float a[], size_t n ) {
    MinIndexFoo mif(a);
    parallel_reduce(blocked_range<size_t>(0,n), mif );
    return mif.index_of_min;
}

Advanced Topic: Other Kinds of Iteration Spaces

到目前为止的示例都使用了类 blocks_range<T> 来指定范围。 这个类在许多情况下都很有用，但它并不适合所有情况。 你可以使用 oneTBB 来定义自己的迭代空间对象。 该对象必须通过提供一个基本的拆分构造函数、一个可选的比例拆分构造函数（伴随着启用其使用的特征值）和两个谓词方法来指定如何将其拆分为子空间。 如果你的类称为 R，则方法和构造函数应如下所示：

class R {
    // True if range is empty
    bool empty() const;
    // True if range can be split into non-empty subranges
    bool is_divisible() const;
    // Splits r into subranges r and *this
    R( R& r, split );
    // Splits r into subranges r and *this in proportion p
    R( R& r, proportional_split p );
    // Allows usage of proportional splitting constructor
    static const bool is_splittable_in_proportion = true;
    ...
};

如果范围为空，方法 empty 应该返回 true。如果范围可以拆分为两个非空子空间，则方法 is_divisible 应该返回 true，这样的拆分值得开销。基本的拆分构造函数应该有两个参数：

• 第一个 R 类型
• 第二个 oneapi::tbb::split 类型

第二个参数没有被实际使用；它仅用于将构造函数与普通的复制构造函数区分开来。基本的拆分构造函数应该尝试将 r 大致拆分为两半，并将 r 更新为前半部分，并将构造的对象设置为后半部分。
与基本拆分构造函数不同，比例拆分构造函数是可选的，并采用 oneapi::tbb::proportional_split 类型的第二个参数。该类型具有返回比例值的 left 和 right 方法。应该使用这些值来相应地拆分 r，使更新后的 r 对应于比例的左侧部分，而构造的对象对应于右侧部分。只有在类中定义了静态常量 is_splittable_in_proportion 并赋值为 true 时，才会使用比例拆分构造函数。
两个拆分构造函数都应该保证更新后的 r 部分和构造的对象不为空。只有当 r.is_divisible 为真时，并行算法模板才会调用 r 上的拆分构造函数。
迭代空间不必是线性的。查看 oneapi/tbb/blocked_range2d.h 以获得二维范围的示例。它的拆分构造函数尝试沿其最长轴拆分范围。当与parallel_for 一起使用时，它会导致循环以提高缓存使用率的方式“递归阻塞”。这种良好的缓存行为意味着，在blocked_range2d<T> 上使用parallel_for 可以使循环比顺序等效的循环运行得更快，即使在单个处理器上也是如此。