并行计算论文

*******************实践教学*******************兰州理工大学理学院2016年春季学期并行计算课程设计专业班级：13级信息与计算科学*名：***学号：********指导教师：孟新友郭秀婷成绩：摘要本文主要设计了MIMD（多指令流多数据）异步并行的Gauss一Seidel迭代法，并利用MPI消息传递模型实现了该算法并实际测试了求解所花费的时间和进行了算法的理论评估。

算法的主要过程如下：首先，主进程在数据域对线性方程组的增广矩阵进行带状化分，并将划分的结果广播到其他处理器，然后，各个处理器按照自己的指令流进行计算，并利用存储转发的方式进行消息传递，最后，满足精度要求后，通过设置同步路障使得各个处理器的结果归约至主进程，由主进程并输出求解结果。

关键词：异步并行Gauss一Seidel迭代法消息传递同步路障目录摘要 (2)一．题目及要求 (4)1.1题目 (4)1.2要求 (4)二．算法设计原理 (4)2.1设计算法 (4)2.2算法原理 (5)三．算法描述及设计流程 (6)3.1算法描述 (6)3.2设计流程 (7)四．源程序代码及运行结果 (9)4.1源程序代码 (9)4.2运行结果 (16)五．算法分析及优缺点 (17)5.1算法分析 (17)5.2优缺点 (17)总结 (18)参考文献 (19)一．题目及要求1.1题目迭代求解的高斯-赛德尔法（MIMD 异步并行算法）求解方程组⎪⎩⎪⎨⎧-=+-=+-=++1352422532132321x x x x x x x x 的解 1.2要求本次课程设计的题目是利用高斯-赛德迭尔求解方程组，在求解过程中要求用MIMD 异步并行算法求解。

二．算法设计原理2.1设计算法求解线性方程组的高斯—赛德尔法（Gauss-Seidel Method ），实际上是迭代法，不仅稠密矩阵性系适用，而且稀疏矩阵性系也适用。

下面是他的主要原理：在求解错误!未找到引用源。

,可以将系数矩阵错误!未找到引用源。

分解为错误!未找到引用源。

,其中错误!未找到引用源。

均为错误!未找到引用源。

的矩阵，具体的定义如下：[2]这样，错误!未找到引用源。

可以换成错误!未找到引用源。

如果给定初始解错误!未找到引用源。

,则第错误!未找到引用源。

次迭代可计算如下：错误!未找到引用源。

(1)由于在迭代的算法中需要判断是否在第错误!未找到引用源。

次收敛，这里利用向量的1范数去判断是否终止迭代过程，计算公式如下：错误!未找到引用源。

(2)(1)式能够加快顺序计算速度，因此当依次计算错误!未找到引用源。

时，第错误!未找到引用源。

次迭代的错误!未找到引用源。

值，一部分可以用本次迭代的值（相应于上三角部分），而一部分可用本次迭代的值（相应于下三角部分）。

2.2算法原理对于以上的分析，该算法无法直接并行化，为此通过分析，将线性方程组错误!未找到引用源。

的增广矩阵错误!未找到引用源。

按行划分如图 1。

图 1 增广矩阵的划分对于高斯-塞德尔迭代，计算i x 的新值时，使用11,,i n x x +-⋯的旧值和01,,i x x -⋯的新值。

计算过程中i x 与01,,i x x -⋯及11,,i n x x +-⋯的新值会在不同的处理器中产生，因此可以考虑采用时间偏移的方法，使各个处理器对新值计算的开始和结束时间产生一定的偏差。

编号为 my_rank 的处理器一旦计算出(_(_1))i x my rank m i my rank m ⨯≤+⨯的新值，就立即广播给其余处理器，以供各处理器对x 的其它分量计算有关i x 的乘积项并求和。

当它计算完 x 的所有分量后，它还要接收其它处理器发送的新的x 分量，并对这些分量进行求和计算，为计算下一轮的i x 作准备。

计算开始时，所有处理器并行地对主对角元素右边的数据项进行求和，此时编号为 0 的处理器 （简称为0p ）计算出0x 然后广播给其余处理器，其余所有的处理器用0x 的新值和其对应项进行求和计算，接着0p 计算出12,,,x x ⋯当0p 完成对1m x -的计算和广播后，1p 计算出m x ，并广播给其余处理器，其余所有的处理器用m x 的新值求其对应项的乘积并作求和计算。

然后1p计算出12,,,m m x x ++⋯当1p 完成对2*1m x -的计算和广播后，2p 计算出2*m x ⋯，如此重复下去，直至1n x - 在1p p - 中被计算出并广播至其余的处理器之后，0p 计算出下一轮的新的0x ，这样逐次迭代下去，直至收敛为止。

[1]三．算法描述及设计流程3.1算法描述通过在第二小节的分析，我们才用类程序语言对MIMD 异步并行算法的Guass 一Seidel 法进行描述：输入：输入矩阵n n A ⨯，常数项1⨯n b ，初始解())0(1)0(1)0(0)0(,...,,-=n x x x x ，精度ε。

输出：解向量()n x x x ,...,1,0。

对所有处理器my_rank(my_rank=)同时执行如for(i=my_rank m;(my_rank+1) m ；i++)//所有处理器并行的对主对角线元素右边的数据求和sum[i]=0.0for(j=i+1;j<n-1;j++)sum[i]=sum[i]+A[i,j]*x[j]total=0endforendforwhile(total<n){//total 为新旧值之差小于 的 分量个数（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)//表示朱行所在的处理器//产生的新值 Temp=x[j],x[j]=(b[j]-sum[j])/a[I,j]EndifIf (|x[j]-temp|<)Iteration++Endif将x[j]的新值广播到其他的n-1个处理器S[j]=0For(i=my_rank*m;i<=(my_rank+1)*m;i++)If(j)Sum[i]=sum[i]+A[I,j]*x[j]Else //其他处理器接受广播来的x[j]新值endif//对所存在各行计算x[j]所对应的内积项累加For（i=my_rank*m;i<=((my_rank+1)*m-1）;i++）Sum[i]=sum[i]+A[I,j]*x[j]endforendfor求出所有处理器中iteration和total的值并广播到所有处理器endforendwhile3.2设计流程为了实现MIMD异步并行算法的Guass一Seidel法，采用的消息传模型进行实现。

具体的实现流程如图2图2.程序设计流程消息传递中所使用的相关函数如表1 函数名 功能MPI_Init(&argc,&argv) 启动MPI MPI_Comm_rank(MPI_COMM_WORLD, int *rank) 获取进程编号 MPI_Comm_size(MPI_COMM_WORLD,int *size) 获取进程数 MPI_Bcast(Address,Count,Datatype,Root ，COMM) 广播 MPI_Send(void *buf, int count, MPI_datatype,intdest,inttag,MPI_COMM_WORLD)消息发送 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)全局规约函数 MPI_Barrier(MPI_COMM_WORLD) 路障 MPI_Finalize() MPI 结束表1 本算法中使用的ＭＰＩ函数四．源程序代码及运行结果4.1源程序代码#include "stdio.h"#include "stdlib.h"#include "mpi.h"#include "math.h"#define E 0.001#define a(x,y) a[x*size+y]#define b(x) b[x]#define v(x) v[x]#define A(x,y) A[x*size+y]#define B(x) B[x]#define V(x) V[x]#define intsize sizeof(int)#define floatsize sizeof(float)#define 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,group_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){starttime=MPI_Wtime();fdA=fopen("dataIn.txt","r");fscanf(fdA,"%d %d", &size, &N);if (size != N-1){printf("the input is wrong\n");exit(1);}A=(float *)malloc(floatsize*size*size);B=(float *)malloc(floatsize*size);V=(float *)malloc(floatsize*size);for(i = 0; i < size; i++){for(j = 0; j < N; j++){if(j==N-1)fscanf(fdA,"%f", B+i);elsefscanf(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("%d\t%d\n", size, N);for(i = 0; i < size; i ++){for(j = 0; j < size; j ++)printf("%f\t",A(i,j));printf("%f\n",B(i));}printf("\n");for(i = 0; i < size; i ++){printf("%f\t", V(i));}printf("\n\n");printf("\nOutput of result\n");}MPI_Bcast(&size,1,MPI_INT,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){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!");elseprintf("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); MPI_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,j));}printf("\n");}for(i=0;i<m;i++){sum[i]=0.0;for(j=0;j<size;j++){if (j>(my_rank*m+i))sum[i]=sum[i]+a(i,j)*v(j);}printf("This is %f\n",sum[i]);}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++)sum[r]=sum[r]+a(r,my_rank*m+i)*v(my_rank*m+i); v(my_rank*m+r)=(1-w)*temp+w*(b(r)-sum[r])/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);sum[r]=0.0;for(i=0;i<r;i++)sum[i]=sum[i]+a(i,my_rank*m+r)*v(my_rank*m+r);}else{MPI_Bcast(&v(q*m+r),1,MPI_FLOAT,q,MPI_COMM_WORLD);[4]for(i=0;i<m;i++)sum[i]=sum[i]+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 = %d\n",loop,total);}time2=MPI_Wtime();if (my_rank==0){for(i = 0; i < size; i ++) printf("x[%d] = %f\n",i,v(i));printf("\n");printf("Iteration num = %d\n",loop);printf("Whole running time = %f seconds\n",time2-starttime); printf("Distribute data time = %f seconds\n",time1-starttime); printf("Parallel compute time = %f seconds\n",time2-time1); }MPI_Barrier(MPI_COMM_WORLD);MPI_Finalize();Environment_Finalize(a,b,v);return 0;}4.2运行结果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.000000 0.000000 0.000000 1.000000Output 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 = 3x[0] = -1.999635x[1] = 1.999954x[2] = 0.999866Iteration num = 7Whole running time = 0.004000 secondsDistribute data time = 0.003000 secondsParallel compute time = 0.001000 seconds五．算法分析及优缺点5.1算法分析若取一次乘法和加法运算时间或一次比较运算时间为一个单位时间。