并行計算論文

1、*實踐教學*蘭州理工大學 理學院 2016年春季學期 并行計算 課程設(shè)計專業(yè)班級： 13級信息與計算科學 姓 名： 田 興 福 學 號： 13540231 指導教師： 孟新友 郭秀婷 成 績： 摘要本文主要設(shè)計了MIMD（多指令流多數(shù)據(jù)）異步并行的Gauss一Seidel迭代法，并利用MPI消息傳遞模型實現(xiàn)了該算法并實際測試了求解所花費的時間和進行了算法的理論評估。算法的主要過程如下：首先，主進程在數(shù)據(jù)域?qū)€性方程組的增廣矩陣進行帶狀化分，并將劃分的結(jié)果廣播到其他處理器，然后，各個處理器按照自己的指令流進行計算，并利用存儲轉(zhuǎn)發(fā)的方式進行消息傳遞，最后，滿足精度要求后，通過設(shè)置同步路障使得各個處

2、理器的結(jié)果歸約至主進程，由主進程并輸出求解結(jié)果。關(guān)鍵詞：異步并行 Gauss一Seidel迭代法 消息傳遞 同步 路障目錄摘要2一題目及要求41.1題目41.2要求4二算法設(shè)計原理42.1設(shè)計算法42.2算法原理5三算法描述及設(shè)計流程63.1算法描述63.2設(shè)計流程7四源程序代碼及運行結(jié)果94.1源程序代碼94.2運行結(jié)果16五算法分析及優(yōu)缺點175.1算法分析175.2優(yōu)缺點17總結(jié)18參考文獻19一題目及要求1.1題目迭代求解的高斯-賽德爾法（MIMD異步并行算法）求解方程組的解1.2要求本次課程設(shè)計的題目是利用高斯-賽德迭爾求解方程組，在求解過程中要求用MIMD異步并行算法求解。二算法設(shè)

3、計原理2.1設(shè)計算法求解線性方程組的高斯賽德爾法（Gauss-Seidel Method），實際上是迭代法，不僅稠密矩陣性系適用，而且稀疏矩陣性系也適用。下面是他的主要原理：在求解,可以將系數(shù)矩陣分解為,其中均為的矩陣，具體的定義如下：2這樣，可以換成。如果給定初始解,則第次迭代可計算如下： (1)由于在迭代的算法中需要判斷是否在第次收斂，這里利用向量的1范數(shù)去判斷是否終止迭代過程，計算公式如下： (2)(1) 式能夠加快順序計算速度，因此當依次計算時，第次迭代的值，一部分可以用本次迭代的值（相應于上三角部分），而一部分可用本次迭代的值（相應于下三角部分）。2.2算法原理對于以上的分析，該算法

4、無法直接并行化，為此通過分析，將線性方程組的增廣矩陣按行劃分如圖 1。圖 1 增廣矩陣的劃分對于高斯-塞德爾迭代，計算的新值時，使用的舊值和的新值。計算過程中與及的新值會在不同的處理器中產(chǎn)生，因此可以考慮采用時間偏移的方法，使各個處理器對新值計算的開始和結(jié)束時間產(chǎn)生一定的偏差。編號為 my_rank 的處理器一旦計算出的新值，就立即廣播給其余處理器，以供各處理器對x 的其它分量計算有關(guān)的乘積項并求和。當它計算完 x的所有分量后，它還要接收其它處理器發(fā)送的新的分量，并對這些分量進行求和計算，為計算下一輪的作準備。計算開始時，所有處理器并行地對主對角元素右邊的數(shù)據(jù)項進行求和，此時編號為 0 的處理

5、器 （簡稱為）計算出然后廣播給其余處理器，其余所有的處理器用的新值和其對應項進行求和計算，接著計算出當完成對的計算和廣播后，計算出，并廣播給其余處理器，其余所有的處理器用的新值求其對應項的乘積并作求和計算。 然后計算出當完成對的計算和廣播后，計算出，如此重復下去，直至 在 中被計算出并廣播至其余的處理器之后，計算出下一輪的新的，這樣逐次迭代下去，直至收斂為止。1三算法描述及設(shè)計流程3.1算法描述通過在第二小節(jié)的分析，我們才用類程序語言對MIMD異步并行算法的Guass一Seidel法進行描述：輸入：輸入矩陣，常數(shù)項，初始解，精度。輸出：解向量。對所有處理器my_rank(my_rank=)同時

6、執(zhí)行如for(i=my_rank m;(my_rank+1) m；i+)/所有處理器并行的對主對角線元素右邊的數(shù)據(jù)求和sumi=0.0for(j=i+1;j<n-1;j+)sumi=sumi+Ai,j*xjtotal=0endforendforwhile(total<n)/total為新舊值之差小于 的 分量個數(shù)（2.1）iteration=0/ iteration為本處理器中新舊值之差小于 的分量個（2.2）for(j=0;j<=n-1;j+)/依次以第0,1,n-1行為主行 (i)q=j/m (ii)if(my_rank=q)/表示朱行所在的處理器 /產(chǎn)生的新值Temp=

7、xj,xj=(bj-sumj)/aI,jEndifIf (|xj-temp|<)Iteration+Endif將xj的新值廣播到其他的n-1個處理器Sj=0For(i=my_rank*m;i<=(my_rank+1)*m;i+)If(j)Sumi=sumi+AI,j*xjElse /其他處理器接受廣播來的xj新值endif/對所存在各行計算xj所對應的內(nèi)積項累加For（i=my_rank*m;i<=(my_rank+1)*m-1）;i+）Sumi=sumi+AI,j*xjendforendfor求出所有處理器中iteration和total的值并廣播到所有處理器endfore

8、ndwhile3.2設(shè)計流程為了實現(xiàn)MIMD異步并行算法的Guass一Seidel法，采用的消息傳模型進行實現(xiàn)。具體的實現(xiàn)流程如圖2圖2.程序設(shè)計流程消息傳遞中所使用的相關(guān)函數(shù)如表1函數(shù)名功能MPI_Init(&argc,&argv)啟動MPIMPI_Comm_rank(MPI_COMM_WORLD, int *rank)獲取進程編號MPI_Comm_size(MPI_COMM_WORLD,int *size)獲取進程數(shù)MPI_Bcast(Address,Count,Datatype,Root，COMM)廣播MPI_Send(void *buf, int count, MPI_

9、datatype,intdest,inttag,MPI_COMM_WORLD)消息發(fā)送MPI_Recv(void *buf, int count, MPI_datatype,intsource,inttag,MPI_COMM_WORLD,MPI_Status *status)消息接收MPI_Allreduce(&iteration,&total,1, MPI_FLOAT,MPI_SUM, MPI_COMM_WORLD)全局規(guī)約函數(shù)MPI_Barrier(MPI_COMM_WORLD)路障MPI_Finalize()MPI結(jié)束表1 本算法中使用的函數(shù)四源程序代碼及運行結(jié)果4.1源

10、程序代碼#include "stdio.h"#include "stdlib.h"#include "mpi.h"#include "math.h"#define E 0.001#define a(x,y) ax*size+y#define b(x) bx#define v(x) vx#define A(x,y) Ax*size+y#define B(x) Bx#define V(x) Vx#define intsize sizeof(int)#define floatsize sizeof(float)#defi

11、ne charsize sizeof(char)int size,N; int m;float *B;float *A;float *V; double starttime;double time1;double time2;int my_rank;int p;MPI_Status status;FILE *fdA,*fdB,*fdB1;void Environment_Finalize(float *a,float *b,float *v) free(a); free(b); free(v);int main(int argc, char *argv) int i,j,my_rank,gro

12、up_size; int k; float *sum; float *b; float *v; float *a; float *differ; float temp,w; int iteration,total,loop; int r,q; loop=0; total=0; w=0.9; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD,&group_size); MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); p=group_size; if(my_rank=0) star

13、ttime=MPI_Wtime(); fdA=fopen("dataIn.txt","r"); fscanf(fdA,"%d %d", &size, &N); if (size != N-1) printf("the input is wrongn"); exit(1); A=(float *)malloc(floatsize*size*size); B=(float *)malloc(floatsize*size); V=(float *)malloc(floatsize*size);for(i

14、= 0; i < size; i+) for(j = 0; j < N; j+) if(j=N-1) fscanf(fdA,"%f", B+i); else fscanf(fdA,"%f", A+i*size+j); for(i = 0; i < size; i +) fscanf(fdA,"%f", V+i); fclose(fdA); printf("Input of file "dataIn.txt"n"); printf("%dt%dn", size

15、, N); for(i = 0; i < size; i +) for(j = 0; j < size; j +) printf("%ft",A(i,j); printf("%fn",B(i); printf("n"); for(i = 0; i < size; i +) printf("%ft", V(i); printf("nn"); printf("nOutput of resultn"); MPI_Bcast(&size,1,MPI_INT,

16、0,MPI_COMM_WORLD); m=size/p;if (size%p!=0) m+; v=(float *)malloc(floatsize*size); a=(float *)malloc(floatsize*m*size); b=(float *)malloc(floatsize*m); sum=(float *)malloc(floatsize*m); if (my_rank=0) for(i=0;i<m;i+) for(j=0;j<size;j+) a(i,j)=A(i,j); for(i=0;i<m;i+) b(i)=B(i); if (my_rank=0)

17、 for(i=0;i<size;i+) v(i)=V(i); MPI_Bcast(v,size,MPI_FLOAT,0,MPI_COMM_WORLD); if (a=NULL|b=NULL|v=NULL) printf("allocate space fail!"); else printf("allocate space successful!n"); if (my_rank=0) for(i=1;i<p;i+) MPI_Send(&(A(m*i,0),m*size,MPI_FLOAT,i,i,MPI_COMM_WORLD); MP

18、I_Send(&(B(m*i),m,MPI_FLOAT,i,i,MPI_COMM_WORLD); free(A); free(B); free(V); else MPI_Recv(a,m*size,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD,&status); MPI_Recv(b,m,MPI_FLOAT,0,my_rank,MPI_COMM_WORLD,&status); time1=MPI_Wtime(); for(i=0;i<m;i+) for(j=0;j<size;j+) printf("%f ",a(i

19、,j); printf("n"); for(i=0;i<m;i+) sumi=0.0; for(j=0;j<size;j+) if (j>(my_rank*m+i) sumi=sumi+a(i,j)*v(j); printf("This is %fn",sumi); while (total<size) iteration=0; total=0; for(j=0;j<size;j+) r=j%m; q=j/m; if (my_rank=q) temp=v(my_rank*m+r); for(i=0;i<r;i+) su

20、mr=sumr+a(r,my_rank*m+i)*v(my_rank*m+i); v(my_rank*m+r)=(1-w)*temp+w*(b(r)-sumr)/a(r,my_rank*m+r); if (fabs(v(my_rank*m+r)-temp)<E) iteration+; MPI_Bcast(&v(my_rank*m+r),1,MPI_FLOAT,my_rank,MPI_COMM_WORLD); sumr=0.0; for(i=0;i<r;i+) sumi=sumi+a(i,my_rank*m+r)*v(my_rank*m+r); else MPI_Bcast

21、(&v(q*m+r),1,MPI_FLOAT,q,MPI_COMM_WORLD);4 for(i=0;i<m;i+) sumi=sumi+a(i,q*m+r)*v(q*m+r); MPI_Allreduce(&iteration,&total,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD);loop+; if (my_rank=0) printf("%-2d th total vaule = %dn",loop,total); time2=MPI_Wtime(); if (my_rank=0) for(i = 0; i &l

22、t; size; i +) printf("x%d = %fn",i,v(i); printf("n"); printf("Iteration num = %dn",loop); printf("Whole running time = %f secondsn",time2-starttime); printf("Distribute data time = %f secondsn",time1-starttime); printf("Parallel compute time = %

23、f secondsn",time2-time1); MPI_Barrier(MPI_COMM_WORLD); MPI_Finalize(); Environment_Finalize(a,b,v);return 0;4.2運行結(jié)果Input of file "dataIn.txt"3 45.000000 2.000000 2.000000 -4.0000000.000000 2.000000 1.000000 5.0000001.000000 -1.000000 3.000000 -1.0000000.000000 0.000000 1.000000Output

24、of resultallocate space successful!This is 2.0000001 th total vaule = 02 th total vaule = 03 th total vaule = 04 th total vaule = 05 th total vaule = 06 th total vaule = 07 th total vaule = 3x0 = -1.999635x1 = 1.999954x2 = 0.999866Iteration num = 7Whole running time = 0.004000 secondsDistribute data time = 0.003000 secondsParallel compute time = 0.001000 seconds五算法分析及優(yōu)缺點5.1算法分析若取一次乘法和加法運算時間或一次比較運算時間為一個單位時間。在算法開始階段，各處理器并行計算主對角線右邊元素相應的乘積項并求和,所需時間為，進入迭代計算后，雖然各個處理器所負責的的分量在每一輪計算中的開始時間是不一樣的，但一輪迭代中的計算量都是相等的，因此不妨取0 號處理器為對象分析一輪迭代計算的時間，容易得出0號處理器計算時間為；另外它在一輪迭代中做廣播操作 次，通信量為