在 openMP 中并行化

Question

我想使用 OpenMP 并行化以下代码

for(m=0; m<r_c; m++)
{
    for(n=0; n<c_c; n++)
    {
        double value = 0.0;
        for(j=0; j<r_b; j++)
            for(k=0; k<c_b; k++)
            {
                double a;
                if((m-j)<0 || (n-k)<0 || (m-j)>r_a || (n-k)>c_a)
                    a = 0.0;
                else
                    a = h_a[((m-j)*c_a) + (n-k)];
                //printf("%lf\t", a);
                value += h_b[(j*c_b) + k] * a;
            }
        h_c[m*c_c + n] = value;
        //printf("%lf\t", h_c[m*c_c + n]);
    }
    //cout<<"row "<<m<<" completed"<<endl;
}

在此我希望每个线程同时执行“for j”和“for k”。 我正在尝试在“for m”循环之前使用 pragma omp parallel for，但没有得到正确的结果。 我怎样才能以优化的方式做到这一点。提前致谢。

Answer 1

根据您要并行化的具体循环，您有三个选项：

#pragma omp parallel
{
#pragma omp for    // Option #1
    for(m=0; m<r_c; m++)
    {
        for(n=0; n<c_c; n++)
        {
            double value = 0.0;
#pragma omp for    // Option #2
            for(j=0; j<r_b; j++)
                for(k=0; k<c_b; k++)
                {
                    double a;
                    if((m-j)<0 || (n-k)<0 || (m-j)>r_a || (n-k)>c_a)
                        a = 0.0;
                    else
                        a = h_a[((m-j)*c_a) + (n-k)];
                    //printf("%lf\t", a);
                    value += h_b[(j*c_b) + k] * a;
                }
                h_c[m*c_c + n] = value;
                //printf("%lf\t", h_c[m*c_c + n]);
        }
        //cout<<"row "<<m<<" completed"<<endl;
    }
}
//////////////////////////////////////////////////////////////////////////
// Option #3
for(m=0; m<r_c; m++)
{
    for(n=0; n<c_c; n++)
    {
#pragma omp parallel
            {
            double value = 0.0;
#pragma omp for
            for(j=0; j<r_b; j++)
                for(k=0; k<c_b; k++)
                {
                    double a;
                    if((m-j)<0 || (n-k)<0 || (m-j)>r_a || (n-k)>c_a)
                        a = 0.0;
                    else
                        a = h_a[((m-j)*c_a) + (n-k)];
                    //printf("%lf\t", a);
                    value += h_b[(j*c_b) + k] * a;
                }
            h_c[m*c_c + n] = value;
            //printf("%lf\t", h_c[m*c_c + n]);
        }
    }
    //cout<<"row "<<m<<" completed"<<endl;
}

测试和配置每个。如果每个线程的工作量不大，您可能会发现选项 #1 最快，或者您可能会发现启用优化后，启用 OMP 时没有区别（甚至减速）。

编辑

我采用了 cmets 中提供的 MCVE，如下所示：

#include <iostream>
#include <chrono>
#include <omp.h>
#include <algorithm>
#include <vector>

#define W_OMP
int main(int argc, char *argv[])
{
    std::vector<double> h_a(9);
    std::generate(h_a.begin(), h_a.end(), std::rand);
    int r_b = 500;
    int c_b = r_b;
    std::vector<double> h_b(r_b * c_b);
    std::generate(h_b.begin(), h_b.end(), std::rand);
    int r_c = 500;
    int c_c = r_c;
    int r_a = 3, c_a = 3;
    std::vector<double> h_c(r_c * c_c);

    auto start = std::chrono::system_clock::now();

#ifdef W_OMP
#pragma omp parallel 
    {
#endif
        int m,n,j,k;
#ifdef W_OMP
#pragma omp for 
#endif
        for(m=0; m<r_c; m++)
        {
            for(n=0; n<c_c; n++)
            {
                double value = 0.0,a;
                for(j=0; j<r_b; j++)
                {
                    for(k=0; k<c_b; k++)
                    {
                        if((m-j)<0 || (n-k)<0 || (m-j)>r_a || (n-k)>c_a)
                            a = 0.0;
                        else a = h_a[((m-j)*c_a) + (n-k)];
                        value += h_b[(j*c_b) + k] * a;
                    }
                }
                h_c[m*c_c + n] = value;
            }
        }
#ifdef W_OMP
    }
#endif
    auto end = std::chrono::system_clock::now();
    auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
    std::cout << elapsed.count() << "ms"
#ifdef W_OMP
        "\t with OMP"
#else
        "\t without OMP"
#endif
        "\n";

    return 0;

}

作为参考，我使用的是 VS2012（OMP 2.0，grrr）。我不确定何时引入 collapse，但显然是在 2.0 之后。优化是 /O2 并在 Release x64 中编译。

基准

使用循环的原始大小 (7,7,5,5) 和数组，结果是 0ms 不带 OMP 和 1ms 带。结论：优化更好，增加的开销是不值得的。此外，测量结果不可靠（太短）。

使用稍大的循环（100、100、100、100）和数组，结果大致相等，约为108ms。结论：仍然不值得天真的努力，调整 OMP 参数可能会影响规模。绝对不是我希望的 x4 加速。

使用更大尺寸的循环（500、500、500、500）和数组，OMP 开始领先。没有 OMP 74.3ms，有 15s。结论：值得。诡异的。我在 i5 上通过四个线程和四个内核获得了 x5 加速。我不会试图弄清楚这是怎么发生的。

总结

正如在 SO 上的无数答案中所述，并行化您遇到的每个 for 循环并不总是一个好主意。可能会破坏您想要的 xN 加速的事情：

1. 每个线程的工作量不足以证明创建额外线程的开销是合理的
2. 工作本身受内存限制。这意味着 CPU 可以以 1petaHz 运行，但您仍然看不到加速。
3. 内存访问模式。我不会去那里。如果需要，请随时在相关信息中进行编辑。
4. OMP 参数。最佳参数选择通常是整个列表的结果（不包括此项，以避免递归问题）。
5. SIMD 操作。根据您的操作和操作方式，编译器可能会向量化您的操作。我不知道 OMP 是否会篡夺 SIMD 操作，但这是可能的。检查您的程序集（对我来说是外语）以确认。