中科院计算流体力学最新讲义CFD2011-第14讲-MPI并行程序设计初步2.ppt

1、计算流体力学讲义第六讲 MPI并行程序设计 （2）李新亮 ；力学所主楼219； 82543801,知识点： 阻塞通信与非阻塞通信 非连续数据的发送与接收OpenMP并行程序设计初步,1,Copyright by Li Xinliang,讲义、课件上传至 （流体中文网） - “流体论坛” -“ CFD基础理论” 也可到如下网址下载：http:/cid- 程序的运行原理：服务器（前端机）编译可执行代码复制 N 份，每个节点运行一份调用MPI库函数 得到每个节点号 my_id根据my_id 不同，程序执行情况不同调用MPI 库函数进行通讯,MPI 编程的基本思想： 主从式，对等式,2,Copyr

2、ight by Li Xinliang,重点：对等式程序设计,知识回顾,Copyright by Li Xinliang,3,a.exe,对等式 设计,“对等式”程序设计思想,如果我是其中一个进程； 我应当做完成我需要完成的任务,站在其中一个进程的角度思考,基本的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)消息发送 MP

3、I_Send(buf,count,datatype,dest,tag,comm, ierr) 消息接收 MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr),4,Copyright by Li Xinliang,MPI的消息发送机制 两步进行MPI_Send( A, ) 发送 MPI_Recv( B, ) 接收,发送 变量A,接收 到变量B,配合使用,5,Copyright by Li Xinliang,一、 阻塞式通信与非阻塞式通信,阻塞式发送与接收,MPI_Send( A, ),MPI_Recv( B , ),6,Copyrigh

4、t by Li Xinliang,MPI_Send( ) 返回后缓冲区可释放sum= call MPI_Send(sum,)sum= 变量可重复利用MPI_Recv() 返回后缓冲区数据可使用Call MPI_Recv(sum1,) Sum=sum0+sum1 ,7,Copyright by Li Xinliang,非阻塞发送,启动发送,立即返回,计算,通信完成,释放发送缓冲区,发送消息,非阻塞接收,启动接收,立即返回,计算,通信完成,引用接收数据,接收消息,计算 与 通信 重叠,非阻塞消息发送与接收,8,Copyright by Li Xinliang,非阻塞消息发送 MPI_ISend(b

5、uf,count,datatype,dest,tag,comm,request,ierr) In buf,count,datatype,dest,tag,comm Out request,ierrRequest (返回的非阻塞通信对象, 整数)非阻塞消息接收 MPI_IRecv(buf,count,datatype,source,tag,comm,request,ierr) In buf,count,datatype,source,tag,comm Out request,ierr非阻塞通信的完成MPI_Wait(request,status,ierr) 等待消息收发完成MPI_Test(re

6、quest, flag,stutus,ierr)MPI_Waitall(const,request_array,status,ierr) 等待多个消息完成 In request Out status， flag (logical型),9,Copyright by Li Xinliang,非阻塞通信调用后立即返回，缓冲区不能立即使用Sum= 计算某变量 MPI_Isend(sum .) 发送该变量 sum= 不能给变量重新赋值 （发送可能尚未完成）MPI_Irecv(sum1, ) sum=sum0+sum1 数据不能立即使用 （接收可能未完成）MPI_Isend(sum, , request,

7、 ) Call MPI_Wait(request,status,ierr) Sum= ,MPI_Irecv(sum1, , request, ) Call MPI_Wait(request,status,ierr) Sum=sum0+sum1 ,10,Copyright by Li Xinliang,利用通信与计算重叠技术提高效率,例： 计算差分串行程序real A(N,N),B(N,N),h Do i=1,N B(I,1)=(A(I,2)-A(I,1)/h B(I,N)=(A(I,N)-A(I,N-1)/h enddo Do j=2,N-1 Do i=1,N B(i,j)=(A(i,j+1)

8、-A(i,j-1)/(2.*h) Enddo Enddo,0,J=1,2,3 . N-1, N,i=1 i=2i=N,11,Copyright by Li Xinliang,并行程序 以两个进程并行为例real A(N,N/2),B(N,N/2),A1(N),h If(myid .eq. 0) thencall MPI_send(A(1,N/2),N,MPI_real,1,99,MPI_Comm_world,ierr)call MPI_recv(A1,N,MPI_real,1,99,MPI_Comm_World,status,ierr) Elsecall MPI_recv(A1,N,MPI_r

9、eal,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),12,Copyright by Li Xinliang,If(myid .eq. 0) thenDo i=1,NB(i,1)=(A(i,2)-A(i,1)/h B(i,N)=(A1(i)-A(i,N-1)/(2.*h)Enddo ElseDo i=1,NB(i,1)=(A(i,2)-A1(i)/(2

10、.*h) B(i,N)=(A(i,N)-A(i,N-1)/hEnddo endif Do j=2,N-1 Do i=1,N B(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h) Enddo Enddo,0,1,J=1,2 N/2,特点：先收发边界信息再进行计算缺点：通信过程中CPU 空闲,13,Copyright by Li Xinliang,“内边界”,通信与计算重叠real A(N,N/2),B(N,N/2),A1(N),h integer myid,ierr, req1, req2,status() If(myid .eq. 0) thencall MPI_ISend(A(1

11、,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) Elsecall 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,0,1,J=1,2 N/2,14,Copyright by Li Xinliang,Do j=2,N-1 Do i=

12、1,N B(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h) Enddo EnddoCall MPI_wait(req2,statue,ierr)If(myid .eq. 0) thenDo i=1,NB(I,1)=(A(I,2)-A(I,1)/h B(I,N)=(A1(i)-A(I,N-1)/(2.*h)Enddo ElseDo i=1,NB(I,1)=(A(I,2)-A1(i)/(2.*h) B(I,N)=(A1(i)-A(I,N-1)/hEnddo endif,0,1,J=1,2 N/2,特点：传递边界信息同时进行计算,内点,读取系统时间doubleprecision tim

13、etime=MPI_Wtime( ),15,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. 多次发送通信开销大、效率低,16,Copyright by Li Xinliang,方法2. 将发送的数据拷贝到连续的数组中dimension A(100,50), B(50) If(myid .eq. 0) thenDo i=1,50B(i)=A(1,i)Enddocall MPI_Send(B,50,MPI_REAL,1,9

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

15、 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（个元素）,Fortran 数组的一行 Real A

16、(3,4) A(1,:),在内存中的排列次序,18,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,

17、2,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）连续分布，无需构造新类型,19,Copyright by Li Xinliang,MPI_TYPE_INDEXED(count, array_of_blocklengths, array_of_displacements, oldtype,newtype，ierr),构造数据类型更灵活的函数 直接指定每块的元素个数及偏移量,块的数量（整数）,每块元素的个数 （整形数组）,每块的偏移量 （整形数组）,例： 数组 real A(

18、N,N), 欲将其上三角元素作为消息发送，试构造其数据类型,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 行优先次序存储；

19、C为列优先次序存储,观察规律： N块； 第k块有k个元素；第k块的偏移为(k-1)*N （从0算起）,Integer: count, blocklengths(N), displacements(N) Integer: Newtype,ierrcount=Ndo k=1,Nblocklengthes(k)=kdisplacements(k)=(k-1)*N enddo call MPI_TYPE_INDEXED(count, blocklengths, & displacements,MPI_REAL,newtype，ierr) Call MPI_TYPE_Commit(Newtype, ie

20、rr)call MPI_Send (A(1,1),1,Newtype, ),20,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

21、,ierr,raw,column Raw=mod(myid,3); column=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排,

22、提交新定义的组 （否则新组无效，不要忘记）,计算行号、列号,21,Copyright by Li Xinliang,例： 计算差分 三维分割 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,22,Copyright by Li Xinliang,Paramet

23、er(M1=M/NM,N1=N/NN,P1=P/NP) Real A(0:M1+1,0:N1+1,0:P1+1) Integer 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) Ca

24、ll MPI_Comm_Split(MPI_Comm_World,myid/(NM*NN),0,Comm_Z,ierr)Call MPI_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),定义三个方向的通信域,23,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,N

25、1+2,(M1+2)*(N1+2),MPI_real,Type_Y,ierr) Call MPI_Type_Commit(Type_X,ierr) 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+NM Id_X_Next=id_X+1, if(id_X_Next .ge. NM) id_X_Next=id_X_Next-NMCall MPI_Isend(A(1,0,0) ,1,TYPE_X, id_X_Pre, 99,Comm_X,request(1),ierr)Ca

26、ll MPI_Isend(A(M1,0,0),1,TYPE_X,id_X_next,99,Comm_X,request(2),ierr)Call MPI_Irecv(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) ,定义新的数据结构,24,Copyright by Li Xinliang,Do k=2,P1-1Do j=2,N1-1Do i=2,M1-1Ax(I,j,k)=(A(i+1,j,

27、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)-A(I,j,k-1)/(2.*hz) Enddo Enddo Enddocall MPI_Wait_All(12,request,status,ierr)do k=1,P1do j=1,N1Ax(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)enddoEnddo.,内点,边界点,25,Copyright by Li Xinl

28、iang,四、分布数组的文件存储分布数组 real A(M/m1,N/n1)存储方式1. 每个进程存储到独立的文件real A(M/m1,N/n1)character(len=50) filenamewrite(filename,”(file-I4.4.dat)”) myidopen(55,file=filename,form=unformatted)write(55) Aclose(55) -file-0000.dat file-0001.dat file-0002.dat 优点：程序简单 缺点： 数据文件多，不易处理； 改变处理器数目时需特殊处理,0,1,2,3,26,Copyright

29、by Li Xinliang,分布数组 real A(M/m1,N/n1)存储方式2： 收集到0节点存储存储到 一个文件缺点： 改变处理器规模时，需要处理存储方式3：收集到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) thendo k=0,m1*n1call MPI_recv(A1, M/m1*N/n1,MPI_real,k,) A0( i_global, j_global ) = A1(i,j ) 把A1 装配到A0enddoWrite(33) A0

30、elsecall MPI_Send(A,)endif,0,1,2,3,0,1,2,3,0,27,Copyright by Li Xinliang,存储方式4. 按列搜集后存储,Real Aj(M) If( myid .eq. 0) thenopen(33,file=“A.dat”,form= “binary”)do j=1,N收集矩阵A0 的第 j 列存储到 Aj(:)write(33) Ajenddo Else endif,第 1列 第 2列 第 3列,优点： 存储的数据形式与内存中A0的存放格式一致。 存储的文件串行程序可直接读取real A(M,N)open(55,file=“A.dat

31、”,form=“binary”)read(55) Aclose(55),28,Copyright by Li Xinliang,存储方式5 并行IO (MPI 2.0)打开文件：MPI_file_open(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,s

32、tatus,ierr)MPI_file_write_at(fileno,offset,buff,const,datatype,status,ierr)offset 偏移， buff 缓冲区，const 数目,29,Copyright by Li Xinliang,Part 3 实例教学 CFD程序的MPI实现实例 （1） 用拟谱方法求解不可压N-S方程实例（2） 用流水线方法计算紧致差分常用的优化方法,30,Copyright by Li Xinliang,回顾 基本的MPI函数（6个） MPI初始化 MPI_Init(ierr) ； MPI结束 MPI_Finalize(ierr) 得到当前

33、进程标识 MPI_Comm_rank(MPI_COMM_WORLD,myid,ierr)得到通信域包含的进程数MPI_Comm_size(MPI_COMM_WORLD,numprocs,ierr)消息发送 MPI_Send(buf,count,datatype,dest,tag,comm, ierr) 消息接收 MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr),31,Copyright by Li Xinliang,非阻塞消息发送 MPI_ISend(buf,count,datatype,dest,tag,comm,request

34、,ierr) In buf,count,datatype,dest,tag,comm Out request,ierrRequest (返回的非阻塞通信对象, 整数)非阻塞消息接收 MPI_IRecv(buf,count,datatype,source,tag,comm,request,ierr) In buf,count,datatype,source,tag,comm Out request,ierr非阻塞通信的完成MPI_Wait(request,status,ierr) 等待消息收发完成MPI_Test(request, flag,stutus,ierr)MPI_Waitall(con

35、st,request_array,status,ierr) 等待多个消息完成 In request Out status， flag (logical型),32,Copyright by Li Xinliang,发送非连续数据构建新的数据结构,MPI_TYPE_VECTOR(count,blocklength,stride,oldtype,newtype,ierr)Count: 块的数量； blocklength: 每块的元素个数 Stride: 跨度 （各块起始元素之间的距离） Oldtype: 旧数据类型， Newtype: 新数据类型 （整数）例：integer MY_TYPE Call

36、 MPI_TYPE_VECTOR(50,1,100,MPI_REAL,MY_TYPE,ierr)Call MPI_TYPE_Commit(MY_TYPE,ierr),33,Copyright by Li Xinliang,通讯域的分割MPI_Comm_Split(comm,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/3 MPI_Comm_Split(MPI_Comm_World, raw, 0,Com

37、m_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_World,34,Copyright by Li Xinliang,实例 1. 用（拟）谱方法求解二维不可压N-S方程,2p,物理模型,周期性边界条件,按照给定能谱布置初始流动研究流动的演化规律,35,Copyright by Li Xinliang,Fourier 变换 (1D),Fourier

38、变换 的特点： 求导数 - 乘积,困难： 非线性项,卷积,计算量巨大,在物理空间计算,Fourier 变换的快速算法FFT,36,Copyright by Li Xinliang,二维 Fourier 变换,两次一维 Fourier 变换,37,Copyright by Li Xinliang,求解步骤：1） 读入初值2） 调用FFT 得到3） 计算 调用FFT 得到 4） 计算 调用FFT 得到 5） 计算6） 积分 求出下一时间步的值7） 调用 FFT 得到8） 循环 3）-7） 直到给定的时间,38,Copyright by Li Xinliang,实际计算中，要采用抑制混淆误差的措施,

39、程序的并行化： 二维 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,nxu1(:)= u(i,:)call FFT1d(ny,u1)Fu(i,:)=u1(:) enddodo j=1,nyu2(:)=Fu(:,j)call FFT1d(nx,u1)u(:,j)=u1(:) enddoend,39,Copyright by Li Xinliang,数据重分布

40、的实现,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) 成 A2,40,Copyright by Li Xinliang,程序：Subroutine Redistibute_ItoJ(A1,A2,M,N,P) Integer M,N,P,k,ierr,status(MPI_Status_Size) real

41、 A1(M/P,N), A2(M,N/P), B1(M/P,N/P,P), B2(M/P,N/P,P)do k=1,PB1(:,:,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,Pcall MPI_Recv(B2,M*N/(P*P),MPI_Real, k-1, .) A2(k-1)*M/P+1: k*M/P) , : )=B2(:,:,P) Enddo end问题： 全部发送， 发送成功后再启动接收。 容易死锁,按行分布 - 按列分布,41,Copyright b

42、y 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,PB1(:,:,P)=A1(:, (k-1)*N/P+1: k*N/P) ) id_send=myid-k mod Pid_recv= myid+k mod Pcall MPI_Send(B1,M*N/(P*P),MPI_Real, id_send, .) call MPI_

43、Recv(B2,M*N/(P*P),MPI_Real, id_recv, .) A2(k-1)*M/P+1: k*M/P) , : )=B2(:,:,P) Enddoend问题： 按顺序发送、接收，不易死锁,42,Copyright by Li Xinliang,数据全交换： MPI_AlltoAll(sendbuf,sendcount,sendtype,recvbuf,recvcount,recvtype,comm,ierr)sendbuf 发送缓冲区（首地址） recvbuf 接收缓冲区（首地址）sendcount 发送数目 recvcount 接收数目sendtype 发送类型 recv

44、type 接收类型Comm 通信域 ierr 整数，返回错误值（0为成功）,To 0,To 1,To 2,To 3,Sendbuf 的数据格式,sendcount,From 0,From 1,From 2,From 3,Recvbuf 的数据格式,recvcount,43,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,

45、N/P,P)do k=1,PB1(:,:,P)=A1(:, (k-1)*N/P+1: k*N/P) ) enddocall MPI_AlltoAll (B1,M*N/(P*P),MPI_Real, B2, M*N/(P*P),MPI_Real,MPI_Comm_World,ierr)do k=1,PA2(k-1)*M/P+1: k*M/P) , : )=B2(:,:,P) Enddo end问题： 无法做到计算与通信重叠,44,Copyright by Li Xinliang,二维 并行FFT 的实现 （输入数据、输出数据均为按列分布）,1） 调用一维FFT实现 i- 方向的变换 u - u1

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

47、接计算导数，并行容易,紧致格式的计算： 递推,递推公式：,计算出 （由边界条件或边界格式给出）2） 由 递推计算 出全部导数,后面的数据必须等待前一步计算完成，无法并行,47,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),缺点： 通信次数过多,48,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),