並行歸並排序——MPI

      並行歸並排序在程序開始時，會將n/comm_comm個鍵值分配給每個進程，程序結束時，所有的鍵值會按順序存儲在進程0中。為了做到這點，它使用了樹形結構通信模式。當進程接收到另一個進程的鍵值時，它將該鍵值合並進自己排序的鍵值列表中。編寫一個程序實現歸並排序。進程0應該讀入n的值，將其廣播給其余進程。每個進程需要使用隨機數生成器來創建n/comm_sz的局部int型數據列表。每個進程先排序各自的局部列表，然后進程0收集並打印這些局部列表。然后，這些進程使用樹形結構通信合並全局列表給進程0，並打印最終結果。

　　

如圖，一般情況下，應由進程讀入數據，然后一層一層的傳到每個進程，每個進程排序完成之后，由上一層的進程進行合並。我這次做的只是輸入數n，由每個進程產生n/comm_sz個隨機數。

這里有一個問題，如何確定一個進程的父節點和子節點？

　　你最初可能認為讓在處理樹中的每個節點是一個單獨的進程中。這樣，你可以簡單地從二叉堆借用一個想法任何父節點的左子節點的下標2* K+1，右子節點的下標是2* K+2，一個節點的父節點是（K-1 ）/ 2，這在一個完全二叉樹中確定了父子關系。因此，一個內部節點將數據分成兩半，並發送給每個子進程進行處理。這樣葉節點，只是做了排序，內部節點等待然后從兩個子節點接收回數據，執行兩半的合並，以及（對於所有內部節點但非根節點本身）將結果發送到父節點。

　　但是這樣帶來一個問題，就是當葉節點進行排序時，父節點在等待子節點，沒有理由讓父節點空閑等待子節點進行排序，我們希望讓父節點充當左子節點進行工作，如圖：

　　這樣每個節點將數據分成兩半，每個節點自己處理數據的左半，右子節點處理數據的右半。這樣之后你就要確定每個節點的父節點和子節點，你可以根據樹形通信結構來確定每個節點的父節點和子節點。

　　我們要先確定每個節點在樹形通信結構的高度，葉節點的高度為0，這樣根節點0的高度是3。節點0需要與節點4進行通信，對於如何計算出節點0的右子節點是4，可以通過將1轉化成二進制並左移節點0所在的高度3減1位，即(myRank|(1<<2)),myRank表示節點的編號，這里是0，這樣就可以計算出高度為3的節點0的右子節點。下一步節點0需要與節點2通信，節點4需要與節點6通信，此時節點0和4的高度為2，因此需要計算他們的子節點時通過公式((myRank|(1<<1))，這樣計算出節點0的右子節點是2，節點4的右子節點時6，一般化的公式是(myRank|(myHeight-1))。

　　計算完每個節點的子節點后，還要計算每個節點的父節點，每個節點的父節點可以通過如下公式進行計算，myRank&~(1<<myHeight)。

　　下面是一個小的demo來展示通信過程，他展示了高度為3的通信樹節點之間的通信過程，

#include <stdio.h>

void communicate ( int myHeight, int myRank )
{  int parent = myRank & ~(1<<myHeight);

   if ( myHeight > 0 )
   {  int nxt     = myHeight - 1;
      int rtChild = myRank | ( 1 << nxt );

      printf ("%d sending data to %d\n", myRank, rtChild);
      communicate ( nxt, myRank );
      communicate ( nxt, rtChild );
      printf ("%d getting data from %d\n", myRank, rtChild);
   }
   if ( parent != myRank )
      printf ("%d transmitting to %d\n", myRank, parent);
}

int main ( void )
{  int myHeight = 3, myRank = 0;

   printf ("Building a height %d tree\n", myHeight);
   communicate(myHeight, myRank);
   return 0;
}

下面是輸出結果

Building a height 3 tree
0 sending data to 4
0 sending data to 2
0 sending data to 1
1 transmitting to 0
0 getting data from 1
2 sending data to 3
3 transmitting to 2
2 getting data from 3
2 transmitting to 0
0 getting data from 2
4 sending data to 6
4 sending data to 5
5 transmitting to 4
4 getting data from 5
6 sending data to 7
7 transmitting to 6
6 getting data from 7
6 transmitting to 4
4 getting data from 6
4 transmitting to 0
0 getting data from 4

了解了通信結構之后，下面是mpi_merge_sort.c文件，利用MPI寫的歸並排序，

#include<stdlib.h>
#include<stdio.h>
#include<time.h>
#include<math.h>
#include<string.h>
#include<mpi.h>
//讀取輸入的要排序的數的個數，並計算每個進程需要排序的個數
void read_n(int* n,int my_rank,MPI_Comm comm,int comm_sz,int* size);
//根據節點在通信樹中的最高度為數組分配空間
void allocate_arrays(int** left,int** right,int** merge,int size,int height,int* myHeight,int my_rank);
//在left數組中生成size個隨機數
void randnum(int* left,int size,int my_rank);
//以下三個函數時用於在通信樹的葉節點對生成在left中的size個隨機數進行快速排序
void swap(int *a,int *b);
int partition(int left[],int lo,int hi);
void QuickSort(int left[], int lo, int hi);
//父節點將接收子節點的數據與自己節點的數據進行歸並
void merge_sort(int* left,int* right,int* merge,int size);
//控制通信，確定從哪個子節點中接收數據或向哪個父節點發送數據
void communicate ( int Height, int my_rank,int size,int* right,int* left,int* merge,MPI_Comm comm);

int main(int argc,char* argv[]){
    int n,comm_sz,my_rank;
    int* left;
    int* right;
    int* merge;
    int size;
    int myHeight=0;
    MPI_Comm comm;
    int height=0;
    
    MPI_Init(NULL,NULL);
    comm=MPI_COMM_WORLD;
    MPI_Comm_size(comm,&comm_sz);
    MPI_Comm_rank(comm,&my_rank);
    read_n(&n,my_rank,comm,comm_sz,&size);
    int tt=comm_sz;
    for(int i=0;i<comm_sz;i++){
        tt=tt/2;
        if(tt==0)
            break;
        height++;    
    }
    allocate_arrays(&left,&right,&merge,size,height,&myHeight,my_rank);
    randnum(left,size,my_rank);
    QuickSort(left,0,size-1);
    
    communicate(height,my_rank,size,right,left,merge,comm);
    if(my_rank==0){
        if(height%2!=0)
            for(int i=0;i<n;i++)
                printf("%d    ",merge[i]);
        else
            for(int i=0;i<n;i++)
                printf("%d    ",left[i]);
        printf("\n");
    }
    MPI_Finalize();
    free(left);
    free(right);
    free(merge);
    
}
void read_n(int* n,int my_rank,MPI_Comm comm,int comm_sz,int* size){
    if(my_rank==0){
        printf("please intput the number of number\n");
        scanf("%d",n);    
    }
    MPI_Bcast (n,1,MPI_INT,0,comm);
    *size=*n/comm_sz;
}
void allocate_arrays(int** left,int** right,int** merge,int size,int height,int* myHeight,int my_rank){
    for(int i=0;i<height;i++){
        int parent=my_rank&~(1<<i);
        if(parent!=my_rank)
            break;
        (*myHeight)++;
    }
    *left=malloc((1<<(*myHeight))*size*sizeof(int));
    
    *right=malloc((1<<(*myHeight-1))*size*sizeof(int));
    *merge=malloc((1<<(*myHeight))*size*sizeof(int));
    
}
void randnum(int* left,int size,int my_rank){
    srand(my_rank);
    for(int i=0;i<size;i++)
        left[i]=rand();
}
void swap(int *a,int *b)
{
 int temp=*a;
 *a=*b;
 *b=temp;
}
int partition(int* left,int lo,int hi)
{
 int key=left[hi];  
 int i=lo-1;
 for(int j=lo;j<hi;j++)   
 {
  if(left[j]<=key)
  {
   i=i+1;
   swap(&left[i],&left[j]);
  }
 }
 swap(&left[i+1],&left[hi]);  
 return i+1;
}
void QuickSort(int *left, int lo, int hi)
{
    if (lo<hi)
    {
        int k = partition(left, lo, hi);
        QuickSort(left, lo, k-1);
        QuickSort(left, k+1, hi);
    }
}
void merge_sort(int* left,int* right,int* merge,int size){
    
    int lp=0,rp=0,mp=0;
    for(int i=0;i<2*size;i++){
        if(left[lp]<=right[rp]){
            merge[mp]=left[lp];
            lp++;
            mp++;    
        }else{
            merge[mp]=right[rp];
            rp++;
            mp++;        
        }
        if(lp==size||rp==size)
            break;
    }
    if(lp==size){
        memcpy(&merge[mp],&right[rp],(size-rp)*sizeof(int));
    }else{
        memcpy(&merge[mp],&left[lp],(size-lp)*sizeof(int));
    }    
}
void communicate ( int Height, int my_rank,int size,int* right,int* left,int* merge,MPI_Comm comm)
{      
    for(int i=0;i<=Height;i++){
         int parent=my_rank&~(1<<i);
        if(parent==my_rank){
            int rtChild=my_rank|(1<<i);
            
            MPI_Recv(right,size,MPI_INT,rtChild,0,comm,MPI_STATUS_IGNORE);
            
            merge_sort(left,right,merge,size);

            int* temp;
            temp=left;
            left=merge;
            merge=temp;
            size*=2;
            
        }else{
            MPI_Send(left,size,MPI_INT,parent,0,comm);
            break;
            
        }
    }
    
}

使用 mpicc -g -Wall -std=c99 -o mpi_merge_sort mpi_merge_sort.c 進行編譯

使用 mpirun -n 4 ./mpi_merge_sort 

注意輸入的進程個數一定是2的整數次冪，且輸入的要排序數的個數要能整除進程個數

下面是一次的運行結果：

please intput the number of number
12
190686788    483147985    844158168    846930886    846930886    1205554746    1505335290    1681692777    1681692777    1738766719    1804289383    1804289383

 參考資料：http://penguin.ewu.edu/~trolfe/ParallelMerge/ParallelMerge.html