教学课件第6讲MPI并行程序设计初步

《教学课件第6讲MPI并行程序设计初步》由会员分享，可在线阅读，更多相关《教学课件第6讲MPI并行程序设计初步（54页珍藏版）》请在金锄头文库上搜索。

1、计算流体力学讲义计算流体力学讲义 第六讲第六讲 MPI并行程序设计并行程序设计 （2）李新亮李新亮 ；力学所主楼；力学所主楼219； 82543801 知识点：知识点： 阻塞通信与非阻塞通信阻塞通信与非阻塞通信 非连续数据的发送与接收非连续数据的发送与接收 1讲义、课件上传至讲义、课件上传至 (流体中文网）流体中文网） - “流体论坛流体论坛” -“ CFD基础理论基础理论 ”Copyright by Li XinliangPart 2 MPI 并行程序设计并行程序设计 1. 阻塞式通信与非阻塞式通信阻塞式通信与非阻塞式通信 2. 非连续数据的传送与接收非连续数据的传送与接收 3. MPI的

2、通信域与组；的通信域与组； 算例计算三维差分算例计算三维差分 4. 分布式数据的存储分布式数据的存储Part 3 实例教学实例教学 1. 谱方法求解不可压缩谱方法求解不可压缩N-S方程方程 2. 利用紧致型差分格式求导数利用紧致型差分格式求导数 2Copyright by Li XinliangMPI的消息的消息发送机制送机制 两步两步进行行MPI_Send( A, ) 发送送MPI_Recv( B, ) 接收接收 发送发送 变量变量A接收接收 到变量到变量B配合使用配合使用3Copyright by Li Xinliang阻塞发送阻塞发送开始开始结束结束消息成功发出消息成功发出缓冲区可释放缓

3、冲区可释放阻塞接收阻塞接收开始开始结束结束消息成功接收消息成功接收缓冲区数据可使用缓冲区数据可使用一、一、 阻塞式通信与非阻塞式通信阻塞式通信与非阻塞式通信阻塞式发送与接收阻塞式发送与接收MPI_Send( A, )MPI_Recv( B , )4Copyright by Li Xinliang MPI_Send( ) MPI_Send( ) 返回后缓冲区可释放返回后缓冲区可释放 sum= sum= call MPI_Send(sum,) call MPI_Send(sum,) sum= sum= 变量可重复利用变量可重复利用 MPI_Recv() MPI_Recv() 返回后缓冲区数据可使用

4、返回后缓冲区数据可使用Call MPI_Recv(sum1,)Call MPI_Recv(sum1,)Sum=sum0+sum1Sum=sum0+sum1 5Copyright by Li Xinliang非阻塞发送非阻塞发送启动发送启动发送立即返回立即返回计计算算通信完成通信完成释放发送缓冲区释放发送缓冲区发发送送消消息息非阻塞接收非阻塞接收启动接收启动接收立即返回立即返回计计算算通信完成通信完成引用接收数据引用接收数据接接收收消消息息计算计算与与通信通信重叠重叠非阻塞消息发送与接收非阻塞消息发送与接收6Copyright by Li Xinliang非阻塞消息非阻塞消息发送送MPI_ISe

5、nd(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,request,ierr)In buf,count,datatype,source,tag,commOut request,ierr非阻塞通信的完成非阻塞通信的完成 MPI_Wait(request,statu

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

7、能尚未完成）（发送可能尚未完成）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 8Copyright by Li Xinliang利用通信与计算重叠技术提高效率利用通信与计算重叠技术提高效率例：例： 计算差分计算差分串行程序串行程序

8、 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(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)EnddoEnddo0J=1,2,3 . N-1, Ni=1i=2 i=N9Copyright 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,M

9、PI_real,1,99,MPI_Comm_world,ierr) call MPI_recv(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)10Copyright by Li XinliangIf(myid .

10、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)=(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 空闲空闲11Copyright by

11、 Li Xinliang“内边界”通信与计算重叠通信与计算重叠real A(N,N/2),B(N,N/2),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,

12、req2,ierr) call MPI_Isend(A(1,1),N,MPI_real,0,99,MPI_Comm_world,req1,ierr)endif01J=1,2 N/212Copyright 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

13、 Do i=1,N B(I,1)=(A(I,2)-A1(i)/(2.*h) B(I,N)=(A1(i)-A(I,N-1)/h Enddoendif01J=1,2 N/2特点：特点： 传递边界信息传递边界信息 同时同时进行计算进行计算内点内点读取系统时间读取系统时间 doubleprecision time time=MPI_Wtime( ) 13Copyright by Li Xinliang二、二、 如何收发非连续数据如何收发非连续数据例如：例如： 发送数组的一行发送数组的一行A(100,50)发送发送 A(1,1),A(1,2) ,A(1,3)A(1,1), A(1,2), A(1,3)

14、方法方法1. 多次发送多次发送 通信开销大、效率低通信开销大、效率低A(1,1), A(2,1), A(1,2), A(2,2) . A(1,3).14Copyright 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,

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

16、ype: 新数据类型新数据类型 （整数）（整数）例：例：integer MY_TYPE Call MPI_TYPE_VECTOR(4,1,3,MPI_REAL,MY_TYPE,ierr) Call MPI_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, str

17、ide ,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(3,4). A(1,:)在内存中的排列次序16Copyright by Li Xinliang例：例： 发送三维数组中的一个面发送三维数组中的一个面 (Fortran) 数组：数组： real A(M,N,P) 通信通信 1) A(i,:,:) ; 2) A(

18、:,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,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、无需构造新类型 17Copyright by Li XinliangMPI_TYPE_INDEXED(count, array_of_blocklengths, array_of_displacements, oldtype,newtype，ierr)构造数据类型更灵活的函数构造数据类型更灵活的函数 直接指定每块的元素个数及偏移量直接指定每块的元素个数及偏移量块的数量块的数量（整数）（整数）每块元素的个数每块元素的个数（整形数组）（整形数组）每块的偏移量每块的偏移量（整形数组）（整形数组）例：例： 数组数组 real A(N,N), 欲将其欲将其上三角元素上三角元素作为消息发送，试构造其数据类

20、型作为消息发送，试构造其数据类型 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为列优先次序存储为列优先次序存储观察规律：观察规律： N块；块； 第第k块

21、有块有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) ca

22、ll MPI_Send (A(1,1),1,Newtype, )18Copyright by Li XinliangN三、三、 MPI的通信域和组的通信域和组1.预定义通讯域预定义通讯域 MPI_Comm_World : 包含所有进程的组包含所有进程的组2.通讯域的分割通讯域的分割 MPI_Comm_Split(comm,color, key,New_Comm ) 02143576891011Color 相同的进程在同一组相同的进程在同一组根据根据key的大小排序的大小排序 （key相同时按原相同时按原ID排序）排序）例如：例如： 12个进程，个进程， 分成分成 3行行4列列Integer m

23、yid, 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_column)Call MPI_Comm_rank(Comm_Raw,myid_raw,ierr)Call MPI_Comm_rank(Comm_line, myid_line,ierr)MPI_Comm_WorldRAWCol

24、umnColor, 分组标准Key, 排序依据如相同，按原ID排 提交新定义的组提交新定义的组（否则新组无效，不要忘记）（否则新组无效，不要忘记）计算行号、列号19Copyright 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 得到组内编号得到组内编号 3)建立三个方向通讯的数据结构建立三个方向通讯的数据结构4

25、） 通信通信 ， 计算内点差分计算内点差分5） 计算边界差分计算边界差分02143576891011MPI_Comm_World20Copyright 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_Z, request(12),.Call MPI_Comm_Rank(MPI_Comm_World,myid,ierr)Call MPI_Comm_Split(MPI_Comm_World, mod(m

26、yid,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_Rank(Comm_Y,id_y,ierr)Call MPI_Comm_Rank(Comm_Z,id_z,ierr)定义三个方向的通信域定义三个方向的通信域21Copyright by Li Xi

27、nliangCall 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). 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) i

28、d_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_X,request(3),ierr)Call MPI_Irecv(A(M1+1,0,0),1,TYPE_X,id_X_Pre,99,Comm_X,request(4),

29、ierr) 定义新的数据结构定义新的数据结构22Copyright 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)EnddoEnddoEnddo call MPI_Wait_All(12,request,status,ierr) do k=1,P1 do j=1,N1 Ax(1,j,k)=(A(2,

30、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 .内点内点边界点边界点23Copyright by Li Xinliang四、分布数组的文件存储四、分布数组的文件存储 分布数组分布数组 real A(M/m1,N/n1)real A(M/m1,N/n1) 存储方式存储方式1. 每个进程存储到独立的文件每个进程存储到独立的文件 real A(M/m1,N/n1) character(len=50) filename write(filename,”(file-I4.4.dat)”) my

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

32、：缺点： 改变处理器规模时，需要处理改变处理器规模时，需要处理存储方式存储方式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(i,j ) 把把A1 装配到装配到A0 enddo Write(33) A0 else call

33、MPI_Send(A,) endif 01230123025Copyright 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列列优点：优点： 存储的数据形式与内存中存储的数据形式与内存中A0的存放格式一致。的存放格式一致。 存储的文件串行程序可直接读

34、取存储的文件串行程序可直接读取 real A(M,N) open(55,file=“A.dat”,form=“binary”) read(55) A close(55)26Copyright 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_fil

35、e_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 数目数目 27Copyright by Li XinliangPart 3 实例教学实例教学 CFD程序的程序的MPI实现实现实例实例 （1） 用拟谱方法求解不可压用拟谱方法求解不可压N-S方程方程

36、实例（实例（2） 用流水线方法计算紧致差分用流水线方法计算紧致差分 常用的优化方法常用的优化方法28Copyright 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(buf,count,data

37、type,dest,tag,comm, ierr)消息接收消息接收MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr) 29Copyright 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,

38、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 (logical型型)30Co

39、pyright 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,100,MP

40、I_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).31Copyright by Li Xinliang通讯域的分割通讯域的分割 MPI_Comm_Split(comm,color, key,New_Comm ) 02143576891011Color 相同的进程在同一组相同的进程在同一组根据根据key的大小排序的大小排序例如：例如： 12个进程，个进程， 分成分成 3行行4列列Line=mod(myid,3); raw=myid/3MPI_Comm_Spl

41、it(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_World32Copyright by Li Xinliang实例实例 1. 用（拟）谱方法求解二维不可压用（拟）谱方法求解二维不可压N-S方程方程2p物理模型物理模型周期性边界条件周期性边界条件按照给定能谱布置初始流动按照给定能谱布置初始流动 研

42、究流动的演化规律研究流动的演化规律33Copyright by Li XinliangFourier 变换变换 (1D)Fourier 变换变换 的特点：的特点： 求导数求导数 - 乘积乘积困难：困难： 非线性项非线性项卷积卷积计算量巨大计算量巨大在物理空间计算在物理空间计算Fourier 变换的快速算法变换的快速算法FFT34Copyright by Li Xinliang二维 Fourier 变换变换两次一维两次一维 Fourier 变换变换35Copyright by Li Xinliang求解步骤：求解步骤： 1） 读入初值读入初值 2） 调用调用FFT 得到得到 3） 计算计算 调用

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

44、y 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 end37Copyright by Li Xinliang数据重分布的实现数据重分布的实现A1(M/P,N)A2 (M,N/P)1234abcd对等式编程思想对等式编程思想 “我我”需要完成的工作需要完成的工作 1) 将数据将数据 A1(M/P,N) 切割成切割成

45、P块块 ，存入数组，存入数组B1(M/P, N/P,P) 2) 将数据将数据 B1(:,:,k) 发到发到 进程进程 k (k=0,1.P-1) 3) 从进程从进程k 接收接收 B2(:,:,k) 4) 组合组合 B2(:,:,k) 成成 A238Copyright 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=

46、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 问题：问题： 全部发送，全部发送， 发送成功发送成功后再启动接收。后再启动接收。 容易死锁容易死锁 按行分布按行分布 - 按列分布按列分布39Copyright by Li XinliangSubroutine R

47、edistibute_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 MPI_Recv(B2,M*N/(P*P),MPI_R

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

49、接收数目接收数目 sendtype 发送类型发送类型 recvtype 接收类型接收类型 Comm 通信域通信域 ierr 整数，返回错误值（整数，返回错误值（0为成功）为成功） To 0To 1To 2To 3Sendbuf 的数据格式的数据格式sendcountFrom 0From 1From 2From 3Recvbuf 的数据格式的数据格式recvcount41Copyright by Li Xinliang程序：程序：Subroutine Redistibute_ItoJ(A1,A2,M,N,P) Integer M,N,P,k,ierr,status(MPI_Status_Size

50、)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 问题：问题： 无法做到计算与通信重叠无法做到计算与通信重叠 42Copyright by Li

51、 Xinliang 二维二维 并行并行FFT 的实现的实现 （输入数据、输出数据均为按列分布）（输入数据、输出数据均为按列分布）1） 调用一维调用一维FFT实现实现 i- 方向的变换方向的变换 u - u12） 重新分布数据重新分布数据 （按列（按列- 按行）按行） u1 - u23) 调用一维调用一维FFT 实现实现j- 方向的变换方向的变换 u2- Fu24) 重新分布数据重新分布数据 （按行（按行 - 按列）按列） Fu2- Fu43Copyright by Li Xinliang实例实例 （2） 利用流水线利用流水线 实现紧致差分的并行化实现紧致差分的并行化紧致型差分格式：紧致型差分格

52、式： 相同网格点上引入更多信息。相同网格点上引入更多信息。 性能更优化。性能更优化。 是是 的差分逼近的差分逼近普通差分格式：普通差分格式： 显式给出显式给出 Fi 的表达式的表达式紧致型差分格式：紧致型差分格式： 隐式给出隐式给出 Fi 的表达式的表达式 6 阶中心6 阶对称紧致阶对称紧致 (Lele)5 阶迎风紧致阶迎风紧致 (Fu) j-2 j-1 j j+1 j+244Copyright by Li Xinliang 普通差分格式：普通差分格式： 直接计算导数，并行容易直接计算导数，并行容易紧致格式的计算：紧致格式的计算： 递推递推递推公式：递推公式：1)计算出计算出 （由边界条件或边

53、界格式给出）（由边界条件或边界格式给出）2） 由由 递推计算递推计算 出全部导数出全部导数 后面的数据必须等待前一步计算完成，无法并行后面的数据必须等待前一步计算完成，无法并行45Copyright 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(

54、k,N1) 缺点：缺点： 通信次数过多通信次数过多46Copyright 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 块块 47Copyright by Li Xinliang对称紧致

55、格式对称紧致格式追赶法追赶法令则代入（代入（1） 得得对比（对比（2） 得得边界处导数可由边界条件或边界格式给出边界处导数可由边界条件或边界格式给出： 则步骤：步骤： 1） 2） 由由 （3）式递推，得到）式递推，得到 3) 4) 由由 （2）式递推，得到）式递推，得到特点：特点： 两次递推。两次递推。 并行方法与前文类似并行方法与前文类似48Copyright by Li Xinliang3.常用的并行优化方法常用的并行优化方法 1） 通信与计算重叠通信与计算重叠 采用非阻塞通信采用非阻塞通信 Isend, Irecv 2） 用重复计算代替通信用重复计算代替通信 3） 拆分长消息、合并短消息

56、拆分长消息、合并短消息 4） 优化通信方式优化通信方式 49Copyright by Li Xinliang 用重复计算代替通信用重复计算代替通信 例如：例如： 计算差分计算差分 u 分布存储分布存储, f(u) 为为 u 的函数的函数 01方法方法 1） 计算出计算出 v=f(u) 通信得到通信得到 uN+1, vN+1 计算差分计算差分方法方法 2） 计算出计算出 v=f(u) 通信得到通信得到 uN+1 （边界外）（边界外） 计算出计算出 vN+1 =f(uN+1) 计算差分计算差分 方法方法 2） 计算量大，通信量小计算量大，通信量小 当函数当函数 f(u)不复杂时，可提高效率不复杂时

57、，可提高效率1, 2 N N+150Copyright 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_SendMPI_RecvMPI_Recv51Copyright by Li Xinliang合并短消息合并短消息 do

58、m=1,100 call MPI_Send(A(1,m),1,MPI_real,1 ) enddo 改为改为 do m=1,100 B(m)=A(1,m) enddo call MPI_Send(B(1),100, MPI_Real, 1, ) 52Copyright 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) 后再散发后再散发 53Copyright by Li Xinliang参考文献：参考文献： 郁志辉：郁志辉： 高性能计算之并行编程技术高性能计算之并行编程技术 MPI并行程序设计并行程序设计 迟学斌：迟学斌： 高性能并行计算高性能并行计算 （) 联系方式：联系方式： 李新亮李新亮 8254380154Copyright by Li Xinliang