Réponses différentes en virgule flottante avec un nombre différent de processus

Question

Je suis nouveau dans l'apprentissage de MPI et j'ai codé le programme simple suivant pour effectuer l'intégration. avec la règle trapézoïdale en utilisant Open MPI sur Ubuntu 10.10. Voici le code :

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

//function to integrate
double f (double x )
{
  return 4.0/(1+x*x);
}

//function which integrates the function defined above on the interval local_a and local_b for a given refinement parameters
double Trap(double local_a , double local_b, int local_n , double h)
{

  double integral ;
  double x;

  integral = ( f(local_a) + f(local_b) )/2.0;

  x = local_a ;

  for (int i = 1; i < local_n - 1; ++i)
    {
      x        += h; 
      integral += f(x);
    }

  integral *= h;
  return integral;
}

int main(int argc, char *argv[])
{

  int my_rank;
  int p;
  double a = 0.0;
  double b = 1.0;
  int n = atoi(argv[1]);//number of subdivisions of the interval
  double h;
  double local_a;
  double local_b;
  int local_n;

  double integral;
  double total;
  int source;
  int dest = 0;
  int tag = 0;
  MPI_Status status;

  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD,&p);//get number pf processes
  MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);//get rank

  double start , finish;
  MPI_Barrier(MPI_COMM_WORLD);
  start = MPI_Wtime();  

////////////////////////////////////////////////////////////////////////////////////////////////////
  h = (b-a)/n;
  local_n = n/p;

  local_a = a + my_rank*local_n*h;
  local_b = local_a + local_n*h;

  integral = Trap(local_a , local_b , local_n , h);

  if (my_rank==0)
    {
      total = integral;

     for (source = 1; source < p; ++source)
       {
     MPI_Recv(&integral, 1, MPI_DOUBLE , source , tag , MPI_COMM_WORLD, &status );
         total+= integral;

       }
     }

  else
    {
      MPI_Send(&integral, 1, MPI_DOUBLE, dest, tag , MPI_COMM_WORLD);
    }

  if (my_rank == 0)
    {
      printf("With n=%d trapezoids our estimate \n", n );
      printf("Of the integral from %f to %f = %f \n" , a ,b , total);

    }

   ////////////////////////////////////////////////////////////////////////////////////////////////////
  MPI_Barrier(MPI_COMM_WORLD);
  finish = MPI_Wtime();

  if(my_rank == 0)  std::cout << "Time taken is " << finish - start << std::endl ; 

  MPI_Finalize();
  return 0;
}

La fonction intégrée est f(x) = 4.0 / 1+x^2 qui, lorsqu'il est intégré sur [0,1] donne pi = 3.14159...

Maintenant, lorsque j'ai exécuté le programme avec un nombre différent de processus, j'obtiens des réponses différentes. Et la différence est assez significative comme vous pouvez le voir ci-dessous.

Desktop: mpirun -np 1 ./a.out 50000
With n=50000 trapezoids our estimate 
Of the integral from 0.000000 to 1.000000 = 3.141553 
Time taken is 0.000718832
Desktop: 
Desktop: 
Desktop: mpirun -np 2 ./a.out 50000
With n=50000 trapezoids our estimate 
Of the integral from 0.000000 to 1.000000 = 3.141489 
Time taken is 0.000422001
Desktop: 
Desktop: 
Desktop: 
Desktop: mpirun -np 3 ./a.out 50000
With n=50000 trapezoids our estimate 
Of the integral from 0.000000 to 1.000000 = 3.141345 
Time taken is 0.000365019
Desktop: 
Desktop: 
Desktop: 
Desktop: mpirun -np 4 ./a.out 50000
With n=50000 trapezoids our estimate 
Of the integral from 0.000000 to 1.000000 = 3.141362 
Time taken is 0.0395319

Answer 1

Vous avez deux problèmes différents dans votre code :

1. Les limites d'intégration dépendent du nombre de processus MPI, et sont fausses quand p ne divise pas n . A savoir, la borne supérieure du dernier processus est

a + p * int(n/p) * (b-a)/n

qui est différent de b . Je pense que c'est l'erreur la plus importante de votre code (sauf s'il y a un autre bug que je n'ai pas vu).

2. Les opérations en virgule flottante ne sont ni associatives ni commutatives. Le résultat de votre algorithme parallèle, qui est agrégé à partir de sommes partielles, dépendra donc du nombre de sommes partielles.

Answer 2

Lorsque l'on fait de l'arithmétique à virgule flottante, l'ordre des opérations est important. En arithmétique réelle (je veux dire l'arithmétique des nombres réels) a+b+c+d==a+c+d+b (et tout autre ordre des ajouts). Ceci n'est pas nécessairement vrai pour l'arithmétique à virgule flottante. Puisque MPI ne garantit pas que les réductions de M processeurs à 1 processeur seront effectuées dans le même ordre à chaque fois, son comportement en virgule flottante est non-déterministe, du moins pour la plupart d'entre nous.

En mettant cela de côté, les différences entre les résultats sur le nombre variable de processeurs semblent plutôt plus importantes que ce à quoi je m'attendais. En regardant votre code, je pense que vous utilisez l'arithmétique des entiers dans cette ligne :

local_n = n/p;

ce qui fait que de petites parties de la surface totale ne sont attribuées à aucun des processus pour le calcul. D'après mes calculs, la ligne

local_b = local_a + local_n*h;

ne fixe pas local_b à 1,0 pour le dernier processus.