如何使用 OpenMP 以正确的方式并行化此数组？

Question

我尝试用openmp并行代码后，数组中的元素是错误的，至于元素的顺序不是很重要。或者使用 c++ std vector 而不是 array 来并行化更方便，你能推荐一个简单的方法吗？

#include <stdio.h>
#include <math.h>

int main()
{
    int n = 100;
    int a[n*(n+1)/2]={0};
    int count=0;

    #pragma omp parallel for reduction(+:a,count)
    for (int i = 1; i <= n; i++) {
        for (int j = i + 1; j <= n; j++) {
            double k = sqrt(i * i + j * j);
            if (fabs(round(k) - k) < 1e-10) {
                a[count++] = i;
                a[count++] = j;
                a[count++] = (int) k;
            }
        }
    }
    
    for(int i=0;i<count;i++)
        printf("%d %s",a[i],(i+1)%3?"":", ");
    printf("\ncount: %d", count);
    return 0;
}

原始输出：

3 4 5 , 5 12 13 , 6 8 10 , 7 24 25 , 8 15 17 , 9 12 15 , 9 40 41 , 10 24 26 , 11 60 61 , 12 16 20 , 12 35 37 , 13 84 85 , 14 48 50 , 15 20 25 , 15 36 39 , 16 30 34 , 16 63 65 , 18 24 30 , 18 80 82 , 20 21 29 , 20 48 52 , 20 99 101 , 21 28 35 , 21 72 75 , 24 32 40 , 24 45 51 , 24 70 74 , 25 60 65 , 27 36 45 , 28 45 53 , 28 96 100 , 30 40 50 , 30 72 78 , 32 60 68 , 33 44 55 , 33 56 65 , 35 84 91 , 36 48 60 , 36 77 85 , 39 52 65 , 39 80 89 , 40 42 58 , 40 75 85 , 40 96 104 , 42 56 70 , 45 60 75 , 48 55 73 , 48 64 80 , 48 90 102 , 51 68 85 , 54 72 90 , 56 90 106 , 57 76 95 , 60 63 87 , 60 80 100 , 60 91 109 , 63 84 105 , 65 72 97 , 66 88 110 , 69 92 115 , 72 96 120 , 75 100 125 , 80 84 116 ,
count: 189

使用 OpenMP(gcc file.c -fopenmp) 后:

411 538 679 , 344 609 711 , 354 533 649 , 218 387 449 , 225 475 534 , 182 283 339 , 81 161 182 , 74 190 204 , 77 138 159 , 79 176 195 , 18 24 30 , 18 80 82 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 , 0 0 0 ,
count: 189

Answer 1

您的线程都在访问共享 count。

您最好消除 count 并让每个循环迭代仅根据 i 和 j 的（每个线程）值来确定将其输出写入何处。

或者，使用向量来累积结果：

#include <cmath>
#include <iostream>
#include <utility>
#include <vector>

#pragma omp declare                                                        \
    reduction(vec_append : std::vector<std::pair<int,int>> :               \
              omp_out.insert(omp_out.end(), omp_in.begin(), omp_in.end()))

int main()
{
    constexpr int n = 100'000;
    std::vector<std::pair<int,int>> result;

#pragma omp parallel for \
            reduction(vec_append:result) \
            schedule(dynamic)
    for (int i = 1;  i <= n;  ++i) {
        for (int j = i + 1;  j <= n;  ++j) {
            auto const h2 = i * i + j * j; // hypotenuse squared
            int h = std::sqrt(h2) + 0.5;   // integer square root
            if (h * h == h2) {
                result.emplace_back(i, j);
            }
        }
    }

    // for (auto const& v: result) {
    //     std::cout << v.first << ' '
    //               << v.second << ' '
    //               << std::hypot(v.first, v.second) << ' ';
    // }
    std::cout << "\ncount: " << result.size() << '\n';
}

Answer 2

count 变量是 a 的索引。 reduction(+:a,count) 运算符正在对数组求和，这不是我认为您正在寻找的串联运算。

计数变量需要用互斥量包围，类似于#pragma omp critical，但我不是 OpenMP 专家。

或者，创建 int a[n][n]，将它们全部设置为 -1（表示“无效”的标记值），然后在 sqrt() 的结果足够接近一个整体时分配它数.

Answer 3

作为使用 critical 部分的替代方案，此解决方案使用原子，因此可能更快。

由于内存消耗，以下代码可能会冻结您的计算机。小心！

#include <cstdio>
#include <cmath>

#include <vector>

int main() {
    int const n = 100;
    // without a better (smaller) upper_bound this is extremely
    // wasteful in terms of memory for big n 
    long const upper_bound = 3L * static_cast<long>(n) *
                             (static_cast<long>(n) - 1L) / 2l; 
    std::vector<int> a(upper_bound, 0);
    int count = 0;

    #pragma omp parallel for schedule(dynamic) shared(a, count)
    for (int i = 1; i <= n; ++i) {
        for (int j = i + 1; j <= n; ++j) {
            double const k = std::sqrt(static_cast<double>(i * i + j * j));

            if (std::fabs(std::round(k) - k) < 1e-10) {
                int my_pos;
                #pragma omp atomic capture
                my_pos = count++;

                a[3 * my_pos] = i;
                a[3 * my_pos + 1] = j;
                a[3 * my_pos + 2] = static_cast<int>(std::round(k));
            }
        }
    }
    count *= 3;
    
    for(int i = 0; i < count; ++i) {
        std::printf("%d %s", a[i], (i + 1) % 3 ? "" : ", ");
    }
    printf("\ncount: %d", count);

    return 0;
}

编辑： 我的回答最初是以次优方式使用 critical 部分对现已删除的答案的反应。在下文中，我将介绍另一种解决方案，它将 critical 部分与使用 std::vector::emplace_back() 相结合，以规避对 upper_bound 的需求，类似于 Toby Speight 的解决方案。一般来说，像 Toby Speight 的解决方案中那样使用 reduce 子句应该优于 critical 部分和 atomics，因为缩减对于大量线程应该更好地扩展。在这种特殊情况下（相对较少的计算将写入 a）并且没有大量的核心到运行，以下代码可能仍然更可取。

#include <cstdio>
#include <cmath>

#include <tuple>
#include <vector>

int main() {
    int const n = 100;

    std::vector<std::tuple<int, int, int>> a{};
    
    // optional, might reduce number of reallocations
    a.reserve(2 * n); // 2 * n is an arbitrary choice

    #pragma omp parallel for schedule(dynamic) shared(a)
    for (int i = 1; i <= n; ++i) {
        for (int j = i + 1; j <= n; ++j) {
            double const k = std::sqrt(static_cast<double>(i * i + j * j));

            if (std::fabs(std::round(k) - k) < 1e-10) {
                #pragma omp critical
                a.emplace_back(i, j, static_cast<int>(std::round(k)));
            }
        }
    }
    long const count = 3L * static_cast<long>(a.size());
    
    for(unsigned long i = 0UL; i < a.size(); ++i) {
        std::printf("%d %d %d\n",
                    std::get<0>(a[i]), std::get<1>(a[i]), std::get<2>(a[i]));
    }
    printf("\ncount: %ld", count);

    return 0;
}

如何使用 OpenMP 以正确的方式并行化此数组？

How to parallelize this array correct way using OpenMP?

c++

openmp