c++ - 红黑高斯赛德尔和 OpenMP

Question

与 MPICH 相比，我试图用 OpenMP 证明一个观点，我编写了以下示例来演示在 OpenMP 中实现一些高性能是多么容易。

Gauss-Seidel 迭代被分成两个独立的运行，这样在每次扫描中，每个操作都可以按任何顺序执行，并且每个任务之间应该没有依赖关系。所以理论上每个处理器都不必等待另一个进程来执行任何类型的同步。

我遇到的问题是，无论问题大小如何，我发现只有 2 个处理器的微弱加速，超过 2 个处理器甚至可能更慢。许多其他线性并行例程我可以获得很好的缩放比例，但是这个很棘手。

我担心的是我无法向编译器“解释”我对数组执行的操作是线程安全的，因此它无法真正有效。

请参阅下面的示例。

有人知道如何使用 OpenMP 使其更有效吗？

void redBlackSmooth(std::vector<double> const & b,
                    std::vector<double> & x,
                    double h)
{
    // Setup relevant constants.
    double const invh2 = 1.0/(h*h);
    double const h2 = (h*h);
    int const N = static_cast<int>(x.size());
    double sigma = 0;

    // Setup some boundary conditions.
    x[0] = 0.0;
    x[N-1] = 0.0;

    // Red sweep.
    #pragma omp parallel for shared(b, x) private(sigma)
    for (int i = 1; i < N-1; i+=2)
    {
        sigma = -invh2*(x[i-1] + x[i+1]);
        x[i] = (h2/2.0)*(b[i] - sigma);
    }

    // Black sweep.
    #pragma omp parallel for shared(b, x) private(sigma)
    for (int i = 2; i < N-1; i+=2)
    {
        sigma = -invh2*(x[i-1] + x[i+1]);
        x[i] = (h2/2.0)*(b[i] - sigma);
    }
}

添加:我现在也尝试过使用原始指针实现，它具有与使用 STL 容器相同的行为，因此可以排除它是来自 STL 的一些伪关键行为。

Answer 1

首先，确保 x vector 与缓存边界对齐。我做了一些测试，如果我强制对齐内存，我的机器(核心双核)上的代码得到了 100% 的改进:

double * x;
const size_t CACHE_LINE_SIZE = 256;
posix_memalign( reinterpret_cast<void**>(&x), CACHE_LINE_SIZE, sizeof(double) * N);

其次，您可以尝试为每个线程分配更多的计算(通过这种方式您可以保持缓存行分离)，但我怀疑 openmp 已经在幕后做了类似的事情，所以它可能对大 N 毫无值(value).

在我的例子中，当 x 不是缓存对齐时，这个实现要快得多。

const int workGroupSize = CACHE_LINE_SIZE / sizeof(double);
assert(N % workGroupSize == 0); //Need to tweak the code a bit to let it work with any N
const int workgroups = N / workGroupSize;
int j, base , k, i;

#pragma omp parallel for shared(b, x) private(sigma, j, base, k, i)
for ( j = 0; j < workgroups; j++ ) {
    base = j * workGroupSize;
    for (int k = 0; k < workGroupSize; k+=2)
    {
        i = base + k + (redSweep ? 1 : 0);
        if ( i == 0 || i+1 == N) continue;
        sigma = -invh2* ( x[i-1] + x[i+1] );
        x[i] = ( h2/2.0 ) * ( b[i] - sigma );
    }
}

总而言之，您肯定有缓存争用的问题，但考虑到 openmp 的工作方式(遗憾的是我不熟悉它)它应该足以处理正确分配的缓冲区。