运行后数字未随机化

Ibr*_*him 2 c openmp

我正在尝试创建一个 openMP 程序,该程序随机化双精度数组并通过以下公式运行值:y[i] = (a[i] * b[i]) + c[i] + (d[i] * e[i]) + (f[i] / 2);

如果我多次运行该程序,我就会意识到 Y[] 值是相同的,即使在第一次初始化数组时它们应该是随机的#pragma omp for。关于为什么会发生这种情况有什么想法吗?

#include<stdio.h>
#include <stdio.h>
#include <stdlib.h>
#include<omp.h>
#define ARRAY_SIZE 10

double randfrom(double min, double max);

double randfrom(double min, double max)
{
    double range = (max - min);
    double div = RAND_MAX / range;
    return min + (rand() / div);
}

int main() {
    int i;
    double a[ARRAY_SIZE], b[ARRAY_SIZE], c[ARRAY_SIZE], d[ARRAY_SIZE], e[ARRAY_SIZE], f[ARRAY_SIZE], y[ARRAY_SIZE];
    double min, max;
    int imin, imax;







    /*A[10] consists of random number in between 1 and 100
    B[10] consists of random number in between 10 and 50
    C[10] consists of random number in between 1 and 10
    D[10] consists of random number in between 1 and 50
    E[10] consists of random number in between 1 and 5
    F[10] consists of random number in between 10 and 80*/

    srand(time(NULL));

#pragma omp parallel 

    {

#pragma omp parallel for
        for (i = 0; i < ARRAY_SIZE; i++) {
            a[i] = randfrom(1, 100);
            b[i] = randfrom(10, 50);
            c[i] = randfrom(1, 50);
            d[i] = randfrom(1, 50);
            e[i] = randfrom(1, 5);
            f[i] = randfrom(10, 80);
        }
    }




    printf("This is the parallel Print\n\n\n");

#pragma omp parallel shared(a,b,c,d,e,f,y) private(i)
    {
        //Y=(A*B)+C+(D*E)+(F/2)
#pragma omp for schedule(dynamic) nowait
        for (i = 0; i < ARRAY_SIZE; i++) {
            /*printf("A[%d]%.2f",i, a[i]);
            printf("\n\n");
            printf("B[%d]%.2f", i, b[i]);
            printf("\n\n");
            printf("C[%d]%.2f", i, c[i]);
            printf("\n\n");
            printf("D[%d]%.2f", i, d[i]);
            printf("\n\n");
            printf("E[%d]%.2f", i, e[i]);
            printf("\n\n");
            printf("F[%d]%.2f", i, f[i]);
            printf("\n\n");*/
            y[i] = (a[i] * b[i]) + c[i] + (d[i] * e[i]) + (f[i] / 2);
            printf("Y[%d]=%.2f\n", i, y[i]);
        }
    }




#pragma omp parallel shared(y, min,imin,max,imax) private(i)
    {
        //min
#pragma omp for schedule(dynamic) nowait
        for (i = 0; i < ARRAY_SIZE; i++) {
            if (i == 0) {
                min = y[i];
                imin = i;
            }
            else {
                if (y[i] < min) {
                    min = y[i];
                    imin = i;
                }
            }
        }

        //max
#pragma omp for schedule(dynamic) nowait
        for (i = 0; i < ARRAY_SIZE; i++) {
            if (i == 0) {
                max = y[i];
                imax = i;
            }
            else {
                if (y[i] > max) {
                    max = y[i];
                    imax = i;
                }
            }
        }
    }
    printf("min y[%d] = %.2f\nmax y[%d] = %.2f\n", imin, min, imax, max);
    return 0;
}
Run Code Online (Sandbox Code Playgroud)

Lac*_*aci 6

    \n
  1. 首先,我想强调的是,OpenMP 具有显着的开销,因此您需要在代码中进行合理的工作量,否则开销会大于并行化带来的收益。在您的代码中就是这种情况,因此最快的解决方案是使用串行代码。但是,您提到您的目标是学习 OpenMP,所以我将向您展示如何做到这一点。

    \n
  2. \n
  3. 在您上一篇文章的评论中,@paleonix 链接了一篇文章(如何并行生成随机数?),它回答了您有关随机数的问题。解决方案之一是使用rand_r.

    \n
  4. \n
  5. 在搜索 array 的最小值和最大值时,您的代码会出现数据争用Y。如果您需要找到最小值/最大值,那么这非常容易,因为您可以像这样使用归约:

    \n
  6. \n
\n
double max=y[0];\n#pragma omp parallel for default(none) shared(y) reduction(max:max) \nfor (int i = 1; i < ARRAY_SIZE; i++) {  \n    if (y[i] > max) {\n        max = y[i];\n    }\n}\n
Run Code Online (Sandbox Code Playgroud)\n

但在你的情况下,你还需要最小值和最大值的索引,所以它有点复杂。您必须使用临界区来确保其他线程在您更新 、 和 值时无法更改max它们min的值imaximin因此,可以通过以下方式完成(例如寻找最小值):

\n
 #pragma omp parallel for\n for (int i = 0; i < ARRAY_SIZE; i++) {  \n    if (y[i] < min) {\n        #pragma omp critical\n        if (y[i] < min) {\n            min = y[i];\n            imin = i;\n        }\n    }\n }\n
Run Code Online (Sandbox Code Playgroud)\n

请注意,if (y[i] < min)出现了两次,因为第一次比较后其他线程可能会更改 的值,因此在更新和值之前,必须min在临界区域内再次检查它。在查找最大值的情况下,您可以采用完全相同的方法。minimin

\n
    \n
  1. 始终在所需的最小范围内使用变量。

    \n
  2. \n
  3. 还建议default(none)在 OpenMP 并行区域中使用子句,因此您必须显式定义所有变量的共享属性。

    \n
  4. \n
  5. 您可以填充数组并在单个循环中查找其最小/最大值,并在不同的串行循环中打印它们的值。

    \n
  6. \n
  7. 如果在循环之前设置minand max,则可以消除if (i == 0)循环内使用的额外比较。

    \n
  8. \n
\n

把它放在一起:

\n
double threadsafe_rand(unsigned int* seed, double min, double max)\n{\n    double range = (max - min);\n    double div = RAND_MAX / range;\n    return min + (rand_r(seed) / div);\n}\n
Run Code Online (Sandbox Code Playgroud)\n

主要内容:

\n
double min=DBL_MAX;\ndouble max=-DBL_MAX;\n\n#pragma omp parallel default(none) shared(a,b,c,d,e,f,y,imin,imax,min,max)  \n{\n    unsigned int seed=omp_get_thread_num();   \n    #pragma omp for\n    for (int i = 0; i < ARRAY_SIZE; i++) {  \n        a[i] = threadsafe_rand(&seed, 1,100);\n        b[i] = threadsafe_rand(&seed,10, 50);                    \n        c[i] = threadsafe_rand(&seed,1, 10);\n        d[i] = threadsafe_rand(&seed,1, 50);\n        e[i] = threadsafe_rand(&seed,1, 5);\n        f[i] = threadsafe_rand(&seed,10, 80);\n        y[i] = (a[i] * b[i]) + c[i] + (d[i] * e[i]) + (f[i] / 2);\n\n        if (y[i] < min) {\n            #pragma omp critical\n            if (y[i] < min) {\n                min = y[i];\n                imin = i;\n            }\n        }\n\n        if (y[i] > max) {\n            #pragma omp critical\n            if (y[i] > max) {\n                max = y[i];\n                imax = i;\n            }\n        }\n    }\n}\n\n// printout \nfor (int i = 0; i < ARRAY_SIZE; i++) {\n    printf("Y[%d]=%.2f\\n", i, y[i]);\n}\nprintf("min y[%d] = %.2f\\nmax y[%d] = %.2f\\n", imin, min, imax, max);\n\n
Run Code Online (Sandbox Code Playgroud)\n

更新: \n我已经根据@Qubit\'s和@J\xc3\xa9r\xc3\xb4meRichard\的建议更新了代码:

\n
    \n
  1. 我使用了“真正最小的 PCG32 代码”/(c) 2014 ME O\'Neill/来自https://www.pcg-random.org/download.html。请注意,我不打算正确处理这个简单随机数生成器的播种。如果您想这样做,请使用完整的随机数生成器库。

    \n
  2. \n
  3. 我已更改代码以使用用户定义的缩减。事实上,它使代码更加高效,但对初学者来说并不友好。这需要很长的文章来解释它,所以如果你对细节感兴趣,请阅读一本关于 OpenMP 的书。

    \n
  4. \n
  5. 我减少了部门的数量threadsafe_rand

    \n
  6. \n
\n

更新后的代码:

\n
#include<stdio.h>\n#include<stdint.h>\n#include<time.h>\n#include<float.h>\n#include<limits.h>\n#include<omp.h>\n#define ARRAY_SIZE 10\n\n// *Really* minimal PCG32 code / (c) 2014 M.E. O\'Neill / pcg-random.org\n// Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)\n\ntypedef struct { uint64_t state;  uint64_t inc; } pcg32_random_t;\n\ninline uint32_t pcg32_random_r(pcg32_random_t* rng)\n{\n    uint64_t oldstate = rng->state;\n    // Advance internal state\n    rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);\n    // Calculate output function (XSH RR), uses old state for max ILP\n    uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;\n    uint32_t rot = oldstate >> 59u;\n    return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));\n}\n\ninline double threadsafe_rand(pcg32_random_t* seed, double min, double max)\n{\n    const double tmp=1.0/UINT32_MAX;  \n    return min + tmp*(max - min)*pcg32_random_r(seed);\n}\n\nstruct v{\n double value;\n int i;   \n};\n\n#pragma omp declare reduction(custom_min: struct v: \\\n omp_out = omp_in.value < omp_out.value ? omp_in : omp_out )\\\n initializer(omp_priv={DBL_MAX,0} )\n\n#pragma omp declare reduction(custom_max: struct v: \\\n omp_out = omp_in.value > omp_out.value ? omp_in : omp_out )\\\n initializer(omp_priv={-DBL_MAX,0} )\n\nint main() {\n    double a[ARRAY_SIZE], b[ARRAY_SIZE], c[ARRAY_SIZE], d[ARRAY_SIZE], e[ARRAY_SIZE], f[ARRAY_SIZE], y[ARRAY_SIZE];\n    struct v max={-DBL_MAX,0};\n    struct v min={DBL_MAX,0}; \n\n    #pragma omp parallel default(none) shared(a,b,c,d,e,f,y) reduction(custom_min:min) reduction(custom_max:max)\n    {\n        pcg32_random_t seed={omp_get_thread_num()*7842 + time(NULL)%2299, 1234+omp_get_thread_num()};   \n        #pragma omp for\n        for (int i=0 ; i < ARRAY_SIZE; i++) {  \n            a[i] = threadsafe_rand(&seed, 1,100);\n            b[i] = threadsafe_rand(&seed,10, 50);                    \n            c[i] = threadsafe_rand(&seed,1, 10);\n            d[i] = threadsafe_rand(&seed,1, 50);\n            e[i] = threadsafe_rand(&seed,1, 5);\n            f[i] = threadsafe_rand(&seed,10, 80);\n            y[i] = (a[i] * b[i]) + c[i] + (d[i] * e[i]) + (f[i] / 2);\n\n            if (y[i] < min.value) {\n                    min.value = y[i];\n                    min.i = i;\n                }\n\n            if (y[i] > max.value) {\n                    max.value = y[i];\n                    max.i = i;\n                }\n        }\n    }\n\n    // printout \n    for (int i = 0; i < ARRAY_SIZE; i++) {\n        printf("Y[%d]=%.2f\\n", i, y[i]);\n    }\n    printf("min y[%d] = %.2f\\nmax y[%d] = %.2f\\n", min.i, min.value, max.i, max.value);\n    return 0;\n}\n
Run Code Online (Sandbox Code Playgroud)\n

  • 我只是想指出,关键部分仍然不是必需的,而且效率低下。该问题仍然可以通过结构体和自定义归约的方式并行解决。 (3认同)
  • AFAIK `rand_r` 是一个仅限 posix 的函数,而不是标准的 C 函数,但据我所知,没有线程安全的 C 函数。使用自定义解决方案可能是一种解决方案(主要用于 SIMD 优化)。请注意,Windows 上有等效的“rand_s”。在 C++ 中,有完整的 STL 部分(遗憾的是 OP 不使用 C++)。`threadsafe_rand` 效率相当低,因为它使用了 2 个昂贵的除法(如果在此特定示例中内联函数,则可以优化其中一个除法,另一个只能通过 fastmath 进行优化)。为什么不使用“range / RAND_MAX”的产品来代替? (2认同)
  • @JérômeRichard 虽然确实如此,但人们可以轻松下载 PRGN 的 C 代码,例如:https://www.pcg-random.org/download.html 作为示例。另请注意,通常建议这样做,而不是 C++ STL 中的 Mersenne Twister 之类的东西,因为 C++ STL 的变化速度有点慢。 (2认同)