Cannon乘法的MPI實現(xiàn)

1、并行算法實踐要求學生在KD60實驗平臺上設(shè)計并行算法并實現(xiàn)。實驗平臺由一塊處理板、一塊監(jiān)控板和一塊背板等組成。處理板邏輯結(jié)構(gòu)如圖1所示。處理板承載4個處理單元，每個處理單元包括一個龍芯3號四核CPU、2GB DDR2內(nèi)存、RTL8110千兆以太網(wǎng)卡芯片、BIOS Flash、串口收發(fā)芯片以及電源變換電路等。四個龍芯3號處理器通過HT總線實現(xiàn)互連。監(jiān)控電路檢測4個處理單元的狀態(tài)，并實現(xiàn)對其控制。 圖1 處理板邏輯結(jié)構(gòu)實驗平臺的系統(tǒng)軟件以開源軟件為主（見圖2），具有兼容性強、易維護、易升級、易使用等特點。處理單元操作系統(tǒng)為Debian GNU/Linux無盤系統(tǒng)，采用穩(wěn)定高效的內(nèi)核。圖2 軟件系統(tǒng)

2、結(jié)構(gòu)要求選修并行算法實踐同學在下面表1中選一個題目，（1）闡述基本原理（包括對算法改進和優(yōu)化方法）；(2)根據(jù)KD60實驗平臺設(shè)計實驗方法和步驟（包括主要調(diào)試過程要求拷屏）。(3) 數(shù)據(jù)及結(jié)果分析：根據(jù)不同的實驗內(nèi)容，記錄具體的實驗數(shù)據(jù)或程序運行結(jié)果（要求拷屏）。實驗數(shù)據(jù)量較大時，最好制成表格形式。附上程序功能、模塊說明、完整源代碼，源代碼中盡量多帶注釋；(4)分析和總結(jié)：對程序與算法性能改進結(jié)論，總結(jié)和體會。表1 并行算法實踐題目 序號題目名稱基本方法和內(nèi)容要求1LU分解的Open MP實現(xiàn)編寫LU分解的Open MP程序2KMP算法的Open MP實現(xiàn)編寫KMP算法的OpenMP程序3高斯

3、消元法解線性方程組的Open MP實現(xiàn)編寫高斯消元法解線性方程組的OpenMP程序4高斯消元法解線性方程組的MPI實現(xiàn)編寫高斯消元法解線性方程組的MPI程序5高斯-塞德爾迭代解線性方程組的MPI實現(xiàn)編寫高斯-塞德爾迭代解線性方程組的MPI程序6Cannon乘法的MPI實現(xiàn)編寫Cannon乘法的MPI程序7LU分解的MPI實現(xiàn)編寫LU分解的MPI程序8隨機串匹配算法的MPI實現(xiàn)編寫隨機串匹配算法的MPI程序9單源最短路徑Dijkstra 算法的MPI實現(xiàn)編寫單源最短路徑Dijkstra算法的MPI程序10快速排序算法的MPI實現(xiàn)編寫快速排序算法的MPI程序11KMP串匹配的MPI實現(xiàn)編寫KMP串

4、匹配算法的MPI程序圖2 軟件系統(tǒng)結(jié)構(gòu)Cannon乘法的MPI實現(xiàn)及性能分析摘要：cannon算法是矩陣的并行乘法，屬于數(shù)值并行算法MPI編程實現(xiàn)一篇，其中關(guān)于數(shù)值并行算法MPI編程由于要處理的數(shù)據(jù)量巨大，程序循環(huán)次數(shù)多，對于串行而言，處理時間將非常長，將其并行化非常必要。本文將矩陣數(shù)據(jù)進行棋盤劃分成多個子矩陣，再分別指派給多個處理器，使個處理器并行運算。關(guān)鍵字：cannon乘法 并行計算 數(shù)據(jù)劃分一、Cannon乘法的MPI實現(xiàn)基本原理Cannon乘法屬于數(shù)值并行算法MPI編程實現(xiàn)一篇，其中關(guān)于數(shù)值并行算法MPI編程由于要處理的數(shù)據(jù)量巨大，程序循環(huán)次數(shù)多，對于串行而言，處理時間將非常長，使其

5、并行化的一般方法有：1）數(shù)據(jù)相關(guān)分析 2）數(shù)據(jù)劃分和處理器指派 3）循環(huán)重構(gòu)對原有程序并行化，首先要分析計算程序中所有語句間的依賴關(guān)系，這稱之為相關(guān)分析。本項目Cannon乘法的mpi實現(xiàn)，是矩陣運算，階往往都很高，而且行列之間數(shù)據(jù)依賴關(guān)系也不強，所以就對矩陣進行劃分，然后指派給不同的處理器進行處理。最常用的矩陣劃分有帶狀劃分和塊狀劃分。1 帶狀劃分方法帶狀劃分又叫行列劃分，就是將矩陣整行或整列地分成若干組，各組指派給一個處理器。也可以將若干行或列指派給一個處理器，而且這些行和列可以是連續(xù)的，也可以是等間距的，前者稱為塊帶狀的，后者稱為循環(huán)帶狀的。2. 塊狀劃分方法塊狀劃分又叫棋盤劃分，就是將

6、矩陣劃分成若干個子矩陣，每個子矩陣指派給一個處理器，此時任意處理器均不含整行或整列。和帶狀劃分類似，棋盤劃分也可分為塊棋盤劃分和循環(huán)棋盤劃分。棋盤劃分比帶狀劃分可開發(fā)更高的并行度，Cannon乘法的mpi實現(xiàn)也正是基于棋盤劃分的并行實現(xiàn)。循環(huán)重構(gòu)是指在數(shù)據(jù)分解之后，相應(yīng)地將串行程序循環(huán)部分進行重構(gòu)，以實現(xiàn)這種劃分所確定的并行計算，主要方法有1）循環(huán)交換 2）拉伸法 3）分裂法 4）輪轉(zhuǎn)法 5）并列法 在三種程序并行化的方法中，數(shù)據(jù)相關(guān)分析和循環(huán)重構(gòu)目的都是挖掘語句間的并行性，而數(shù)據(jù)劃分和處理器指派則重在策略，宏觀上挖掘并行性。Cannon算法是一種存儲有效的算法，設(shè)矩陣和相乘。為了使兩矩陣下標

7、滿足相乘的要求，和帶狀的并行分塊乘法不同，不是僅僅讓B矩陣的各列塊循環(huán)移動，而是有目的地讓A的各行塊以及B的各列塊皆施行循環(huán)移位，從而實現(xiàn)對C的子塊的計算。將矩陣A和B分成p個方塊Aij和Bij，每塊大小為，并將它們分配給個處理器。開始時處理器Pij存放塊Aij和Bij，并負責計算塊Cij，然后算法開始執(zhí)行： 將塊Aij向左循環(huán)移動i步；將塊Bij向上循環(huán)移動j步；Pij執(zhí)行乘加運算后將塊Aij向左循環(huán)移動1步，塊Bij向上循環(huán)移動1步；重復(fù)第步，總共執(zhí)行次乘加運算和次塊Aij和Bij的循環(huán)單步移位。二、Cannon乘法的MPI實現(xiàn)內(nèi)容和步驟實驗涉及內(nèi)容主要有：1）數(shù)據(jù)劃分和指派處理器最常用的

8、矩陣數(shù)據(jù)劃分有帶狀劃分和塊狀劃分。棋盤劃分比帶狀劃分可開發(fā)更高的并行度，Cannon乘法的mpi實現(xiàn)也正是基于棋盤劃分的并行實現(xiàn)。設(shè)有P個處理器，將矩陣A和B分成p個方塊Aij和Bij，每塊大小為，并將它們分配給個處理器。2） 子矩陣的循環(huán)移動處理器Pij存放塊Aij和Bij，并負責計算塊Cij，在使A矩陣的左右循環(huán)移動和B矩陣的上下循環(huán)移動時，為了避免在通信過程中發(fā)生死鎖，奇數(shù)號及偶數(shù)號處理器的收發(fā)順序被錯開，使偶數(shù)號處理器先發(fā)送后接收；而奇數(shù)號處理器先將子矩陣塊存于緩沖區(qū)Buffer中，然后接收編號在其后面的處理器所發(fā)送的子矩陣塊，最后再將緩沖區(qū)中子矩陣塊發(fā)送給編號在其前面的處理器。基本算

9、法如下：Begin(1)if (j=0) then /*最左端的子塊*/(1.1)將所存的A的子塊發(fā)送到同行最右端子塊所在的處理器中(1.2)接收其右鄰處理器中發(fā)來的A的子塊end if (2)if (j = sqrt(p)-1) and (j mod 2 = 0) then /*最右端子塊處理器且塊列號為偶數(shù)*/(2.1)將所存的A的子塊發(fā)送到其左鄰處理器中(2.2)接收其同行最左端子塊所在的處理器發(fā)來的A的子塊end if (3)if (j = sqrt(p)-1) and (j mod 2 0) then /*最右端子塊處理器且塊列號為奇數(shù)*/(3.1)將所存的A的子塊在緩沖區(qū)buffe

10、r中做備份(3.2)接收其同行最左端子塊所在的處理器發(fā)來的A的子塊(3.3)將在緩沖區(qū)buffer中所存的A的子塊發(fā)送到其左鄰處理器中end if (4)if (j sqrt(p)-1) and (j mod 2 = 0) and (j 0) then /*其余的偶數(shù)號處理器*/(4.1)將所存的A的子塊發(fā)送到其左鄰處理器中(4.2)接收其右鄰處理器中發(fā)來的A的子塊end if (5)if (j sqrt(p)-1) and (j mod 2 = 1) and (j 0) then /*其余的奇數(shù)號處理器*/(5.1)將所存的A的子塊在緩沖區(qū)buffer中做備份(5.2)接收其右鄰處理器中發(fā)來

11、的A的子塊(5.3)將在緩沖區(qū)buffer中所存的A的子塊發(fā)送到其左鄰處理器中end if End實驗步驟1） 登陸KD-60圖 2.1 KD-60登陸界面2） 轉(zhuǎn)至node80節(jié)點，上傳程序輸入命令：ssh loongsonnode80 和密碼 進入圖界面圖2.2 轉(zhuǎn)到節(jié)點80的界面再命令vim，進入vim編輯器加入程序,保存為cannon.c 3） 編譯程序輸入命令：mpicc cannon.c o cannon lm 在目錄中查看，已成功。如下圖圖2.3 將程序保存并編譯后界面4） 運行程序輸入：mpirun np 4 cannon 4 ,其中第一個4是指定的處理器個數(shù)，第二個4是產(chǎn)生隨

12、機矩陣的維數(shù)，這兩個參數(shù)在實驗過程中可以調(diào)整，但要求第一個參數(shù)即處理器的個數(shù)必須是一個數(shù)的平方數(shù)。輸出：圖2.4 cannon乘法運行結(jié)果圖 2.4 并行程序運行界面 兩個參數(shù)都是4， 分別輸出兩個隨機矩陣和矩陣的乘積三、數(shù)據(jù)及結(jié)果1.下面列出了兩組數(shù)據(jù)，分別是用一個處理器進行串行運算和四個處理器進行并行運算矩陣維數(shù)為200的計算時間比較。四個處理器處理階數(shù)為200的矩陣相乘時，所花時間為：1.705844秒。單個處理器處理階數(shù)為200的矩陣相乘時，所花時間為：3.727210秒。如圖3.1 和圖3.2所示。圖3.1 四個處理器并行執(zhí)行結(jié)果圖圖3.2 單個處理器串行執(zhí)行結(jié)果圖附：1. 程序模塊

13、偽代碼：輸入：Ann，Bnn輸出：CnnBegin對所有處理器my_rank(my_rank=0,p-1)同時執(zhí)行如下的算法:(1)計算子塊的行號 i=my_rank/sqrt(p)計算子塊的列號 j=my_rank mod sqrt(p)(2)for k=0 to -1 doif (ik) then Leftmoveonestep(a) end if /* a循環(huán)左移至同行相鄰處理器中*/if (jk) then Upmoveonestep(b) end if /* b循環(huán)上移至同列相鄰處理器中*/ end for (3)for i=0 to m-1 dofor j=0 to m-1 doc

14、i,j=0end forend for (4)for k=0 to -1 dofor i=0 to m-1 dofor j=0 to m-1 dofor k1=0 to m-1 doci,j= ci,j+ ai,k1* bk1,jend forend forend for Leftmoveonestep(a) /*子塊a循環(huán)左移至同行相鄰的處理器中*/Upmoveonestep(b) /*子塊b循環(huán)上移至同列相鄰的處理器中*/end for EndLeftmoveonestep(a) 見實驗內(nèi)容處附2.程序源碼#include #include #include #include #inclu

15、de #include /* 全局變量聲明 */float *A, *B, *C; /* 總矩陣,C = A * B */float *a, *b, *c, *tmp_a, *tmp_b; /* a、b、c表分塊，tmp_a、tmp_b表緩沖區(qū) */int dg, dl, dl2,p, sp; /* dg:總矩陣維數(shù);dl:矩陣塊維數(shù);dl2=dl*dl;p:處理器個數(shù);spsqrt(p) */int my_rank, my_row, my_col; /* my_rank:處理器ID;(my_row,my_col):處理器邏輯陣列坐標 */MPI_Status status;float sta

16、rttime;float time1;/* *函數(shù)名: get_index *功能：處理器邏輯陣列坐標至rank號的轉(zhuǎn)換 *輸入：坐標、邏輯陣列維數(shù) *輸出：rank號 */int get_index(int row, int col, int sp) return (row+sp)%sp)*sp + (col+sp)%sp;/* *函數(shù)名：random_A_B *功能：隨機生成矩陣A和B */void random_A_B() int i,j; srand(unsigned int)time(NULL); /*設(shè)隨機數(shù)種子*/*隨機生成A,B,并初始化C*/ for(i=0; idg ; i

17、+) for(j=0; jdg ; j+) Aij = rand(); Bij = rand(); Cij = 0.0; /* 函數(shù)名：scatter_A_B * 功能：rank為0的處理器向其他處理器發(fā)送A、B矩陣的相關(guān)塊 */void scatter_A_B() int i,j,k,l; int p_imin,p_imax,p_jmin,p_jmax; for(k=0; kp; k+) /*計算相應(yīng)處理器所分得的矩陣塊在總矩陣中的坐標范圍*/ p_jmin = (k % sp ) * dl; p_jmax = (k % sp + 1) * dl-1; p_imin = (k - (k %

18、sp)/sp * dl; p_imax = (k - (k % sp)/sp +1) *dl -1; l = 0; /*rank=0的處理器將A,B中的相應(yīng)塊拷至tmp_a,tmp_b，準備向其他處理器發(fā)送*/ for(i=p_imin; i=p_imax; i+) for(j=p_jmin; j=p_jmax; j+) tmp_al = Aij; tmp_bl = Bij; l+; /*rank=0的處理器直接將自己對應(yīng)的矩陣塊從tmp_a,tmp_b拷至a,b*/ if(k=0) memcpy(a, tmp_a, dl2 * sizeof(float); memcpy(b, tmp_b,

19、dl2 * sizeof(float); else /*rank=0的處理器向其他處理器發(fā)送tmp_a,tmp_b中相關(guān)的矩陣塊*/ MPI_Send(tmp_a, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD); MPI_Send(tmp_b, dl2, MPI_FLOAT, k, 2, MPI_COMM_WORLD); /* *函數(shù)名:init_alignment *功能:矩陣A和B初始對準 */void init_alignment() /*將A中坐標為(i,j)的分塊A(i,j)向左循環(huán)移動i步*/ MPI_Sendrecv(a, dl2, MPI_FLOA

20、T, get_index(my_row,my_col-my_row,sp), 1, tmp_a, dl2, MPI_FLOAT, get_index(my_row,my_col+my_row,sp), 1, MPI_COMM_WORLD, &status); memcpy(a, tmp_a, dl2 * sizeof(float) ); /*將B中坐標為(i,j)的分塊B(i,j)向上循環(huán)移動j步*/ MPI_Sendrecv(b, dl2, MPI_FLOAT, get_index(my_row-my_col,my_col,sp), 1, tmp_b, dl2, MPI_FLOAT, get

21、_index(my_row+my_col,my_col,sp), 1, MPI_COMM_WORLD, &status); memcpy(b, tmp_b, dl2 * sizeof(float) );/* *函數(shù)名：main_shift *功能：分塊矩陣左移和上移，并計算分塊c */void main_shift() int i,j,k,l; for(l=0; lsp; l+) /*矩陣塊相乘，c+=a*b */ for(i=0; idl; i+) for(j=0; jdl; j+) for(k=0; kdl; k+) ci*dl+j += ai*dl+k*bk*dl+j; /* 將分塊a左

22、移1位 */ MPI_Send(a , dl2, MPI_FLOAT, get_index(my_row, my_col-1, sp), 1, MPI_COMM_WORLD); MPI_Recv(a , dl2, MPI_FLOAT, get_index(my_row, my_col+1, sp), 1, MPI_COMM_WORLD, &status); /* 將分塊b上移1位 */ MPI_Send(b , dl2, MPI_FLOAT, get_index(my_row-1, my_col, sp), 1, MPI_COMM_WORLD); MPI_Recv(b , dl2, MPI_F

23、LOAT, get_index(my_row+1, my_col, sp), 1, MPI_COMM_WORLD, &status); /* *函數(shù)名：collect_c *功能：rank為0的處理器從其余處理器收集分塊矩陣c */void collect_C() int i,j,i2,j2,k; int p_imin,p_imax,p_jmin,p_jmax; /* 分塊矩陣在總矩陣中頂點邊界值 */ /* 將rank為0的處理器中分塊矩陣c結(jié)果賦給總矩陣C對應(yīng)位置 */ for (i=0;idl;i+) for(j=0;jdl;j+) Cij=ci*dl+j; for (k=1;kp;k+

24、) /*將rank為0的處理器從其他處理器接收相應(yīng)的分塊c*/ MPI_Recv(c, dl2, MPI_FLOAT, k, 1, MPI_COMM_WORLD, &status); p_jmin = (k % sp ) *dl; p_jmax = (k % sp + 1) *dl-1; p_imin = (k - (k % sp)/sp *dl; p_imax = (k - (k % sp)/sp +1) *dl -1; i2=0; /*將接收到的c拷至C中的相應(yīng)位置,從而構(gòu)造出C*/ for(i=p_imin; i=p_imax; i+) j2=0; for(j=p_jmin; j=p_j

25、max; j+) Cij=ci2*dl+j2; j2+; i2+; /*函數(shù)名：print *功能：打印矩陣 *輸入：指向矩陣指針的指針，字符串 */void print(float *m,char *str) int i,j; printf(%s,str); /*打印矩陣m*/ for(i=0;idg;i+) for(j=0;jdg;j+) printf(%15.0f ,mij); printf(n); printf(n);/* *函數(shù)名：main *功能：主過程，Cannon算法，矩陣相乘 *輸入：argc為命令行參數(shù)個數(shù)，argv為每個命令行參數(shù)組成的字符串數(shù)組 */int main(i

26、nt argc, char *argv) int i; MPI_Init(&argc, &argv); /* 啟動MPI計算 */ MPI_Comm_size(MPI_COMM_WORLD, &p); /* 確定處理器個數(shù) */ MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); /* 確定各自的處理器標識符 */ sp = sqrt(p); /* 確保處理器個數(shù)是完全平方數(shù)，否則打印錯誤信息，程序退出 */ if (sp*sp != p) if (my_rank = 0) printf(Number of processors is not a quadrati

27、c number!n); MPI_Finalize(); exit(1); if (argc != 2) if (my_rank = 0) printf(usage: mpirun -np ProcNum cannon MatrixDimensionn); MPI_Finalize(); exit(1); dg = atoi(argv1); /* 總矩陣維數(shù) */ dl = dg / sp; /* 計算分塊矩陣維數(shù) */ dl2 = dl * dl; /* 計算處理器在邏輯陣列中的坐標 */ my_col = my_rank % sp ; my_row = (my_rank-my_col) /

28、 sp ; /* 為a、b、c分配空間 */ a = (float *)malloc( dl2 * sizeof(float) ); b = (float *)malloc( dl2 * sizeof(float) ); c = (float *)malloc( dl2 * sizeof(float) ); /* 初始化c */ for(i=0; idl2 ; i+) ci = 0.0; /* 為tmp_a、tmp_b分配空間 */ tmp_a = (float *)malloc( dl2 * sizeof(float) ); tmp_b = (float *)malloc( dl2 * si

29、zeof(float) ); if (my_rank = 0) /* rank為0的處理器為A、B、C分配空間 */ starttime = MPI_Wtime(); A = (float *)malloc( dg * sizeof(float*) ); B = (float *)malloc( dg * sizeof(float*) ); C = (float *)malloc( dg * sizeof(float*) ); for(i=0; idg; i+) Ai = (float *)malloc( dg * sizeof(float) ); Bi = (float *)malloc(

30、dg * sizeof(float) ); Ci = (float *)malloc( dg * sizeof(float) ); random_A_B(); /* rank為0的處理器隨機化生成A、B矩陣 */ scatter_A_B(); /* rank為0的處理器向其他處理器發(fā)送A、B矩陣的相關(guān)塊 */ else /* rank不為0的處理器接收來自rank為0的處理器的相應(yīng)矩陣分塊 */ MPI_Recv(a, dl2, MPI_FLOAT, 0 , 1, MPI_COMM_WORLD, &status); MPI_Recv(b, dl2, MPI_FLOAT, 0 , 2, MPI_

31、COMM_WORLD, &status); init_alignment(); /* A、B矩陣的初始對準 */ main_shift(); /* 分塊矩陣左移、上移, cannon算法的主過程 */ if(my_rank = 0) collect_C(); /* rank為0的處理器從其余處理器收集分塊矩陣c */ print(A,random matrix A : n); /* 打印矩陣A */ print(B,random matrix B : n); /* 打印矩陣B */ print(C,Matrix C = A * B : n); /* 打印矩陣C */ else MPI_Send

32、(c,dl2,MPI_FLOAT,0,1,MPI_COMM_WORLD); /* rank不為0的處理器向rank為0的處理器發(fā)送矩陣塊c */ MPI_Barrier(MPI_COMM_WORLD); /* 同步所有處理器 */if(my_rank = 0) time1 = MPI_Wtime();printf(%s,time1-starttime); /* 運行時間*/ MPI_Finalize(); /* 結(jié)束MPI計算 */ return 0;四、結(jié)果分析和總結(jié)由第三部分數(shù)據(jù)結(jié)果可以看出，多個處理器并行算法比單個處理器串行程序在計算大量數(shù)據(jù)時要快得多。分析如下：設(shè)，相乘，結(jié)果矩陣：，其中，記一次乘法和加法為一個運算時間單位，則這種情況下串行運行時間為個時間單位。在簡單的帶狀劃分中，A