Call function when running MPI in C++

Question

I wanted to try to use OpenMPI with C++, so I wrote a small code to do numerical integration. My problem is that it deos not seem to execute the line where it all happens correctly:

integral = trapezintegration(local_a, local_b, local_n);

Now I am quite confident that the MPI works correctly beside this line. When printing out the local_a, local_b, local_n and rank_world, I get:

0 3.75 2.5e+09 0  
3.75 7.5 2.5e+09 1   
7.5 11.25 2.5e+09 2  
11.25 15 2.5e+09 3 

Which is excactly what I expected. When I print out integral, rank_world. I get:

17.5781 2  
17.5781 3  
17.5781 1  
17.5781 0    

This is the part that seems wierd to me, only rank_world = = 0, should have the value of integral = 17.5781. My question is, how do I make a function call in MPI, so the ranks dosen't all get the value of rank_world==0?

The complete code can be seen below:

#include <mpi.h>
#include <iostream>

double f(const double x){
    return x*x;
}

double trapezintegration(const double a, const double b, const double n){
    "a = start value of integration range";
    "b = end value of integration range";
    "n = number of integration slices";
    double integral=0.0, h=(b-a)/n;
    long loopbound = (long)n;
    "integral = the value of the numeric integral";
    "h        = width of the numeric integration";
    integral = -(f(a)+f(b))/2.0;
    for (long i=1;i<=loopbound;i++){
        integral = integral + f(i*h);
    }
    integral = integral*(b-a)/n;
    return integral;
}


int main(){
    // The MPI enviroment need to be initialized
    MPI_Init(NULL, NULL);

    // The program need to know how many processors that are avaible
    int world_size;
    MPI_Comm_size(MPI_COMM_WORLD, &world_size);

    // The processors index is also needed to be known
    int world_rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &world_rank);

    // Now the execution of the program can be done

    // If no processor index (rank) is specfied the code will be 
    // executed for all the processes

    const double a=0.0, b=15.0, n=1e+10;
    double integral, total_integral;

    // Right now all of the processes have the same a and b
    // now different a and b will be assigned to the processes

    // The rank is the index of the processor going from
    // 0 to WORLD_SIZE-1, all of the processes will now get
    // different local_a and local_b

    double local_a = (b - a)/world_size*world_rank;
    double local_b = (b - a)/world_size*(world_rank+1);
    double local_n = n/world_size;

    std::cout << local_a << ' '<< local_b << ' ' << local_n << ' '  << world_rank << '\n'; 
    integral = trapezintegration(local_a, local_b, local_n);

    // All of the processes have now run the numerical integration
    // for their given interval. All of the integrated parts need
    // to be collected to get the total integration.
    // Lets collect the result in Rank 0
    std::cout << integral << ' ' << world_rank << '\n'; 

    if (world_rank != 0){
        MPI_Send(&integral,1,MPI_DOUBLE,0,555+world_rank,MPI_COMM_WORLD);
    }
    if (world_rank == 0){
        total_integral = integral;
        for (int i=1; i<world_size; i++){
                MPI_Recv(&integral,1,MPI_DOUBLE,i,555+i,MPI_COMM_WORLD,MPI_STATUS_IGNORE);
                total_integral = total_integral + integral;
        }

    }

    // if rank is different from rank 0, the result need to be send
    if (world_rank == 0){
        std::cout << total_integral << '\n';
    }

    // The MPI enviroment need to be closed when the calculation is finished
    MPI_Finalize();
}

Answer 1

Instead of

integral = integral + f(i*h);

Should it be

integral = integral + f(a+i*h);

Instead ?

As a side note, the natural way of computing total_integral is via a single call to MPI_Reduce()