【並行計算】用MPI進行分佈式內存編程（一）

    經過上一篇關於並行計算準備部分的介紹，咱們知道MPI（Message-Passing-Interface 消息傳遞接口）實現並行是進程級別的，經過通訊在進程之間進行消息傳遞。MPI並非一種新的開發語言，它是一個定義了能夠被C、C++和Fortran程序調用的函數庫。這些函數庫裏面主要涉及的是兩個進程之間通訊的函數。MPI能夠在Windows和linux環境中都有相應的庫，本篇以Windows10做爲演示開發環境。

一、Windows10+VS 2015上搭建MPI開發環境

• 下載 mpi for windows

Windows爲了兼容MPI，本身作了一套基於通常我的電腦的MPI實現。若是要安裝正真意義上的MPI的話，請直接去www.mpich.org下載，裏面根據對應的系統下載相應的版本。

 

我因爲是須要在我的筆記本上進行實驗，就用的是微軟的 HPC Pack 2008 R2 MS-MPI Redistributable Package with Service Pack 4 - 中文(簡體)，下載地址： http://www.microsoft.com/zh-cn/download/details.aspx?id=14737 。

• 安裝 mpi

   個人電腦是64位的，因此安裝的是 mpi_x64.msi ，默認安在C:\Program Files\Microsoft HPC Pack 2008 R2，在此，爲了以後調試代碼方便，最好設置一下環境變量：在用戶變量PATH中，加入：C:\Program Files\Microsoft HPC Pack 2008 R2\Bin\。

• 配置2015

配置目錄，即加載Include和Lib庫

 

加載依賴項

• 編譯

幾乎全部人的第一個程序是從「hello，world」程序開始學習的，我也寫了一個這樣測試例子：

#include "mpi.h"  
#include <stdio.h>  

int main(int argc, char* argv[])
{
    int rank, numproces;
    int namelen;
    char processor_name[MPI_MAX_PROCESSOR_NAME];

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);//得到進程號
    MPI_Comm_size(MPI_COMM_WORLD, &numproces);//返回通訊子的進程數

    MPI_Get_processor_name(processor_name, &namelen);
    fprintf(stderr, "hello world! process %d of %d on %s\n", rank, numproces, processor_name);
    MPI_Finalize();

    return 0;
}

上述代碼中，第1行中的#include "mpi.h" 頭文件必須包含，在VS2015下編譯生成exe文件（生成在debug文件中），經過cmd命令，進入debug文件夾目錄中，敲入：mpiexec –n 4 TestForMPI.exe。其中命令中-n 4 表示使用4個進程進行並行計算，具體結果如圖所示：

 

二、理論知識

    經過上述的例子，咱們對MPI編寫的並行計算有了一個初步的認識，可是咱們還不知道如何去真正編寫一個MPI的並行程序，這須要咱們學習必定的理論知識。在上面的例子中，有幾個函數對於初學MPI的人來講並不明白是什麼意思，下面就從這些函數入手。

MPI_Init：告知MPI系統進行全部必要的初始化設置。它是寫在啓動MPI並行計算的最前面的。具體的語法結構爲：

MPI_Init(
     int* argc_p,
     char*** argv_p
    );

    參數argc_p和argv_p分別指向main函數中的指針參數，爲了弄明白這部分，還得從main函數的參數提及：Ｃ語言規定main函數的參數只能有兩個，習慣上這兩個參數寫爲argc和argv。所以，main函數的函數頭可寫爲： main (argc,argv)。Ｃ語言還規定argc(第一個形參)必須是整型變量,argv( 第二個形參)必須是指向字符串的指針數組。其中argc參數表示了命令行中參數的個數(注意：文件名自己也 算一個參數)，argc的值是在輸入命令行時由系統按實際參數的個數自動賦予的。例若有命令行爲： C:">E6 24 BASIC dbase FORTRAN因爲文件名E6 24自己也算一個參數，因此共有4個參數，所以argc取得的值爲4。argv參數是字符串指針數組，其各元素值爲命令行中各字符串(參數均按字符串處 理)的首地址。

然而在MPI_Init函數中，並不必定都須要設置argc_p和argv_p這兩個參數的，不須要的時候，將它們設置爲NULL便可。

通信子（communicator）：MPI_COMM_WORLD表示一組能夠互相發送消息的進程集合。

MPI_Comm_rank:用來獲取正在調用進程的通訊子中的進程號的函數。

MPI_Comm_size:用來獲得通訊子的進程數的函數。

這兩個函數的具體結構以下：

int MPIAPI MPI_Comm_rank(
    __in MPI_Comm comm,
    __out int* rank
    );

int MPIAPI MPI_Comm_size(
    __in MPI_Comm comm,
    __out int* size
    );

它們的第一個參數都傳入通訊子做爲參數，第二參數都用傳出參數分別把正在調用通訊子的進程號和通訊的個數。

MPI_Finalize：告知MPI系統MPI已經使用完畢。它老是放到作並行計算的功能塊的最後面，在此函數以後就不能再出現任何有關MPI相關的東西了。

以上只是表達了做爲一個MPI並行計算的基本結構，並無真正涉及進程之間的通訊，爲了更好的進行並行，必然須要進程間的通訊，下面介紹兩個進程間通訊的函數，它們就是MPI_Send和MPI_Recv，分別用於消息的發送和接收。

MPI_Send：阻塞型消息發送。其結構爲：

int MPI_Send (void *buf, int count, MPI_Datatype datatype,int dest, int tag,MPI_Comm comm)

參數buf爲發送緩衝區；count爲發送的數據個數；datatype爲發送的數據類型；dest爲消息的目的地址(進程號)，其取值範圍爲0到np－1間的整數(np表明通訊器comm中的進程數) 或MPI_PROC_NULL；tag爲消息標籤，其取值範圍爲0到MPI_TAG_UB間的整數；comm爲通訊器。

MPI_Recv：阻塞型消息接收。

int MPI_Recv (void *buf, int count, MPI_Datatype datatype,int source, int tag, MPI_Comm comm,MPI_Status *status)

參數buf爲接收緩衝區；count爲數據個數，它是接收數據長度的上限，具體接收到的數據長度可經過調用MPI_Get_count函數獲得；datatype爲接收的數據類型；source爲消息源地址(進程號)，其取值範圍爲0到np－1 間的整數(np表明通訊器comm 中的進程數)，或MPI_ANY_SOURCE，或MPI_PROC_NULL；tag爲消息標籤，其取值範圍爲0到MPI_TAG_UB間的整數或MPI_ANY_TAG；comm爲通訊器；status返回接收狀態。

MPI_Status：返回消息傳遞的完成狀況。數據結構的相關變量的意義就比較多了，具體能夠參數使用手冊。

typedef struct {
... ...
int MPI_SOURCE;             /*消息源地址*/
int MPI_TAG;                /*消息標籤*/
int MPI_ERROR;              /*錯誤碼*/
... ...
} MPI_Status;

三、舉例

介紹了最基本的進程間通訊的函數，咱們就能編寫一個更爲複雜也更爲有意義的程序，經過編寫程序來實現數據積分中的梯形積分法。

梯形積分法的基本思想是：將x軸上的區間劃分爲n個等長的子區間。而後計算子區間的和。

假設子區間的端點爲x_i和x_i+1，那麼子區間的長度爲：h=x_i+1-x_i。那麼梯形的面積就爲：

因爲n個子區間是等分的，邊界分別爲xi=a和x=b，則：

這片區域的全部梯形的面積和爲：

變換爲：

所以，串行的程序代碼就能夠這樣寫：

h = (b - a) / h;
approx = (f(a) + f(b)) / 2;
for (int i = 1; i < n-1; i++)
{
    x_i = a + i*h;
    approx += f(x_i);
}
approx = h*approx;

經過對串行程序的分析，對於這個例子，咱們能夠識別出兩種任務：第一種獲取單個矩形區域的面積，另外一種是計算這些區域的面積和。

假設求f(x)=x³將梯形劃分爲1024個子區域計算[0,3]區域內的積分。

#include "mpi.h"  
#include <stdio.h>  
#include <cmath>

double Trap(double left_endpt, double right_endpt, double trap_count, double base_len);
double f(double x);

int main(int argc, char* argv[])
{
    int my_rank = 0, comm_sz = 0, n = 1024, local_n = 0;
    double a = 0.0, b = 3.0, h = 0, local_a = 0, local_b = 0;
    double local_int = 0, total_int = 0;
    int source;

    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
    MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

    h = (b - a) / n;       /*  h is the same for all processes  */
    local_n = n / comm_sz; /*  So is the number of trapezoids */

    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;
}
//子區域的積分函數
double Trap(double left_endpt, double right_endpt, double trap_count, double base_len)
{
    double estimate = 0, x = 0;
    int i;

    estimate = (f(left_endpt) + f(right_endpt)) / 2.0;
    for (i = 1; i <= trap_count - 1; i++)
    {
        x = left_endpt + base_len;
        estimate += f(x);
    }
    estimate = estimate*base_len;
    return estimate;
}
//數學函數
double f(double x)
{
    return pow(x, 3);
}

上述代碼中，運行結果

 

這段程序代碼中的意思是，經過輸入的進程數，將1024個劃分的子區間等分的分配到控制檯輸入的進程（100個）進行子任務求解，求解完成以後，1-99進程計算的結果值經過MPI_Send函數發送出去，而0號進程使用MPI_Recv函數接收匯總，將每一個進程的結果求和，獲得區間的積分值。

本次並行計算消息之間的通訊以下圖：

至此，咱們已經能使用MPI_Send消息發送函數和MPI_Recv消息接收函數進行簡單的並行程序計算了，但咱們想一下，最後的求和都是用0號進程去作，爲了更加的提升性能，還需更進一步使用集合通訊，下一章節將深的進行講解。