關於mpi的理論知識以及編寫程序來實現數據積分中的梯形積分法。

 幾乎所有人的第一個程序是從“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)=x3將梯形划分為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消息接收函數進行簡單的並行程序計算了,但最后的求和都是用號進程去做。