计算流体力学中科院力学所第6讲-MPI并行程序设计初步.ppt

《计算流体力学中科院力学所第6讲-MPI并行程序设计初步.ppt》由会员分享，可在线阅读，更多相关《计算流体力学中科院力学所第6讲-MPI并行程序设计初步.ppt（54页珍藏版）》请在三一办公上搜索。

1、计算流体力学讲义 第六讲 MPI并行程序设计（2）李新亮；力学所主楼219；82543801,知识点：阻塞通信与非阻塞通信 非连续数据的发送与接收,1,讲义、课件上传至(流体中文网）-“流体论坛”-“CFD基础理论”,Copyright by Li Xinliang,Part 2 MPI 并行程序设计 1.阻塞式通信与非阻塞式通信 2.非连续数据的传送与接收 3.MPI的通信域与组；算例计算三维差分 4.分布式数据的存储Part 3 实例教学 1.谱方法求解不可压缩N-S方程 2.利用紧致型差分格式求导数,2,Copyright by Li Xinliang,MPI的消息发送机制 两步进行MP

2、I_Send(A,)发送MPI_Recv(B,)接收,发送 变量A,接收 到变量B,配合使用,3,Copyright by Li Xinliang,一、阻塞式通信与非阻塞式通信,阻塞式发送与接收,MPI_Send(A,),MPI_Recv(B,),4,Copyright by Li Xinliang,MPI_Send()返回后缓冲区可释放 sum=call MPI_Send(sum,)sum=变量可重复利用 MPI_Recv()返回后缓冲区数据可使用Call MPI_Recv(sum1,)Sum=sum0+sum1,5,Copyright by Li Xinliang,非阻塞发送,启动发送,立

3、即返回,计算,通信完成,释放发送缓冲区,发送消息,非阻塞接收,启动接收,立即返回,计算,通信完成,引用接收数据,接收消息,计算与通信重叠,非阻塞消息发送与接收,6,Copyright by Li Xinliang,非阻塞消息发送MPI_ISend(buf,count,datatype,dest,tag,comm,request,ierr)In buf,count,datatype,dest,tag,commOut request,ierr Request(返回的非阻塞通信对象,整数)非阻塞消息接收MPI_IRecv(buf,count,datatype,source,tag,comm,requ

4、est,ierr)In buf,count,datatype,source,tag,commOut request,ierr非阻塞通信的完成 MPI_Wait(request,status,ierr)等待消息收发完成 MPI_Test(request,flag,stutus,ierr)MPI_Waitall(const,request_array,status,ierr)等待多个消息完成 In requestOut status，flag(logical型),7,Copyright by Li Xinliang,非阻塞通信调用后立即返回，缓冲区不能立即使用Sum=计算某变量MPI_Isend(

5、sum.)发送该变量 sum=不能给变量重新赋值（发送可能尚未完成）MPI_Irecv(sum1,)sum=sum0+sum1 数据不能立即使用（接收可能未完成）MPI_Isend(sum,request,)Call MPI_Wait(request,status,ierr)Sum=,MPI_Irecv(sum1,request,)Call MPI_Wait(request,status,ierr)Sum=sum0+sum1,8,Copyright by Li Xinliang,利用通信与计算重叠技术提高效率,例：计算差分串行程序 real A(N,N),B(N,N),h.Do i=1,NB(

6、I,1)=(A(I,2)-A(I,1)/h B(I,N)=(A(I,N)-A(I,N-1)/henddoDo j=2,N-1Do i=1,NB(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h)EnddoEnddo,0,J=1,2,3.N-1,N,i=1i=2 i=N,9,Copyright by Li Xinliang,并行程序 以两个进程并行为例real A(N,N/2),B(N,N/2),A1(N),hIf(myid.eq.0)then call MPI_send(A(1,N/2),N,MPI_real,1,99,MPI_Comm_world,ierr)call MPI_rec

7、v(A1,N,MPI_real,1,99,MPI_Comm_World,status,ierr)Else call MPI_recv(A1,N,MPI_real,0,99,MPI_Comm_World,status,ierr)call MPI_send(A(1,1),N,MPI_real,0,99,MPI_Comm_world,ierr)endif,0,1,J=1,2 N/2,A(1,N/2)A(2,N/2)A(3,N/2)A(N,N/2),10,Copyright by Li Xinliang,If(myid.eq.0)then Do i=1,N B(i,1)=(A(i,2)-A(i,1)/

8、h B(i,N)=(A1(i)-A(i,N-1)/(2.*h)EnddoElse Do i=1,N B(i,1)=(A(i,2)-A1(i)/(2.*h)B(i,N)=(A(i,N)-A(i,N-1)/h EnddoendifDo j=2,N-1Do i=1,NB(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h)EnddoEnddo,0,1,J=1,2 N/2,特点：先收发边界信息 再进行计算缺点：通信过程中CPU 空闲,11,Copyright by Li Xinliang,“内边界”,通信与计算重叠real A(N,N/2),B(N,N/2),A1(N),hinteger my

9、id,ierr,req1,req2,status()If(myid.eq.0)then call MPI_ISend(A(1,N/2),N,MPI_real,1,99,MPI_Comm_world,req1,ierr)call MPI_Irecv(A1,N,MPI_real,1,99,MPI_Comm_World,req2,ierr)Else call MPI_Irecv(A1,N,MPI_real,0,99,MPI_Comm_World,req2,ierr)call MPI_Isend(A(1,1),N,MPI_real,0,99,MPI_Comm_world,req1,ierr)endif

10、,0,1,J=1,2 N/2,12,Copyright by Li Xinliang,Do j=2,N-1Do i=1,NB(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h)EnddoEnddoCall MPI_wait(req2,statue,ierr)If(myid.eq.0)then Do i=1,N B(I,1)=(A(I,2)-A(I,1)/h B(I,N)=(A1(i)-A(I,N-1)/(2.*h)EnddoElse Do i=1,N B(I,1)=(A(I,2)-A1(i)/(2.*h)B(I,N)=(A1(i)-A(I,N-1)/h Enddoendif,0,1,

11、J=1,2 N/2,特点：传递边界信息 同时进行计算,内点,读取系统时间 doubleprecision time time=MPI_Wtime(),13,Copyright by Li Xinliang,二、如何收发非连续数据例如：发送数组的一行A(100,50)发送 A(1,1),A(1,2),A(1,3),A(1,1),A(1,2),A(1,3),方法1.多次发送 通信开销大、效率低,14,Copyright by Li Xinliang,方法2.将发送的数据拷贝到连续的数组中dimension A(100,50),B(50)If(myid.eq.0)then Do i=1,50 B(i

12、)=A(1,i)Enddo call MPI_Send(B,50,MPI_REAL,1,99,MPI_COMM_WORLD,ierr)Else call MPI_Recv(B,50,MPI_Real,0,99,)Do i=1,50 A(1,i)=B(i)Enddoendif,不足：额外的内存占用 额外的拷贝操作通信不复杂的情况，内存拷贝工作量不大，该方法也可以采用。,效果还可以,15,Copyright by Li Xinliang,方法3：构建新的数据结构,Count:块的数量；blocklength:每块的元素个数Stride:跨度（各块起始元素之间的距离）Oldtype:旧数据类型，Ne

13、wtype:新数据类型（整数）例：integer MY_TYPE Call MPI_TYPE_VECTOR(4,1,3,MPI_REAL,MY_TYPE,ierr)Call MPI_TYPE_Commit(MY_TYPE,ierr),Stride=3,固定间隔（跨度）的非连续数据 MPI_TYPE_VECTOR(count,blocklength,stride,oldtype,newtype,ierr),A(1,1)A(1,2)A(1,3)A(1,4)A(2,1)A(2,2)A(2,3)A(2,4)A(3,1)A(3,2)A(3,3)A(3,4),4块，每块1个元素，跨度为3（个元素）,For

14、tran 数组的一行Real A(3,4).A(1,:),在内存中的排列次序,16,Copyright by Li Xinliang,例：发送三维数组中的一个面(Fortran)数组：real A(M,N,P)通信 1)A(i,:,:);2)A(:,j,:);3)A(:,:,k)通信1）A(1,1,1),A(2,1,1),A(3,1,1),A(M,1,1),A(1,2,1),A(2,2,1).,MPI_Type_Vector(N*P,1,M,MPI_Real,My_Type,ierr)通信2)A(1,1,1),A(2,1,1),A(3,1,1).,A(1,2,1),A(2,2,1),A(3,2

15、,1),A(1,1,2),A(2,1,2),A(3,1,2),MPI_Type_Vector(P,M,M*N,MPI_Real,My_Type,ierr)通信3）连续分布，无需构造新类型,17,Copyright by Li Xinliang,MPI_TYPE_INDEXED(count,array_of_blocklengths,array_of_displacements,oldtype,newtype，ierr),构造数据类型更灵活的函数 直接指定每块的元素个数及偏移量,块的数量（整数）,每块元素的个数（整形数组）,每块的偏移量（整形数组）,例：数组 real A(N,N),欲将其上三角

16、元素作为消息发送，试构造其数据类型,A(1,1),A(1,2),A(1,3),A(1,4),A(2,2),A(2,3),A(2,4),A(4,4),A(3,3),A(3,4),A(2,1),A(3,1),A(3,2),A(4,1),A(4,2),A(4,3),A(1,1),A(2,1),A(1,2),A(2,2),A(3,1),A(4,1),A(3,2),A(4,2),A(1,3),A(2,3),A(3,3),A(4,3),A(1,4),A(2,4),A(3,4),A(4,4),内存中的存储次序(Fortran),N列,N行,注意：Fortran 行优先次序存储；C为列优先次序存储,观察规律：

17、N块；第k块有k个元素；第k块的偏移为(k-1)*N（从0算起）,Integer:count,blocklengths(N),displacements(N)Integer:Newtype,ierr count=N do k=1,N blocklengthes(k)=k displacements(k)=(k-1)*N enddocall MPI_TYPE_INDEXED(count,blocklengths,&displacements,MPI_REAL,newtype，ierr)Call MPI_TYPE_Commit(Newtype,ierr)call MPI_Send(A(1,1),1

18、,Newtype,),18,Copyright by Li Xinliang,N,三、MPI的通信域和组预定义通讯域 MPI_Comm_World:包含所有进程的组通讯域的分割 MPI_Comm_Split(comm,color,key,New_Comm),0,2,1,4,3,5,7,6,8,9,10,11,Color 相同的进程在同一组根据key的大小排序（key相同时按原ID排序）例如：12个进程，分成 3行4列Integer myid,Comm_Raw,Comm_column,myid_raw,myid_line,ierr,raw,columnRaw=mod(myid,3);column

19、=int(myid/3)MPI_Comm_Split(MPI_Comm_World,raw,0,Comm_Raw)MPI_Comm_Split(MPI_Comm_World,column,0,Comm_column)Call MPI_Comm_rank(Comm_Raw,myid_raw,ierr)Call MPI_Comm_rank(Comm_line,myid_line,ierr),MPI_Comm_World,RAW,Column,Color,分组标准,Key,排序依据如相同，按原ID排,提交新定义的组（否则新组无效，不要忘记）,计算行号、列号,19,Copyright by Li Xi

20、nliang,例：计算差分 三维分割 A(M1,N1,P1)(M1=M/NM,N1=N/NN,P1=P/NP)基本思路：1）“扩大”的数组 A(0:M1+1,0:N1+1,0:P1+1)2）分割成三个组 Comm_X,Comm_Y,Comm_Z 得到组内编号 建立三个方向通讯的数据结构4）通信，计算内点差分5）计算边界差分,0,2,1,4,3,5,7,6,8,9,10,11,MPI_Comm_World,20,Copyright by Li Xinliang,Parameter(M1=M/NM,N1=N/NN,P1=P/NP)Real A(0:M1+1,0:N1+1,0:P1+1)Intege

21、r myid,Comm_X,Comm_Y,Comm_Z,id_X,id_Y,id_Z,request(12),.Call MPI_Comm_Rank(MPI_Comm_World,myid,ierr)Call MPI_Comm_Split(MPI_Comm_World,mod(myid,NM),0,Comm_X,ierr)Call MPI_Comm_Split(MPI_Comm_World,mod(myid,NM*NN)/NM,0,Comm_Y,ierr)Call MPI_Comm_Split(MPI_Comm_World,myid/(NM*NN),0,Comm_Z,ierr)Call MPI

22、_Comm_Rank(Comm_X,id_x,ierr)Call MPI_Comm_Rank(Comm_Y,id_y,ierr)Call MPI_Comm_Rank(Comm_Z,id_z,ierr),定义三个方向的通信域,21,Copyright by Li Xinliang,Call MPI_Type_Vector(N1+2)*(P1+2),1,M1+2,MPI_real,Type_X,ierr)Call MPI_Type_Vector(P1+2,N1+2,(M1+2)*(N1+2),MPI_real,Type_Y,ierr)Call MPI_Type_Commit(Type_X,ierr

23、)Call MPI_Type_Commit(Type_Y,ierr).id_X_Pre=id_X-1,if(id_X_Pre.le.0)id_X_pre=id_X_Pre+NMId_X_Next=id_X+1,if(id_X_Next.ge.NM)id_X_Next=id_X_Next-NM Call MPI_Isend(A(1,0,0),1,TYPE_X,id_X_Pre,99,Comm_X,request(1),ierr)Call MPI_Isend(A(M1,0,0),1,TYPE_X,id_X_next,99,Comm_X,request(2),ierr)Call MPI_Irecv(

24、A(0,0,0),1,TYPE_X,id_X_next,99,Comm_X,request(3),ierr)Call MPI_Irecv(A(M1+1,0,0),1,TYPE_X,id_X_Pre,99,Comm_X,request(4),ierr),定义新的数据结构,22,Copyright by Li Xinliang,Do k=2,P1-1 Do j=2,N1-1 Do i=2,M1-1 Ax(I,j,k)=(A(i+1,j,k)-A(i-1,j,k)/(2.*hx)Ay(I,j,k)=(A(I,j+1,k)-A(I,j-1,k)/(2.*hy)Az(I,j,k)=(A(I,j,k+1)

25、-A(I,j,k-1)/(2.*hz)EnddoEnddoEnddo call MPI_Wait_All(12,request,status,ierr)do k=1,P1 do j=1,N1 Ax(1,j,k)=(A(2,j,k)-A(0,j,k)/(2.*hx)Ax(M1,j,k)=(A(M1+1,j,k)-A(M1-1,j,k)/(2.*hx)enddo Enddo.,内点,边界点,23,Copyright by Li Xinliang,四、分布数组的文件存储 分布数组 real A(M/m1,N/n1)存储方式1.每个进程存储到独立的文件 real A(M/m1,N/n1)charact

26、er(len=50)filename write(filename,”(file-I4.4.dat)”)myid open(55,file=filename,form=unformatted)write(55)A close(55)-file-0000.dat file-0001.dat file-0002.dat 优点：程序简单缺点：数据文件多，不易处理；改变处理器数目时需特殊处理,0,1,2,3,24,Copyright by Li Xinliang,分布数组 real A(M/m1,N/n1)存储方式2：收集到0节点存储 存储到 一个文件 缺点：改变处理器规模时，需要处理存储方式3：收集

27、到0节点，重新装配成大数组 收集 A(M/m1,N/n1)组成 A0(M,N)real A0(M,N),A(M/m1,N/n1),A1(M/m1,N/n1)if(myid.eq.0)then do k=0,m1*n1 call MPI_recv(A1,M/m1*N/n1,MPI_real,k,.).A0(i_global,j_global)=A1(i,j)把A1 装配到A0 enddo Write(33)A0 else call MPI_Send(A,)endif,0,1,2,3,0,1,2,3,0,25,Copyright by Li Xinliang,存储方式4.按列搜集后存储,Real

28、Aj(M)If(myid.eq.0)then open(33,file=“A.dat”,form=“binary”)do j=1,N 收集矩阵A0 的第 j 列存储到 Aj(:)write(33)Aj enddoElse endif,第 1列 第 2列 第 3列,优点：存储的数据形式与内存中A0的存放格式一致。存储的文件串行程序可直接读取 real A(M,N)open(55,file=“A.dat”,form=“binary”)read(55)A close(55),26,Copyright by Li Xinliang,存储方式5 并行IO(MPI 2.0)打开文件：MPI_file_op

29、en(Comm,filename,mode,info,fileno,ierr)mode 打开类型：MPI_Mode_RDONLY,MPI_Mode_RDWR,fileno 文件号，info 整数（信息）关闭文件：MPI_file_close(fileno,ierr)指定偏移位置读写 MPI_file_read_at(fileno,offset,buff,const,datatype,status,ierr)MPI_file_write_at(fileno,offset,buff,const,datatype,status,ierr)offset 偏移，buff 缓冲区，const 数目,27,

30、Copyright by Li Xinliang,Part 3 实例教学 CFD程序的MPI实现实例（1）用拟谱方法求解不可压N-S方程 实例（2）用流水线方法计算紧致差分 常用的优化方法,28,Copyright by Li Xinliang,回顾 基本的MPI函数（6个）MPI初始化 MPI_Init(ierr)；MPI结束 MPI_Finalize(ierr)得到当前进程标识 MPI_Comm_rank(MPI_COMM_WORLD,myid,ierr)得到通信域包含的进程数 MPI_Comm_size(MPI_COMM_WORLD,numprocs,ierr)消息发送MPI_Send(

31、buf,count,datatype,dest,tag,comm,ierr)消息接收MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr),29,Copyright by Li Xinliang,非阻塞消息发送MPI_ISend(buf,count,datatype,dest,tag,comm,request,ierr)In buf,count,datatype,dest,tag,commOut request,ierr Request(返回的非阻塞通信对象,整数)非阻塞消息接收MPI_IRecv(buf,count,datatype,

32、source,tag,comm,request,ierr)In buf,count,datatype,source,tag,commOut request,ierr非阻塞通信的完成 MPI_Wait(request,status,ierr)等待消息收发完成 MPI_Test(request,flag,stutus,ierr)MPI_Waitall(const,request_array,status,ierr)等待多个消息完成 In requestOut status，flag(logical型),30,Copyright by Li Xinliang,发送非连续数据构建新的数据结构,MPI_

33、TYPE_VECTOR(count,blocklength,stride,oldtype,newtype,ierr)Count:块的数量；blocklength:每块的元素个数Stride:跨度（各块起始元素之间的距离）Oldtype:旧数据类型，Newtype:新数据类型（整数）例：integer MY_TYPE Call MPI_TYPE_VECTOR(50,1,100,MPI_REAL,MY_TYPE,ierr)Call MPI_TYPE_Commit(MY_TYPE,ierr),31,Copyright by Li Xinliang,通讯域的分割 MPI_Comm_Split(comm

34、,color,key,New_Comm),0,2,1,4,3,5,7,6,8,9,10,11,Color 相同的进程在同一组根据key的大小排序例如：12个进程，分成 3行4列Line=mod(myid,3);raw=myid/3MPI_Comm_Split(MPI_Comm_World,raw,0,Comm_Raw)MPI_Comm_Split(MPI_Comm_World,line,Comm_Line)Call MPI_Comm_rank(Comm_Raw,myid_raw,ierr)Call MPI_Comm_rank(Comm_line,myid_line,ierr),MPI_Comm

35、_World,32,Copyright by Li Xinliang,实例 1.用（拟）谱方法求解二维不可压N-S方程,2p,物理模型,周期性边界条件,按照给定能谱布置初始流动 研究流动的演化规律,33,Copyright by Li Xinliang,Fourier 变换(1D),Fourier 变换 的特点：求导数-乘积,困难：非线性项,卷积,计算量巨大,在物理空间计算,Fourier 变换的快速算法FFT,34,Copyright by Li Xinliang,二维 Fourier 变换,两次一维 Fourier 变换,35,Copyright by Li Xinliang,求解步骤：1

36、）读入初值 2）调用FFT 得到 3）计算 调用FFT 得到 4）计算 调用FFT 得到 5）计算 6）积分 求出下一时间步的值 7）调用 FFT 得到 8）循环 3）-7）直到给定的时间,36,Copyright by Li Xinliang,实际计算中，要采用抑制混淆误差的措施,程序的并行化：二维 FFT,二维FFT:调用两次一维FFT一维 FFT 算法复杂，并行化难度大二维 FFT 的并行：重新分布,Subroutine FFT2d(nx,ny,u)integer nx,ny Complex u(nx,ny),Fu(nx,ny),u1(ny),u2(nx),do i=1,nx u1(:)

37、=u(i,:)call FFT1d(ny,u1)Fu(i,:)=u1(:)enddo do j=1,ny u2(:)=Fu(:,j)call FFT1d(nx,u1)u(:,j)=u1(:)enddo end,37,Copyright by Li Xinliang,数据重分布的实现,A1(M/P,N),A2(M,N/P),1,2,3,4,a,b,c,d,对等式编程思想“我”需要完成的工作 1)将数据 A1(M/P,N)切割成P块，存入数组B1(M/P,N/P,P)2)将数据 B1(:,:,k)发到 进程 k(k=0,1.P-1)3)从进程k 接收 B2(:,:,k)4)组合 B2(:,:,k)

38、成 A2,38,Copyright by Li Xinliang,程序：Subroutine Redistibute_ItoJ(A1,A2,M,N,P)Integer M,N,P,k,ierr,status(MPI_Status_Size)real A1(M/P,N),A2(M,N/P),B1(M/P,N/P,P),B2(M/P,N/P,P)do k=1,P B1(:,:,P)=A1(:,(k-1)*N/P+1:k*N/P)call MPI_Send(B1,M*N/(P*P),MPI_Real,k-1,.)Enddodo k=1,P call MPI_Recv(B2,M*N/(P*P),MPI

39、_Real,k-1,.)A2(k-1)*M/P+1:k*M/P),:)=B2(:,:,P)Enddoend 问题：全部发送，发送成功后再启动接收。容易死锁,按行分布-按列分布,39,Copyright by Li Xinliang,Subroutine Redistibute_ItoJ(A1,A2,M,N,P)Integer M,N,P,k,ierr,status(MPI_Status_Size)real A1(M/P,N),A2(M,N/P),B1(M/P,N/P,P),B2(M/P,N/P,P)do k=1,P B1(:,:,P)=A1(:,(k-1)*N/P+1:k*N/P)id_sen

40、d=myid-k mod P id_recv=myid+k mod P call MPI_Send(B1,M*N/(P*P),MPI_Real,id_send,.)call MPI_Recv(B2,M*N/(P*P),MPI_Real,id_recv,.)A2(k-1)*M/P+1:k*M/P),:)=B2(:,:,P)Enddoend 问题：按顺序发送、接收，不易死锁,40,Copyright by Li Xinliang,数据全交换：MPI_AlltoAll(sendbuf,sendcount,sendtype,recvbuf,recvcount,recvtype,comm,ierr)se

41、ndbuf 发送缓冲区（首地址）recvbuf 接收缓冲区（首地址）sendcount 发送数目 recvcount 接收数目 sendtype 发送类型 recvtype 接收类型 Comm 通信域 ierr 整数，返回错误值（0为成功）,To 0,To 1,To 2,To 3,Sendbuf 的数据格式,sendcount,From 0,From 1,From 2,From 3,Recvbuf 的数据格式,recvcount,41,Copyright by Li Xinliang,程序：Subroutine Redistibute_ItoJ(A1,A2,M,N,P)Integer M,N,

42、P,k,ierr,status(MPI_Status_Size)real A1(M/P,N),A2(M,N/P),B1(M/P,N/P,P),B2(M/P,N/P,P)do k=1,P B1(:,:,P)=A1(:,(k-1)*N/P+1:k*N/P)enddo call MPI_AlltoAll(B1,M*N/(P*P),MPI_Real,B2,M*N/(P*P),MPI_Real,MPI_Comm_World,ierr)do k=1,P A2(k-1)*M/P+1:k*M/P),:)=B2(:,:,P)Enddoend 问题：无法做到计算与通信重叠,42,Copyright by Li X

43、inliang,二维 并行FFT 的实现（输入数据、输出数据均为按列分布）,1）调用一维FFT实现 i-方向的变换 u-u12）重新分布数据（按列-按行）u1-u2 调用一维FFT 实现j-方向的变换 u2-Fu2 重新分布数据（按行-按列）Fu2-Fu,43,Copyright by Li Xinliang,实例（2）利用流水线 实现紧致差分的并行化,紧致型差分格式：相同网格点上引入更多信息。性能更优化。,是 的差分逼近,普通差分格式：显式给出 Fi 的表达式,紧致型差分格式：隐式给出 Fi 的表达式,6 阶中心,6 阶对称紧致(Lele),5 阶迎风紧致(Fu),j-2 j-1 j j+1

44、 j+2,44,Copyright by Li Xinliang,普通差分格式：直接计算导数，并行容易,紧致格式的计算：递推,递推公式：,计算出（由边界条件或边界格式给出）2）由 递推计算 出全部导数,后面的数据必须等待前一步计算完成，无法并行,45,Copyright by Li Xinliang,二维问题：流水线法求解,流水线示意图,步骤：1)计算 d(:,:)2)for k=1,M 如果 myid=0,计算 F(k,0),否则 从myid-1接收 F(k,0);for i=1,N1(N1=N/P)计算 F(k,i);如果myid P-1 向 myid+1 发送 F(k,N1),缺点：通信

45、次数过多,46,Copyright by Li Xinliang,通信次数过于频繁解决方法：分块流水线,步骤：1)计算 d(:,:)2)for kp=1,MP 如果 myid=0,计算 F(kp,0),否则 从myid-1接收 F(kp,0);for j=1,N1(N1=N/P)计算 F(kp,j);如果myid P-1 向 myid+1 发送 F(kp,N1),F(kp,i)表示第kp 块,47,Copyright by Li Xinliang,对称紧致格式,追赶法,令,则,代入（1）得,对比（2）得,边界处导数可由边界条件或边界格式给出：,则,步骤：1）2）由（3）式递推，得到 3)4)由

46、（2）式递推，得到,特点：两次递推。并行方法与前文类似,48,Copyright by Li Xinliang,常用的并行优化方法 1）通信与计算重叠 采用非阻塞通信 Isend,Irecv 2）用重复计算代替通信 3）拆分长消息、合并短消息 4）优化通信方式,49,Copyright by Li Xinliang,用重复计算代替通信 例如：计算差分 u 分布存储,f(u)为 u 的函数,方法 1）计算出 v=f(u)通信得到 uN+1,vN+1 计算差分,方法 2）计算出 v=f(u)通信得到 uN+1（边界外）计算出 vN+1=f(uN+1)计算差分,方法 2）计算量大，通信量小 当函数

47、f(u)不复杂时，可提高效率,1,2 N N+1,50,Copyright by Li Xinliang,长消息切割成多个短消息发送、接收 call MPI_Send(A(1),100000,MPI_Real,1,)改为：do m=1,10 call MPI_Send(A(m-1)*10000+1),10000,MPI_real,1)enddo 长消息：非缓冲；短消息：缓冲,缓冲区,MPI_Send,缓冲区,MPI_Send,MPI_Recv,MPI_Recv,51,Copyright by Li Xinliang,合并短消息,do m=1,100 call MPI_Send(A(1,m),1

48、,MPI_real,1)enddo 改为 do m=1,100 B(m)=A(1,m)enddo call MPI_Send(B(1),100,MPI_Real,1,),52,Copyright by Li Xinliang,优化通信方式,例：数据散发0号 进程：数据 A(100)，散发给 0-99方式1）0 进程执行100次 MPI_Send 其他进程执行 MPI_Recv MPI_Scatter()采用该算法方式 2）0 进程 把 A(100)切割成10份，发送给10个进程 10个进程接收A1(10)后再散发,53,Copyright by Li Xinliang,参考文献：郁志辉：高性能计算之并行编程技术 MPI并行程序设计 迟学斌：高性能并行计算（),联系方式：李新亮 82543801,54,Copyright by Li Xinliang,