1

This is my solution using openmp which i used to parallelize the code which calculates Pi. The floating point value of Pi changes each time this is executed. Could someone explain why?

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

#define THREAD_NUM 20

static long num_steps = 100000;
double step;

int main(){

    int i;
    double x; 
    double pi;
    double sum = 0.0;
    double t1 = 0.0;
    double t2 = 0.0;

    step = 1.0/(double) num_steps;
    omp_set_num_threads(THREAD_NUM);
    t1 = omp_get_wtime();

    #pragma omp parallel
    {
        double p_sum = 0.0;

        #pragma omp for
        for(i=0; i<num_steps; i++){
    
            x = (i+0.5)*step;
            p_sum = p_sum + 4.0/(1.0+x*x);
        }

        #pragma omp atomic
        sum += p_sum;
    }

    t2 = omp_get_wtime();
    pi = step*sum;

    printf("value of pi = %lf\n", pi);
    printf("time = %lf ms\n", (t2-t1)*1000);
}
vidath007
  • 21
  • 4
  • 1
    I haven't checked your code closely, but this - https://stackoverflow.com/questions/20993327/using-openmp-critical-and-ordered/20996615 - explains why many parallel numeric computations differ in (usually) insignificant digits from run to run. Oh, and I will comment that you seem to have rolled your own *reduction*. Is that because you haven't learned yet about OpenMP's inbuilt reductions ? I don't see your use of a home-grown variant as a likely source of the issue you report. – High Performance Mark Sep 02 '20 at 08:27
  • 1
    Your `x` variable needs to be declared `private` – Gilles Sep 02 '20 at 08:48
  • @Gilles ok yeah, that was the issue. Thank you so much, everything works fine now. If you can provide this as an answer i can mark it as the answer? – vidath007 Sep 02 '20 at 09:31
  • @HighPerformanceMark Thank you, no I didn't know about reductions till you mentioned it. I am new to OpenMP. You mean using, "#pragma omp parallel reduction (+: sum)" instead of an atomic operation right? – vidath007 Sep 02 '20 at 10:15
  • Yes, a reduction instead of a summation guarded by an atomic directive. – High Performance Mark Sep 02 '20 at 11:15
  • @HighPerformanceMark Could you recommend any sources to learn OpenMP for beginners? – vidath007 Sep 02 '20 at 11:27
  • No, that's not allowed on SO, I might have my bagel privileges terminated with extreme prejudice. But both your favourite search engine, and a well-known online bookseller will be able to help you with your query. – High Performance Mark Sep 02 '20 at 12:11
  • @HighPerformanceMark oh okay, sorry i didn't know that. Haha sure. Thank you so much for your time and support. – vidath007 Sep 02 '20 at 13:45

2 Answers2

3

Floating point addition is neither associative nor commutative! This means that the exact value you obtain depends on the order in which the components of p_sum/sum are added up. To understand precisely why you have to understand how floating point addition works in practice. I would recommend reading What Every Computer Scientist should Know About Floating-Point Arithmetic.

Peter
  • 2,919
  • 1
  • 16
  • 35
  • 1
    I would only recommend *What every computer scientist ...* to people developing f-p arithmetic libraries and suchlike. For most programmers *using* the f-p capabilities of their platforms *What Every Programmer Should Know About Floating-Point Arithmetic* (https://floating-point-gui.de) is a much better place to start learning. – High Performance Mark Sep 02 '20 at 08:34
  • @HighPerformanceMark: I didn't even know the latter existed, I just remembered the former from University. That does indeed look more suitable for beginners. – Peter Sep 02 '20 at 08:51
1

As @Gilles pointed out in the comments under the question, the problem was with the x variable which was declared as a shared variable. it should be declared as a private variable.

...

#pragma omp parallel
    {
        double x;
        double p_sum = 0.0;

        #pragma omp for
        for(i=0; i<num_steps; i++){

...
vidath007
  • 21
  • 4