並行計算學習之用MPI實現梯形積分法

• 梯形積分法

基本思想是，將x軸上區間划分成n個等長的子區間。估計介於函數圖像以及每個子區間內梯形區域的面積。

　　設子區間端點為x_i和x_i+1 ，長度h=x_i+1 - x_i， 同樣的兩條垂直線的長度為f(x_i)和f(x_i+1)

　　　　那么面積為：h/2[ f(x_i) + f(x_i+1) ]

n個區間是等分的，如果兩條垂直線包圍區域的邊界分別為a和b，那么

　　　　　　h = (b-a) / n

• 設計並行程序的四個步驟

1. 將問題的解決方案划分成多個任務
2. 在任務間識別出需要的通信信道（識別出聯系，是否需要交流）
3. 將任務聚合成復合任務 --->減少通信
4. 在核上分配復合任務

• 串行代碼

/*  Input  a ,b , n*/
h = (b-a)/n
approx  = (f(a) + f(b))/2.0

for(i =1; i <= n-1; i++){
    x_i = a + i * h;
    approx = h * approx
}

approx = h*approx

• 　並行程序

     偽代碼：

Get a,b,n;
h = (b-a)/n;
local_n = n/comm_sz;
local_a = a + my_rank * local_n * h;
local_b = local_a + local_n*h;
local_integral = Trap(local_a, local_b, local_n, h)  /*Trap函數作用*/

if( my_rank != 0 )
     Send local_integral to process 0;

else /*my_rank == 0*/
      total_integral = local_integral;
      for(proc = 1;proc < comm_sz; proc++){
           Receive local_integral from proc;
           total_integral += local_integral;
     }

if(my_rank == 0)
      print result;

• 　第一個MPI代碼----輸出

int main(void){
     int my_rank , comm_sz, n = 1024, local_n;
     double a = 0.0, b=3.0, h, local_a,local_b;
     double local_int,total_int;
     int source;    

     MPI_Init(NULL,NULL);
     MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
     MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
   
     h = (b-a)/n;
     local_n = n/comm_sz;
     
     local_a = a + my_rank * local_n * h;
     local_b = local_a +local_n * h;
     local_int = Trap(local_a , local_b , local_n, h )

    if (my_rank != 0){
            MPI_Send(&local_int , 1 , MPI_DOUBLE , 0, 0 , MPI_COMM_WORLD);
    }else{
             total_int = local_int;
             for(source = 1 ; source < comm_sz ; source++){
                    MPI_Recv(&local_int, 1, MPI_DOUBLE, source, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                    total_int += local_int;
             }
    }

   if (my_rank == 0 ){
           printf("With n = %d trapezoids, our estimate\n",n);
           printf("of the integral from %f to %f = %.15e\n",a,b,total_int);
   }

   MPI_Finalize();
   return 0;
}                    

　　其中Trap函數

double  Trap(
            double    left_endpt .
            double    right_endpt,
            int          trap_count,
            double    base_len
){
            double estimate , x;
            int i;
             
            estimate = (f(left_endpt) + f(right_endpt))/2.0
            for(i = 1;i <= trap_count-1 ; i++){
                        x = left_endpt + i*base_len;
                        estimate += f(x);
            }
            estimate = estimate * base_len;
            
             return estimate;
}

• 　　輸入

void Get_input(
       int          my_rank,
       int          comm_sz,
       double*   a_p,
       double*   b_p,
       int*         n_p
){
        int dest;
         
        if(my_rank == 0){
            printf("Enter a, b,and n \n");
            scanf("%lf  %lf  %d",a_p,b_p,n_p);
            for(dest = 1; dest < comm_sz;dest++){
                MPI_Send(a_p , 1 , MPI_DOUBLE, dest , 0, MPI_COMM_WORLD);
                MPI_Send(b_p , 1 , MPI_DOUBLE, dest , 0, MPI_COMM_WORLD);
                MPI_Send(n_p , 1 , MPI_INT, dest , 0, MPI_COMM_WORLD);
                 }

        }else{
                 MPI_Recv(a_p , 1, MPI_DOUBLE， 0， 0，MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                 MPI_Recv(b_p , 1, MPI_DOUBLE， 0， 0，MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                 MPI_Recv(n_p , 1, MPI_INT， 0， 0，MPI_COMM_WORLD, MPI_STATUS_IGNORE);
                }
}

• 　　樹形通信