Различные ответы с плавающей запятой с разным количеством процессов - PullRequest
2 голосов
/ 21 марта 2012

Я новичок в изучении MPI и написал следующую простую программу для интеграции с трапециевидным правилом, использующим Open MPI в Ubuntu 10.10. Вот код:

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

Интегрируемая функция - f(x) = 4.0 / 1+x^2, которая при интеграции на [0,1] дает pi = 3.14159...

Теперь, когда я запустил программу с разным количеством процессов, я получил разные ответы. И разница довольно значительна, как вы можете видеть ниже.

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

Ответы [ 2 ]

4 голосов
/ 21 марта 2012

В вашем коде 2 разные проблемы:

1. Границы интеграции зависят от количества процессов MPI и неверны, когда p не делит n. А именно, верхняя граница последнего процесса равна

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

, что отличается от b. Я ожидаю, что это самая важная ошибка в вашем коде (кроме случаев, когда есть еще одна ошибка, которую я не видел)

2. Операции с плавающей точкой не являются ни ассоциативными, ни коммутативными. Таким образом, результат вашего параллельного алгоритма, который агрегируется из частичных сумм, будет зависеть от количества частичных сумм.

2 голосов
/ 21 марта 2012

При выполнении арифметики с плавающей точкой порядок операций имеет значение.В реальной арифметике (я имею в виду арифметику действительных чисел) a+b+c+d==a+c+d+b (и любой другой порядок сложений).Это не обязательно верно для арифметики с плавающей точкой.Поскольку MPI не гарантирует сокращение M процессоров до 1 процессора в одном и том же порядке каждый раз, когда его поведение с плавающей запятой недетерминировано, по крайней мере, насколько это необходимо большинству из нас.

что, с одной стороны, различия между результатами на различном числе процессоров выглядят довольно большими, чем я ожидал.Глядя на ваш код, я думаю, что вы используете целочисленную арифметику в этой строке:

local_n = n/p;

, что приводит к тому, что небольшие части общей площади не назначаются ни одному из процессов для расчета.По моим расчетам, строка

local_b = local_a + local_n*h;

не устанавливает для local_b значение 1.0 для последнего процесса.

...