並行計算——基於MPI實現埃拉托斯特尼篩法及性能優化

具體代碼可以在我碼雲里獲得：https://gitee.com/wu_ji666/mpi_Eratosthenes.git

一.步驟及操作：

1.1, VS2017上MPI的安裝和配置：

1.1.1安裝：

a.mpi官方下載地址：http://www.mpich.org/downloads/,下載后按照安裝指導操作即可

1.1.2配置：

配置具體步驟（每次新建mpi工程時都要重新對工程進行配置）：

右擊項目-->>屬性，進行配置： 

右上角-->>配置管理器-->>活動解決方案平台，選擇：x86（我是64位系統，如果是32位，選擇x64）； 

VC++目錄-->>包含目錄，添加：“D:\Program Files (x86)\Microsoft SDKs\MPI\Include;” 

VC++目錄-->>庫目錄，添加：“D:\Program Files (x86)\Microsoft SDKs\MPI\Lib\x64;”  

C/C++ -->> 預處理器-->>預處理器定義，添加：“MPICH_SKIP_MPICXX;” 

C/C++ -->> 代碼生成 -->> 運行庫，選擇：多線程調試（/MTd）； 

鏈接器 -->> 輸入 -->> 附加依賴項，添加：“msmpi.lib;”

1.1.3調試：

#include<stdio.h>

#include<mpi.h>

int main(int argc, char *argv[]) {

          int myid, numprocs;

    MPI_Init(&argc,&argv);

    MPI_Comm_rank(MPI_COMM_WORLD, &myid);

    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

 

    printf("%d Hello world from process %d \n",numprocs, myid);

 

    MPI_Finalize();

    return 0;

}

對這段代碼編譯並生成exe,用命令行執行mpiexec得到：

 

 圖1.1.1

表明成功！

 

 

1.2，   埃拉托斯特尼篩法基准代碼調試：

1.2.1.   關鍵代碼

就是基准代碼，做了一些修改：

a.加入#include <stdlib.h>

b.int改為long long int

c.加入結果寫入文件模塊（其他優化代碼里一直沿用）：

void write_txtfile(double k, int num_pro, long long int num)

{

    FILE * pFile = fopen("E:/VS_pro/文件名", "a+");

 

    if (NULL == pFile)

    {

         printf("error");

         return;

    }

    fprintf(pFile, "Thread num:%d,Arrange:%lld, Time:%10.6f\n", num_pro, num, k);

    fclose(pFile);

    return;

}

 

1.3. 埃拉托斯特尼篩法去偶數優化：

1.3.1. 思想和方法：

   利用“大於２的質數都是奇數”這一知識，首先去掉所有偶數，偶數必然不是素數，這樣相當於所需要篩選的數減少了一半，存儲和計算性能都得到提高。 

1.3.2.   關鍵代碼

 

int N = (n - 1) / 2;//just half of all

    low_value = (id * (N / p) + MIN(id, N % p)) * 2 + 3; //irst ele. in this pross.

    high_value = ((id + 1) * (N / p) + MIN(id + 1, N % p)) * 2 + 1;//last ele...

    size = (high_value - low_value) / 2 + 1;    //the size in this pross...

……

    // marked all

    for (int i = 0; i < size; i++) marked[i] = 0;

 

    if (!id) index = 0;

    prime = 3;

    do {

 

         if (prime * prime > low_value) {

             first = (prime * prime - low_value) / 2;

         }

         else {

 

             if (!(low_value % prime)) first = 0;

             //kicked out the even

             else if (low_value % prime % 2 == 0)

                  first = prime - ((low_value % prime) / 2);

             else first = (prime - (low_value % prime)) / 2;

         }

……

         if (!id) {

             while (marked[++index]);

             prime = index * 2 + 3; // start of 3

        }

 

事實上，去偶數重要的是索引的確定！

1.4.         埃拉托斯特尼篩法消除廣播：

1.4.1     思想和方法

基准代碼是通過進程0廣播下一個篩選倍數的素數。進程之間需要通過MPI_Bcast函數進行通信。通信就一定會有開銷，因此我們讓每個進程都各自找出它們的前sqrt(n)個數中的素數，在通過這些素數篩選剩下的素數，這樣一來進程之間就不需要每個循環廣播素數了，性能得到提高。

1.4.2關鍵代碼

/*

 廣播優化

    */

    //先找前sqrt(n)內的素數,再通過這些素數篩選后續素數

    int sqrt_N = (((int)sqrt((double)n)) - 1) / 2;

 

    low_NewArray_value = MIN(0, sqrt_N % p) * 2 + 3; //first ele. in this pross.（in NewArray）

    high_NewArray_value = ((sqrt_N / p) + MIN(1, sqrt_N % p)) * 2 + 1;//last ele...（in NewArray）

 

    // Bail out if all the primes used for sieving are not all held by process 0

    proc0_size = (sqrt_N - 1) / p;

 

    if ((2 + proc0_size) < (int)sqrt((double)sqrt_N)) {

         if (!id) printf("Too many processes \n");

         MPI_Finalize();

         exit(1);

    }

 

    NewMarked = (char *)malloc(sqrt_N);

    if (NewMarked == NULL) {

         printf("Cannot allocate enough memory \n");

         MPI_Finalize();

         exit(1);

    }

 

    //all ele. in NewArray was be marked

    for (int i = 0; i < sqrt_N; i++)

         NewMarked[i] = 0;

 

    index = 0;

 

    prime = 3;

    do {

         first = (prime * prime - low_NewArray_value) / 2;

 

         // 該素數的倍數排除

         for (int i = first; i < sqrt_N; i += prime) {

             NewMarked[i] = 1;

         }

 

         while (NewMarked[++index]);

         prime = index * 2 + 3;

 

    } while (prime * prime <= sqrt_N);

 

 

 

1.5.         埃拉托斯特尼篩法cache優化：

1.5.1.   思想和方法：

（0）Cache配置如圖：

 

圖1.5.1

（1）已知Cache是以CacheLine大小為單位交換數據的，考慮到CacheLine共用時的偽共享問題，首先進行Cacheline對齊，使每一個線程需要的數據在不同的Cacheline。我的電腦Cache配置已知CacheLine大小為64B,那么我們設置新的數據結構來完成數據Cacheline對齊：;並計算每個處理器占有的Cacheline數；

（2）考慮到Cache整體與內存的交換行為，我們進行分塊。已知三級緩存為6M(之所以用三級緩存是為了避免L2\L1級緩存太小導致Cache命中率過低，反而降低優化性能)。在已知每個處理器占有的CacheLine數基礎上，把這些CacheLIne集中成塊。

列如：每個處理器可以獲得的Cache大小為：Cache大小/核心數

占有CacheLIne數為：每個處理器可以獲得的Cache大小/數據結構占有空間大小

1.5.2.   關鍵代碼

代碼注釋詳細，思路可見一斑。

int Cache_linenum_pro = CACHE_SIZE / (CACHELINE_SIZE*p);//每個進程占有的cacheline

    int CacheBlock_size = Cache_linenum_pro * 8;//每個進程獲得的用於存取long long的塊大小

    int Block_N = CacheBlock_size - 1;

    int line_need = size / CacheBlock_size;//每個進程一共需要多少塊

    int line_rest = size % CacheBlock_size;//多出來的cacheline

    int time_UseCache = 0;

 

 

    // allocate this process 's share of the array

    marked = (char *)malloc(CacheBlock_size);

    if (marked == NULL) {

         printf("Cannot allocate enough memory \n");

         MPI_Finalize();

         exit(1);

    }

 

    //就是塊內去除偶數的判定過程

    count = 0;

 

    while (time_UseCache <= line_need) {

 

         //cache更新；

         Block_pos_last.value = (time_UseCache + 1) * Block_N + MIN(time_UseCache + 1, line_rest) - 1 + (id * (N / p) + MIN(id, N % p));

         Block_pos_first.value = time_UseCache * Block_N + MIN(time_UseCache, line_rest) + (id * (N / p) + MIN(id, N % p));

         Block_low_value.value = Block_pos_first.value * 2 + 3;

         if (time_UseCache == line_need) {

             Block_high_value.value = high_value;

             Block_pos_last.value = (id + 1) * (N / p) + MIN(id + 1, N % p) - 1;

             CacheBlock_size = ((Block_high_value.value - Block_low_value.value) >> 1) + 1;

         }

         else {

             Block_high_value.value = (Block_pos_last.value + 1) * 2 + 1;

         }

 

 

         index = 0;

         prime.value = 3;

         count_cacheBlock = 0;

 

         for (int i = 0; i < CacheBlock_size; i++) marked[i] = 0;

 

         // 塊內素數

         do {

             if (prime.value * prime.value > Block_low_value.value) {

                  first.value = (prime.value * prime.value - Block_low_value.value) >> 1;

             }

             else {

                  if (!(Block_low_value.value % prime.value)) first.value = 0;

                  else if (Block_low_value.value % prime.value % 2 == 0) first.value = prime.value - ((Block_low_value.value % prime.value) >> 1);

                  else first.value = (prime.value - (Block_low_value.value % prime.value)) >> 1;

             }

             for (int i = first.value; i < CacheBlock_size; i += prime.value) {

                  marked[i] = 1;

             }

             while (NewMarked[++index]);

 

 

             prime.value = index * 2 + 3;

         } while (prime.value * prime.value <= Block_high_value.value);

 

 

         // 塊內計數

         for (int i = 0; i < CacheBlock_size; i++) {

             if (marked[i] == 0) {

                  count_cacheBlock++;

             }

         }

 

         //all blockes togather

         count += count_cacheBlock;

         // 處理下一個塊

         time_UseCache++;

    }

 

 

二、實驗數據及結果分析：

2.1 實驗結果

2.1.1.正確性展示（部分）：

 

a.1e8:共5761455個素數；

1e8結果

b, 1e9:共50847534個素數；

 1e9結果

 

c, 1e10:共105097565個素數。

1e10結果 

2.1.2．結果評價：

a.1e8

執行時間：

圖2.1.1

圖2.1.1中可以看到在1-16個核心數范圍內，所有程序的執行時間都下降明顯，當核心數超過16后，執行時間逐漸平緩。

加速比圖：

 

圖2.1.2

圖2.1.2（加速比圖中）易於發現，去偶數與原基准程序成線性，而去廣播和Cache優化則在核心數16左右達到做大，之后趨於平緩。

並行效率圖：

 圖2.1.3

圖中顯示了所有加速程序並行效率都隨着核心數增加而下降；並且去偶數優化和基准程序基本一致，Cache優化並行效率同等核心數處於最高，去廣播優化前期與基准程序一致，隨着核心數增大，其並行效率明顯高於基准程序。

縱比加速圖：

同核心數下各優化方式相對於基准代碼的加速效果：

 圖2.1.4

可以看出去偶數相較於基准程序加速比接近2倍，去廣播可接近4倍，Cache優化效果最好，在核心數16時達到最大加速倍數為10-11倍；其中去偶數優化優化效果可以看做常數2，去廣播和Cache優化則成峰分布，在16核心數左右達到最大。

  

b. 在此將數據擴大到1e10：

執行時間：

圖2.1.5、

結果特征與1e8時一致。

加速比圖：

 

圖2.1.6

有趣的是，在1e8時加速比最大出現在16核心數，在1e10則為8；同時線性的特征消失（無論是基准程序還是去偶數）。

阿姆達爾效應驗證：

 

 

 圖2.1.7

這個就更有趣了，按照阿姆達爾效應，隨着問題規模增加，同一並行算法的加速比應該上升，然而結果是1e8的加速比遙遙領先，而1e9和1e10則比較接近。

2.1.3. 結果分析：

（1）.去偶數優化其實只是數量上線性的減少了問題規模（接近一半），通信和廣播模式沒有改變，所有優化效果是線性，接近2倍，並行效率也相似；

(2).核心數少時廣播代價並不大，隨着核心數增加廣播代價增大，所以在核心數很大時去廣播優化就明顯優於去偶數。

(3).顯然存在計算代價、通信代價（或者廣播代價、計算代價、通信代價），而通信（通信和廣播）會隨着核心數增加增大，所以Cache優化和去廣播優化在加速計算時，如果通信代價過大，計算加速帶來的收益就下降，加速比曲線存在峰；而通信代價與直接計算代價的相對比例在不同問題規模不同，使不同問題規模下最大加速比對應核心數不同；

(4).通信代價存在，導致並行效率隨核心數增大而下降；

(5).問題規模溢出本人計算機性能，出現大規模加速比下降，反阿姆達爾效應。