【MPI學習4】MPI並行程式設計模式：非阻塞通訊MPI程式設計
這一章講了MPI非阻塞通訊的原理和一些函式介面，最後再用非阻塞通訊方式實現Jacobi迭代，記錄學習中的一些知識。
（1）阻塞通訊與非阻塞通訊
阻塞通訊呼叫時，整個程式只能執行通訊相關的內容，而無法執行計算相關的內容；
非阻塞呼叫的初衷是儘量讓通訊和計算重疊進行，提高程式整體執行效率。
整體對比見下圖：
（2）非阻塞通訊的要素
非阻塞通訊呼叫返回意味著通訊開始啟動；而非阻塞通訊完成則需要呼叫其他的介面來查詢。
要素1：非阻塞通訊的呼叫介面
要素2：非阻塞通訊的完成查詢介面
理想的非阻塞通訊設計應該如下：
非阻塞通訊的 傳送 和 接受 過程都需要同時具備以上兩個要素，“呼叫+完成”
“呼叫”按照通訊方式的不同（標準、快取、同步、就緒），有各種函式介面，具體用到哪個就查手冊的性質。
這裡“完成”是重點，因為程式設計師需要知道非阻塞呼叫是否執行完成了，來做下一步的操作。
MPI為“完成”定義了一個內部變數MPI_Request request，每個request與一個在非阻塞呼叫發生時與該呼叫發生關聯（這裡的呼叫包括髮送和接收）。
“完成”不區分通訊方式的不同，統一用MPI_Wait系列函式來完成，這裡對MPI_Wait函式做一點說明：
1）MPI_Wait(MPI_Request *request)，均等著request執行完畢了，再往下進行
2）對於非重複非阻塞通訊，MPI_Wait系列函式呼叫的返回，還意味著request物件被釋放了，程式設計師不用再顯式釋放request變數。
3）對於重複非阻塞通訊，MPI_Wait系列函式呼叫的返回，意味著將於request物件關聯的非阻塞通訊處於不啟用狀態，並不釋放request
關於2）3）看後面的程式碼示例就瞭解了
（3）非阻塞呼叫實現Jacobi迭代
有了非阻塞呼叫的技術，可以再將Jacobi迭代的程式效率提升，其總體的實現思路如下：
1）先計算Jacobi迭代下次計算所需要的邊界資料，這些資料與每個計算節點中的計算無關，可以先獨立計算好
2）啟動非阻塞通訊，將邊界資料在程序間傳遞
3）計算每個計算節點可以獨立計算的部分；此時，2）中啟動的非阻塞通訊也在進行中，這時通訊和計算就重疊了
4）等著非阻塞通訊完成，再進行下一次迭代
再回顧一下之前用阻塞通訊實現Jacobi迭代的思路：
1）先傳遞邊界資料
2）等著資料都傳遞完了，再進行計算
3）等著計算完成了，進行下一次迭代
可以看到阻塞通訊中實現Jacobi迭代的程式中，在同一計算節點下，通訊和計算是分別進行的，效率不如非阻塞通訊。
總結起來：
“單機程式 → 阻塞通訊MPI程式” 實現單機計算到多機計算，用並行代替序列提高效率。
“阻塞MPI程式 → 非阻塞MPI程式” 不僅將多臺機器之間的並行，而且還能將每臺機器的通訊與計算過程並行，實現更高效的並行。
（4）非阻塞通訊實現Jacobi迭代的程式碼
書上的原始碼是Fortan的，資料儲存是列優先的，矩陣按列分塊；下面的程式碼是我翻譯的C的程式碼，資料儲存是行優先的，矩陣按行分塊。
1 #include "mpi.h" 2 #include <stdio.h> 3 #include <stdlib.h> 4 5 #define N 8 6 #define SIZE N/4 7 #define T 2 8 9 void print_matrix(int myid, float myRows[][N]); 10 11 int main(int argc, char *argv[]) 12 { 13 float matrix1[SIZE+2][N], matrix2[SIZE+2][N]; 14 int myid; 15 MPI_Status status[4]; 16 MPI_Request request[4]; 17 18 MPI_Init(&argc, &argv); 19 MPI_Comm_rank(MPI_COMM_WORLD, &myid); 20 21 // 初始化 22 int i,j; 23 for(i=0; i<SIZE+2; i++) 24 { 25 for(j=0; j<N; j++) 26 { 27 matrix1[i][j] = matrix2[i][j] = 0; 28 } 29 } 30 if(0==myid) // 按行劃分 上面第一分塊矩陣 上邊界 31 { 32 for(j=0; j<N; j++) matrix1[1][j] = matrix2[1][j] = N; 33 } 34 if (3==myid) { // 按行劃分 最下面一分塊矩陣 下邊界 35 for(j=0; j<N; j++) matrix1[SIZE][j] = matrix2[SIZE][j] = N; 36 } 37 for(i=1; i<SIZE+1; i++) // 每個矩陣的兩側邊界 38 { 39 matrix1[i][0] = matrix1[i][N-1] = matrix2[i][0] = matrix2[i][N-1] = N; 40 } 41 // 引入虛擬程序 並計算每個程序上下相鄰程序 42 int up_proc_id = myid==0 ? MPI_PROC_NULL : myid-1; 43 int down_proc_id = myid==3 ? MPI_PROC_NULL : myid+1; 44 // jacobi迭代過程 45 int t,row,col; 46 for(t=0; t<T; t++) 47 { 48 // 1 計算邊界資料 49 if(0==myid) // 最上的矩陣塊 50 { 51 for (col=1; col<N-1; col++) 52 { 53 matrix2[SIZE][col] = (matrix1[SIZE][col-1]+matrix1[SIZE][col+1]+matrix1[SIZE+1][col]+matrix1[SIZE-1][col])*0.25; 54 } 55 } 56 else if (3==myid) { // 最下的矩陣塊 57 for (col=1; col<N-1; col++) 58 { 59 matrix2[1][col] = (matrix1[1][col-1]+matrix1[1][col+1]+matrix1[2][col]+matrix1[0][col])*0.25; 60 } 61 } 62 else { 63 for(col=1; col<N-1; col++) // 中間的矩陣塊 64 { 65 matrix2[SIZE][col] = (matrix1[SIZE][col-1]+matrix1[SIZE][col+1]+matrix1[SIZE+1][col]+matrix1[SIZE-1][col])*0.25; 66 matrix2[1][col] = (matrix1[1][col-1]+matrix1[1][col+1]+matrix1[2][col]+matrix1[0][col])*0.25; 67 } 68 } 69 // 2 利用非阻塞函式傳遞邊界資料 為下一次計算做準備 70 int tag1 = 1, tag2 = 2; 71 MPI_Isend(&matrix2[1][0], N, MPI_FLOAT, up_proc_id, tag1, MPI_COMM_WORLD, &request[0]); 72 MPI_Isend(&matrix2[SIZE][0], N, MPI_FLOAT, down_proc_id, tag2, MPI_COMM_WORLD, &request[1]); 73 MPI_Irecv(&matrix1[SIZE+1][0], N, MPI_FLOAT, down_proc_id, tag1, MPI_COMM_WORLD, &request[2]); 74 MPI_Irecv(&matrix1[0][0], N, MPI_FLOAT, up_proc_id, tag2, MPI_COMM_WORLD, &request[3]); 75 // 3 計算中間資料 76 int begin_row = 0==myid ? 2 : 1; 77 int end_row = 3==myid ? (SIZE-1) : SIZE; 78 for (row=begin_row; row<end_row; row++) 79 { 80 for (col=1; col<N-1; col++) 81 { 82 matrix2[row][col] = (matrix1[row][col-1]+matrix1[row][col+1]+matrix1[row+1][col]+matrix1[row-1][col])*0.25; 83 } 84 } 85 // 4 更新矩陣 並等待各個程序間資料傳遞完畢 86 for (row=begin_row; row<=end_row; row++) 87 { 88 for (col=1; col<N-1; col++) 89 { 90 matrix1[row][col] = matrix2[row][col]; 91 } 92 } 93 MPI_Waitall(4, &request[0], &status[0]); 94 } 95 MPI_Barrier(MPI_COMM_WORLD); 96 print_matrix(myid, matrix1); 97 MPI_Finalize(); 98 } 99 100 101 void print_matrix(int myid, float myRows[][N]) 102 { 103 int i,j; 104 int buf[1]; 105 MPI_Status status; 106 buf[0] = 1; 107 if ( myid>0 ) { 108 MPI_Recv(buf, 1, MPI_INT, myid-1, 0, MPI_COMM_WORLD, &status); 109 } 110 printf("Result in process %d:\n", myid); 111 for ( i = 0; i<SIZE+2; i++) 112 { 113 for ( j = 0; j<N; j++) 114 printf("%1.3f\t", myRows[i][j]); 115 printf("\n"); 116 } 117 if ( myid<3 ) { 118 MPI_Send(buf, 1, MPI_INT, myid+1, 0, MPI_COMM_WORLD); 119 } 120 MPI_Barrier(MPI_COMM_WORLD); 121 }
程式的執行結果如下：
上述程式設計的邏輯如下：
1）各個分塊矩陣的邊界資料是可以需要通訊交換的
2）先計算邊界資料，儘量把需要通訊交換而且又相對獨立的資料先計算出來
3）用非阻塞通訊傳遞分塊矩陣的邊界資料；同時每個節點內計算內部的資料；計算與通訊並行
4）等到每個計算節點的2個傳送、2個接收，總共4個非阻塞呼叫都完成了，進行下一輪迭代
（5）重複非阻塞通訊
上面實現Jacobi迭代的程式碼中，以程序1和程序2為例：
1）迭代一輪二者之間就需要互相通訊一次
2）每次互相通訊，隨著MPI_Wait的執行，request通訊物件釋放，兩個程序通訊完全被切斷了
3）兩個程序之間每次通訊，有一些通訊連線操作都是重複的，最好不用每次通訊都重新執行這些連線操作，以此提高效率
4）因此，比上面實現jacobi迭代更優化一些的做法是：每次不完全掐斷兩個程序的非阻塞通訊，保持那些基礎的通用的操作，每次迭代只需要更新需要傳輸的資料，再啟用兩個程序之間的非阻塞通訊
依照上面的思路，MPI給出了重複非阻塞的通訊呼叫實現。用重複非阻塞的通訊再實現一次Jacobi迭代，程式碼如下：
1 #include "mpi.h" 2 #include <stdio.h> 3 #include <stdlib.h> 4 5 #define N 8 6 #define SIZE N/4 7 #define T 2 8 9 void print_matrix(int myid, float myRows[][N]); 10 11 int main(int argc, char *argv[]) 12 { 13 float matrix1[SIZE+2][N], matrix2[SIZE+2][N]; 14 int myid; 15 MPI_Status status[4]; 16 MPI_Request request[4]; 17 18 MPI_Init(&argc, &argv); 19 MPI_Comm_rank(MPI_COMM_WORLD, &myid); 20 21 // 初始化 22 int i,j; 23 for(i=0; i<SIZE+2; i++) 24 { 25 for(j=0; j<N; j++) 26 { 27 matrix1[i][j] = matrix2[i][j] = 0; 28 } 29 } 30 if(0==myid) // 按行劃分 上面第一分塊矩陣 上邊界 31 { 32 for(j=0; j<N; j++) matrix1[1][j] = matrix2[1][j] = N; 33 } 34 if (3==myid) { // 按行劃分 最下面一分塊矩陣 下邊界 35 for(j=0; j<N; j++) matrix1[SIZE][j] = matrix2[SIZE][j] = N; 36 } 37 for(i=1; i<SIZE+1; i++) // 每個矩陣的兩側邊界 38 { 39 matrix1[i][0] = matrix1[i][N-1] = matrix2[i][0] = matrix2[i][N-1] = N; 40 } 41 // 引入虛擬程序 並計算每個程序上下相鄰程序 42 int up_proc_id = myid==0 ? MPI_PROC_NULL : myid-1; 43 int down_proc_id = myid==3 ? MPI_PROC_NULL : myid+1; 44 // 初始化重複非阻塞通訊 45 int tag1 = 1, tag2 = 2; 46 MPI_Send_init(&matrix2[1][0], N, MPI_FLOAT, up_proc_id, tag1, MPI_COMM_WORLD, &request[0]); 47 MPI_Send_init(&matrix2[SIZE][0], N, MPI_FLOAT, down_proc_id, tag2, MPI_COMM_WORLD, &request[1]); 48 MPI_Recv_init(&matrix1[SIZE+1][0], N, MPI_FLOAT, down_proc_id, tag1, MPI_COMM_WORLD, &request[2]); 49 MPI_Recv_init(&matrix1[0][0], N, MPI_FLOAT, up_proc_id, tag2, MPI_COMM_WORLD, &request[3]); 50 // jacobi迭代過程 51 int t,row,col; 52 for(t=0; t<T; t++) 53 { 54 // 1 計算邊界資料 55 if(0==myid) // 最上的矩陣塊 56 { 57 for (col=1; col<N-1; col++) 58 { 59 matrix2[SIZE][col] = (matrix1[SIZE][col-1]+matrix1[SIZE][col+1]+matrix1[SIZE+1][col]+matrix1[SIZE-1][col])*0.25; 60 } 61 } 62 else if (3==myid) { // 最下的矩陣塊 63 for (col=1; col<N-1; col++) 64 { 65 matrix2[1][col] = (matrix1[1][col-1]+matrix1[1][col+1]+matrix1[2][col]+matrix1[0][col])*0.25; 66 } 67 } 68 else { 69 for(col=1; col<N-1; col++) // 中間的矩陣塊 70 { 71 matrix2[SIZE][col] = (matrix1[SIZE][col-1]+matrix1[SIZE][col+1]+matrix1[SIZE+1][col]+matrix1[SIZE-1][col])*0.25; 72 matrix2[1][col] = (matrix1[1][col-1]+matrix1[1][col+1]+matrix1[2][col]+matrix1[0][col])*0.25; 73 } 74 } 75 // 2 啟動重複非阻塞通訊 76 MPI_Startall(4, &request[0]); 77 // 3 計算中間資料 78 int begin_row = 0==myid ? 2 : 1; 79 int end_row = 3==myid ? (SIZE-1) : SIZE; 80 for (row=begin_row; row<end_row; row++) 81 { 82 for (col=1; col<N-1; col++) 83 { 84 matrix2[row][col] = (matrix1[row][col-1]+matrix1[row][col+1]+matrix1[row+1][col]+matrix1[row-1][col])*0.25; 85 } 86 } 87 // 4 更新矩陣 並等待各個程序間資料傳遞完畢 88 for (row=begin_row; row<=end_row; row++) 89 { 90 for (col=1; col<N-1; col++) 91 { 92 matrix1[row][col] = matrix2[row][col]; 93 } 94 } 95 MPI_Waitall(4, &request[0], &status[0]); 96 } 97 int n; 98 for(n = 0; n < 4; n++) MPI_Request_free(&request[n]); // 釋放非阻塞通訊物件 99 MPI_Barrier(MPI_COMM_WORLD); 100 print_matrix(myid, matrix1); 101 MPI_Finalize(); 102 } 103 104 105 void print_matrix(int myid, float myRows[][N]) 106 { 107 int i,j; 108 int buf[1