1、计算流体力学讲义计算流体力学讲义2020 第第15讲讲 MPI并行程序设计并行程序设计（2）李新亮李新亮知识点：知识点：阻塞通信与非阻塞通信阻塞通信与非阻塞通信 非连续数据的发送与接收非连续数据的发送与接收 OpenMPOpenMP并行程序设计初步并行程序设计初步 1Copyright by Li Xinliang服务器服务器/前端机前端机计算节点计算节点a.exea.exea.exeMPI 程序的运行原理：程序的运行原理：服务器（前端机）编译服务器（前端机）编译 可执行代码可执行代码复制复制 N 份，份，每个节点运行一份每个节点运行一份 调用调用MPI库函数库函数 得到每个节点号得到每个节点
2、号 my_id 根据根据my_id 不同，程序执行情况不同不同，程序执行情况不同 调用调用MPI 库函数进行通讯库函数进行通讯MPI 编程的基本思想：编程的基本思想：主从式，主从式，对等式对等式2Copyright by Li Xinliang重点：对等式程序设计重点：对等式程序设计知识回顾知识回顾Copyright by Li Xinliang3推荐学习网站：推荐学习网站：“超算习堂超算习堂”（）easyhpc.Net 教材资料、在线课程；在线实训Copyright by Li Xinliang4计算节点计算节点a.exea.exea.exea.exe对等式对等式设计设计“对等式对等式”程序
3、设计思想程序设计思想如果我是其中一个进程；如果我是其中一个进程；我应当做我应当做完成我需要完成的任务完成我需要完成的任务站在其中一个进程的角度思考站在其中一个进程的角度思考基本的基本的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(buf,
4、count,datatype,dest,tag,comm,ierr)消息接收消息接收MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr)5Copyright by Li XinliangMPI的消息发送机制的消息发送机制 两步进行两步进行MPI_Send(A,)发送发送MPI_Recv(B,)接收接收 发送发送 变量变量A接收接收 到变量到变量B配合使用配合使用6Copyright by Li Xinliang阻塞发送阻塞发送开始开始结束结束消息成功发出消息成功发出缓冲区可释放缓冲区可释放阻塞接收阻塞接收开始开始结束结束消息成功接收消
5、息成功接收缓冲区数据可使用缓冲区数据可使用一、一、阻塞式通信与非阻塞式通信阻塞式通信与非阻塞式通信阻塞式发送与接收阻塞式发送与接收MPI_Send(A,)MPI_Recv(B,)7Copyright by Li Xinliang MPI_SendMPI_Send()()返回后缓冲区可释放返回后缓冲区可释放 sum=sum=call call MPI_SendMPI_Send(sum,)(sum,)sum=sum=变量可重复利用变量可重复利用 MPI_RecvMPI_Recv()()返回后缓冲区数据可使用返回后缓冲区数据可使用Call Call MPI_RecvMPI_Recv(sum1,)(s
6、um1,)Sum=sum0+sum1Sum=sum0+sum1 8Copyright by Li Xinliang非阻塞发送非阻塞发送启动发送启动发送立即返回立即返回计计算算通信完成通信完成释放发送缓冲区释放发送缓冲区发发送送消消息息非阻塞接收非阻塞接收启动接收启动接收立即返回立即返回计计算算通信完成通信完成引用接收数据引用接收数据接接收收消消息息计算计算与与通信通信重叠重叠非阻塞消息发送与接收非阻塞消息发送与接收9Copyright by Li Xinliang非阻塞消息发送非阻塞消息发送MPI_ISend(buf,count,datatype,dest,tag,comm,request,i
7、err)In buf,count,datatype,dest,tag,commOut request,ierr Request (返回的非阻塞通信对象返回的非阻塞通信对象,整数整数)非阻塞消息接收非阻塞消息接收MPI_IRecv(buf,count,datatype,source,tag,comm,request,ierr)In buf,count,datatype,source,tag,commOut request,ierr非阻塞通信的完成非阻塞通信的完成 MPI_Wait(request,status,ierr)等待消息收发完成等待消息收发完成 MPI_Test(request,flag
8、,stutus,ierr)MPI_Waitall(const,request_array,status,ierr)等待多个消息完成等待多个消息完成 In requestOut status，flag (logical型型)10Copyright by Li Xinliang非阻塞通信调用后立即返回，缓冲区不能非阻塞通信调用后立即返回，缓冲区不能立即使用立即使用Sum=计算某变量计算某变量MPI_Isend(sum.)发送该变量发送该变量 sum=不能给变量重新赋值不能给变量重新赋值 （发送可能尚未完成）（发送可能尚未完成）MPI_Irecv(sum1,)sum=sum0+sum1 数据不能立即
9、使用数据不能立即使用（接收可能未完成）（接收可能未完成）MPI_Isend(sum,request,)Call MPI_Wait(request,status,ierr)Sum=MPI_Irecv(sum1,request,)Call MPI_Wait(request,status,ierr)Sum=sum0+sum1 11Copyright by Li Xinliang利用通信与计算重叠技术提高效率利用通信与计算重叠技术提高效率例：例：计算差分计算差分串行程序串行程序 real A(N,N),B(N,N),h.Do i=1,NB(I,1)=(A(I,2)-A(I,1)/h B(I,N)=(A
10、(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)EnddoEnddo0)2/()()(1,1,hfffjijijiyJ=1,2,3.N-1,Ni=1i=2 i=N12Copyright 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
11、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)endif01J=1,2 N/2A(1,N/2)A(2,N/2)A(3,N/2)A(N,N/2)13Copyright by Li XinliangIf(myid.eq.0)then Do i=1,N B(i,1)=(A(i,2)-A(i,1)/h B(i,N
12、)=(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)EnddoEnddo01J=1,2 N/2特点：特点：先收发边界信息先收发边界信息 再进行计算再进行计算缺点：缺点：通信过程中通信过程中CPU 空闲空闲14Copyright by Li Xinliang“内边界”通信与计算重叠通信与计算重叠real A(N,N/2),B(N,N/2),
13、A1(N),hinteger myid,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_worl
14、d,req1,ierr)endif01J=1,2 N/215Copyright by Li XinliangDo 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 Endd
15、oendif01J=1,2 N/2特点：特点：传递边界信息传递边界信息 同时同时进行计算进行计算内点内点读取系统时间读取系统时间 doubleprecision time time=MPI_Wtime()16Copyright 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.多次发送多次发送 通信开销大、效率低通信开销大、效率低A(1,1),A(2,1),A(1,2),A(2,2).A(1,3).17Copyri
16、ght by Li Xinliang方法方法2.将发送的数据拷贝到连续的数组中将发送的数据拷贝到连续的数组中dimension A(100,50),B(50)If(myid.eq.0)then Do i=1,50 B(i)=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不足：不足：额外的内存占用额外的内存占用 额外的拷贝操作额外的拷贝操作通信不复杂的情况，内存拷贝工作量不大
17、，该方法也可以采用。通信不复杂的情况，内存拷贝工作量不大，该方法也可以采用。效果还效果还可以可以18Copyright by Li Xinliang方法方法3：构建新的数据结构构建新的数据结构 Count:块的数量；块的数量；blocklength:每块的元素个数每块的元素个数Stride:跨度跨度（各块起始元素之间的距离）（各块起始元素之间的距离）Oldtype:旧数据类型，旧数据类型，Newtype:新数据类型新数据类型 （整数）（整数）例：例：integer MY_TYPE Call MPI_TYPE_VECTOR(4,1,3,MPI_REAL,MY_TYPE,ierr)Call MPI
18、_TYPE_Commit(MY_TYPE,ierr)A(1,1),A(2,1),A(3,1),A(1,2),A(2,2),A(3,2),A(1,3),A(2,3),A(3,3),A(1,4),A(2,4),A(3,4)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个元素，跨度为个元素
19、，跨度为3（个元素）（个元素）Fortran 数组的一行数组的一行Real A(3,4).A(1,:)在内存中的排列次序19Copyright 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),
20、A(2,1,1),A(3,1,1).,A(1,2,1),A(2,2,1),A(3,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）连续分布，无需构造新类型连续分布，无需构造新类型 20Copyright by Li XinliangMPI_TYPE_INDEXED(count,array_of_blocklengths,array_of_displacements,oldtype,newtype，ierr)构造数据类型更灵活的函数构造数据类型更灵活的函数 直接指定每块的元素个
21、数及偏移量直接指定每块的元素个数及偏移量块的数量块的数量（整数）（整数）每块元素的个数每块元素的个数（整形数组）（整形数组）每块的偏移量每块的偏移量（整形数组）（整形数组）例：例：数组数组 real A(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(
22、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为列优先次序存储为列优先次序存储观察规律：观察规律：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-
23、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,Newtype,)21Copyright by Li XinliangN三、三、MPI的通信域和组的通信域和组1.预定义通讯域预定义通讯域 MPI_Comm_World:包含所有进程的组包含所有进程的组2.通讯域的分割通讯域的分割 MPI_Comm_Split(comm,color,key,New_Comm)
24、02143576891011Color 相同的进程在同一组相同的进程在同一组根据根据key的大小排序的大小排序（key相同时按原相同时按原ID排序）排序）例如：例如：12个进程，个进程，分成分成 3行行4列列Integer myid,Comm_Raw,Comm_column,myid_raw,myid_line,ierr,raw,columnRaw=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_colum
25、n)Call MPI_Comm_rank(Comm_Raw,myid_raw,ierr)Call MPI_Comm_rank(Comm_line,myid_line,ierr)MPI_Comm_WorldRAWColumnColor,分组标准Key,排序依据如相同，按原ID排 提交新定义的组提交新定义的组（否则新组无效，不要忘记）（否则新组无效，不要忘记）计算行号、列号22Copyright by Li Xinliang例：例：计算差分计算差分 三维分割三维分割 A(M1,N1,P1)(M1=M/NM,N1=N/NN,P1=P/NP)基本思路：基本思路：1）“扩大扩大”的数组的数组 A(0:M
26、1+1,0:N1+1,0:P1+1)2）分割成三个组）分割成三个组 Comm_X,Comm_Y,Comm_Z 得到组内编号得到组内编号 3)建立三个方向通讯的数据结构建立三个方向通讯的数据结构4）通信通信 ，计算内点差分计算内点差分5）计算边界差分计算边界差分zfyfxf,02143576891011MPI_Comm_World23Copyright by Li XinliangParameter(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_
27、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_Comm_Rank(Comm_X,id_x,ierr)Call MPI_Comm
28、_Rank(Comm_Y,id_y,ierr)Call MPI_Comm_Rank(Comm_Z,id_z,ierr)定义三个方向的通信域定义三个方向的通信域24Copyright by Li XinliangCall 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)Call MPI_Type_Commit(Type_Y,ierr).
29、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(A(0,0,0),1,TYPE_X,id_X_next,99,Comm
30、_X,request(3),ierr)Call MPI_Irecv(A(M1+1,0,0),1,TYPE_X,id_X_Pre,99,Comm_X,request(4),ierr)定义新的数据结构定义新的数据结构25Copyright 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)-A(I,j,k-1)/(2.*hz)EnddoEnddoE
31、nddo 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 .内点内点边界点边界点26Copyright by Li Xinliang四、分布数组的文件存储四、分布数组的文件存储 分布数组分布数组 real A(M/m1,N/n1)real A(M/m1,N/n1)存储方式存储方式1.每个进程存储到独立的文件每个进程存储到独立的文件 real
32、A(M/m1,N/n1)character(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 优点：程序简单优点：程序简单缺点：缺点：数据文件多，不易处理；数据文件多，不易处理；改变处理器数目时需特殊处理改变处理器数目时需特殊处理 012327Copyright by Li Xinliang 分布数组分布数组 real A(M/m1,N
33、/n1)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)then do k=0,m1*n1 call MPI_recv(A1,M/m1*N/n1,MPI_real,k,.).A0(i_global,j_global )=A1
34、(i,j)把把A1 装配到装配到A0 enddo Write(33)A0 else call MPI_Send(A,)endif 01230123028Copyright by Li Xinliang存储方式存储方式4.按列搜集后存储按列搜集后存储 Real 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列列优点：优点：存储的数据形式与内存中存储的数据形式
35、与内存中A0的存放格式一致。的存放格式一致。存储的文件串行程序可直接读取存储的文件串行程序可直接读取 real A(M,N)open(55,file=“A.dat”,form=“binary”)read(55)A close(55)29Copyright 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 整数整数（信息
36、）（信息）关闭文件关闭文件：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 数目数目 30Copyright by Li XinliangPart 3 实例教学实例教学 CFD程序的程序的MPI实现实现 实例实例（1）用拟谱方法求解不可压用拟
37、谱方法求解不可压N-S方程方程 实例（实例（2）用流水线方法计算紧致差分用流水线方法计算紧致差分 常用的优化方法常用的优化方法31Copyright 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_S
38、end(buf,count,datatype,dest,tag,comm,ierr)消息接收消息接收MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr)32Copyright by Li Xinliang非阻塞消息发送非阻塞消息发送MPI_ISend(buf,count,datatype,dest,tag,comm,request,ierr)In buf,count,datatype,dest,tag,commOut request,ierr Request (返回的非阻塞通信对象返回的非阻塞通信对象,整数整数)非阻塞消息接收非阻塞消
39、息接收MPI_IRecv(buf,count,datatype,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 (l
40、ogical型型)33Copyright by Li Xinliang发送非连续数据发送非连续数据构建新的数据结构构建新的数据结构MPI_TYPE_VECTOR(count,blocklength,stride,oldtype,newtype,ierr)Count:块的数量；块的数量；blocklength:每块的元素个数每块的元素个数Stride:跨度跨度（各块起始元素之间的距离）（各块起始元素之间的距离）Oldtype:旧数据类型，旧数据类型，Newtype:新数据类型新数据类型 （整数）（整数）例：例：integer MY_TYPE Call MPI_TYPE_VECTOR(50,1,1
41、00,MPI_REAL,MY_TYPE,ierr)Call MPI_TYPE_Commit(MY_TYPE,ierr)A(1,1),A(2,1),A(1,2),A(2,2).A(1,3).34Copyright by Li Xinliang通讯域的分割通讯域的分割 MPI_Comm_Split(comm,color,key,New_Comm)02143576891011Color 相同的进程在同一组相同的进程在同一组根据根据key的大小排序的大小排序例如：例如：12个进程，个进程，分成分成 3行行4列列Line=mod(myid,3);raw=myid/3MPI_Comm_Split(MPI_
42、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_World35Copyright by Li Xinliang实例实例 1.用（拟）谱方法求解二维不可压用（拟）谱方法求解二维不可压N-S方程方程jjijjijiiixxuxpxuutuxuRe102p物理模型物理模型周期性边界条件周期性边界条件按照给定能谱布置初始流
43、动按照给定能谱布置初始流动 研究流动的演化规律研究流动的演化规律36Copyright by Li XinliangFourier 变换变换 (1D)NjxekuNjuuFujNNkkxij/2,)(1)(:)(12/2/p101)()(:)(NjkxijejukuuFuFourier 变换变换 的特点：的特点：求导数求导数-乘积乘积)/(1dxduFuik0iixu02211ukuk困难：困难：非线性项非线性项)(jijxuuF)(*)(jijxuFuF卷积卷积计算量巨大计算量巨大在物理空间计算在物理空间计算jijxuuFourier 变换的快速算法变换的快速算法FFT37Copyright
44、 by Li Xinliang二维 Fourier 变换变换12/2/12/2/22),(1),(:)(NNkmyikxiNNmjieemkuNjiuuFu)()(,2jixymkuFFuF两次一维两次一维 Fourier 变换变换jjijjijiiixxuxpxuutuxuRe10)2(Re1)1(02iiiiiukpktuuk)(jijxuuF0)2(iiiukconsiderk2/kkpiiiiukkktuRe1)/1(2238Copyright by Li XinliangiiiukkktuRe1)/1(22)(jijxuuF求解步骤：求解步骤：1）读入初值读入初值 2）调用调用FFT
45、 得到得到 3）计算计算 调用调用FFT 得到得到 4）计算计算 调用调用FFT 得到得到 5）计算计算 6）积分积分 求出下一时间步的值求出下一时间步的值 7）调用调用 FFT 得到得到 8）循环循环 3）-7）直到给定的时间直到给定的时间 ),(vuuiiu ijukijixujijxuu)(jijxuuFiiiukkkQRe1)/1(22iiQtu)1(niu)1(niu39Copyright by Li Xinliang实际计算中，要采用抑制混淆误差的措施程序的并行化：程序的并行化：二维二维 FFT)()(,2jixymkuFFuF二维二维FFT:调用两次一维调用两次一维FFT一维一维
46、 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(:)=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 end40Copyright by Li Xinliang数据重分布的实现数据重分布的实现A1(M/P
47、,N)A2(M,N/P)1234abcd对等式编程思想对等式编程思想“我我”需要完成的工作需要完成的工作 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)成成 A241Copyright by Li Xinliang程序：程序：Subroutine Redistibute_ItoJ(A1,A2,M,N,P)Integer M,N,P,k,ierr,status(MPI_Status_
48、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_Real,k-1,.)A2(k-1)*M/P+1:k*M/P),:)=B2(:,:,P)Enddoend 问题：问题：全部发送，全部发送，发送成功发送成功后再启动接收。后再启动接收。容易死锁容易死锁 按行分布按行分布 -按列分布
49、按列分布42Copyright by Li XinliangSubroutine 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_send=myid-k mod P id_recv=myid+k mod P call MPI_Send(B1,M*N/(P*P),MPI_Real,id_send,.)call
50、 MPI_Recv(B2,M*N/(P*P),MPI_Real,id_recv,.)A2(k-1)*M/P+1:k*M/P),:)=B2(:,:,P)Enddoend 问题：问题：按顺序发送、接收，不易死锁按顺序发送、接收，不易死锁43Copyright by Li Xinliang数据全交换：数据全交换：MPI_AlltoAll(sendbuf,sendcount,sendtype,recvbuf,recvcount,recvtype,comm,ierr)sendbuf 发送缓冲区（首地址）发送缓冲区（首地址）recvbuf 接收缓冲区（首地址）接收缓冲区（首地址）sendcount 发送数
