Numbers not randomized after runs

Question

I'm trying to create an openMP program that randomizes double arrays and run the values through the formula: y[i] = (a[i] * b[i]) + c[i] + (d[i] * e[i]) + (f[i] / 2);

If I run the program multiple times I've realised that the Y[] values are the same even though they are supposed to be randomized when the arrays are initialized in the first #pragma omp for . Any Ideas as to why this might be happening?

#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;
}

Answer 1

1. First of all, I would like to emphasize that OpenMP has significant overheads, so you need a reasonable amount of work in your code, otherwise the overhead is bigger than the gain by parallelization. In your code this is the case, so the fastest solution is to use serial code. However, you mentioned that your goal is to learn OpenMP, so I will show you how to do it.

2. In your previous post's comments @paleonix linked a post ( How to generate random numbers in parallel? ) which answers your question about random numbers. One of the solutions is to use rand_r.

3. Your code has a data race when searching for minimum and maximum values of array Y. If you need to find the minimum/maximum value only it is very easy, because you can use reduction like this:

double max=y[0];
#pragma omp parallel for default(none) shared(y) reduction(max:max) 
for (int i = 1; i < ARRAY_SIZE; i++) {  
    if (y[i] > max) {
        max = y[i];
    }
}

But in your case you also need the indices of minimum and maximum value, so it is a bit more complicated. You have to use a critical section to be sure that other threads can not change the max, min, imax and imin values while you updating their values. So, it can be done the following way (e.g. for finding minimum value):

 #pragma omp parallel for
 for (int i = 0; i < ARRAY_SIZE; i++) {  
    if (y[i] < min) {
        #pragma omp critical
        if (y[i] < min) {
            min = y[i];
            imin = i;
        }
    }
 }

Note that the if (y[i] < min) appears twice, because after the first comparison other threads may change the value of min, so inside the critical region before updating min and imin values you have to check it again. You can do it exactly the same way in the case of finding the maximum value.

3. Always use your variables at their minimum required scope.

4. It is also recommend to use default(none) clause in your OpenMP parallel region so, you have to explicitly define the sharing attributes all of your variables.

5. You can fill the array and find its minimum/maximum values in a single loop and print their values in a different serial loop.

6. If you set min and max before the loop, you can get rid of the extra comparison if (i == 0) used inside the loop.

Putting it together:

double threadsafe_rand(unsigned int* seed, double min, double max)
{
    double range = (max - min);
    double div = RAND_MAX / range;
    return min + (rand_r(seed) / div);
}

In main:

double min=DBL_MAX;
double max=-DBL_MAX;

#pragma omp parallel default(none) shared(a,b,c,d,e,f,y,imin,imax,min,max)  
{
    unsigned int seed=omp_get_thread_num();   
    #pragma omp for
    for (int i = 0; i < ARRAY_SIZE; i++) {  
        a[i] = threadsafe_rand(&seed, 1,100);
        b[i] = threadsafe_rand(&seed,10, 50);                    
        c[i] = threadsafe_rand(&seed,1, 10);
        d[i] = threadsafe_rand(&seed,1, 50);
        e[i] = threadsafe_rand(&seed,1, 5);
        f[i] = threadsafe_rand(&seed,10, 80);
        y[i] = (a[i] * b[i]) + c[i] + (d[i] * e[i]) + (f[i] / 2);

        if (y[i] < min) {
            #pragma omp critical
            if (y[i] < min) {
                min = y[i];
                imin = i;
            }
        }

        if (y[i] > max) {
            #pragma omp critical
            if (y[i] > max) {
                max = y[i];
                imax = i;
            }
        }
    }
}

// printout 
for (int i = 0; i < ARRAY_SIZE; i++) {
    printf("Y[%d]=%.2f\n", i, y[i]);
}
printf("min y[%d] = %.2f\nmax y[%d] = %.2f\n", imin, min, imax, max);

Update: I have updated the code according to @Qubit's and @JérômeRichard's suggestions:

1. I used the 'Really minimal PCG32 code' / (c) 2014 M.E. O'Neill / from https://www.pcg-random.org/download.html. Note that I do not intend to properly handle the seeding of this simple random number generator. If you would like to do so, please use a complete random number generator library.

2. I have changed the code to use user defined reductions. Indeed, it makes the code much more efficient, but not really beginner friendly. It would require a very long post to explain it, so if you are interested in the details, please read a book about OpenMP.

3. I have reduced the number of divisions in threadsafe_rand

The updated code:

#include<stdio.h>
#include<stdint.h>
#include<time.h>
#include<float.h>
#include<limits.h>
#include<omp.h>
#define ARRAY_SIZE 10

// *Really* minimal PCG32 code / (c) 2014 M.E. O'Neill / pcg-random.org
// Licensed under Apache License 2.0 (NO WARRANTY, etc. see website)

typedef struct { uint64_t state;  uint64_t inc; } pcg32_random_t;

inline uint32_t pcg32_random_r(pcg32_random_t* rng)
{
    uint64_t oldstate = rng->state;
    // Advance internal state
    rng->state = oldstate * 6364136223846793005ULL + (rng->inc|1);
    // Calculate output function (XSH RR), uses old state for max ILP
    uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u;
    uint32_t rot = oldstate >> 59u;
    return (xorshifted >> rot) | (xorshifted << ((-rot) & 31));
}

inline double threadsafe_rand(pcg32_random_t* seed, double min, double max)
{
    const double tmp=1.0/UINT32_MAX;  
    return min + tmp*(max - min)*pcg32_random_r(seed);
}

struct v{
 double value;
 int i;   
};

#pragma omp declare reduction(custom_min: struct v: \
 omp_out = omp_in.value < omp_out.value ? omp_in : omp_out )\
 initializer(omp_priv={DBL_MAX,0} )

#pragma omp declare reduction(custom_max: struct v: \
 omp_out = omp_in.value > omp_out.value ? omp_in : omp_out )\
 initializer(omp_priv={-DBL_MAX,0} )

int main() {
    double a[ARRAY_SIZE], b[ARRAY_SIZE], c[ARRAY_SIZE], d[ARRAY_SIZE], e[ARRAY_SIZE], f[ARRAY_SIZE], y[ARRAY_SIZE];
    struct v max={-DBL_MAX,0};
    struct v min={DBL_MAX,0}; 

    #pragma omp parallel default(none) shared(a,b,c,d,e,f,y) reduction(custom_min:min) reduction(custom_max:max)
    {
        pcg32_random_t seed={omp_get_thread_num()*7842 + time(NULL)%2299, 1234+omp_get_thread_num()};   
        #pragma omp for
        for (int i=0 ; i < ARRAY_SIZE; i++) {  
            a[i] = threadsafe_rand(&seed, 1,100);
            b[i] = threadsafe_rand(&seed,10, 50);                    
            c[i] = threadsafe_rand(&seed,1, 10);
            d[i] = threadsafe_rand(&seed,1, 50);
            e[i] = threadsafe_rand(&seed,1, 5);
            f[i] = threadsafe_rand(&seed,10, 80);
            y[i] = (a[i] * b[i]) + c[i] + (d[i] * e[i]) + (f[i] / 2);

            if (y[i] < min.value) {
                    min.value = y[i];
                    min.i = i;
                }

            if (y[i] > max.value) {
                    max.value = y[i];
                    max.i = i;
                }
        }
    }

    // printout 
    for (int i = 0; i < ARRAY_SIZE; i++) {
        printf("Y[%d]=%.2f\n", i, y[i]);
    }
    printf("min y[%d] = %.2f\nmax y[%d] = %.2f\n", min.i, min.value, max.i, max.value);
    return 0;
}