CFD2020第15讲MPI并行程序设计初步课件

计算流体力学讲义计算流体力学讲义2020 第第15讲讲 MPI并行程序设计并行程序设计（2）李新亮李新亮知识点：知识点：阻塞通信与非阻塞通信阻塞通信与非阻塞通信 非连续数据的发送与接收非连续数据的发送与接收 OpenMP OpenMP并行程序设计初步并行程序设计初步 1Copyright by Li Xinliang计算流体力学讲义2020知识点：1Copyright b服务器服务器/前端机前端机计算节点计算节点a.exea.exea.exeMPI 程序的运行原理：程序的运行原理：服务器（前端机）编译服务器（前端机）编译 可执行代码可执行代码复制复制 N 份，份，每个节点运行一份每个节点运行一份 调用调用MPI库函数库函数 得到每个节点号得到每个节点号 my_id 根据根据my_id 不同，程序执行情况不同不同，程序执行情况不同 调用调用MPI 库函数进行通讯库函数进行通讯MPI 编程的基本思想：编程的基本思想：主从式，主从式，对等式对等式2Copyright by Li Xinliang重点：对等式程序设计重点：对等式程序设计知识回顾知识回顾服务器/前端机计算节点a.exea.exea.exeMPI Copyright by Li Xinliang3推荐学习网站：推荐学习网站：“超算习堂超算习堂”（）easyhpc.Net 教材资料、在线课程；在线实训Copyright by Li Xinliang3推荐学习网Copyright by Li Xinliang4计算节点计算节点a.exea.exea.exea.exe对等式对等式设计设计“对等式对等式”程序设计思想程序设计思想如果我是其中一个进程；如果我是其中一个进程；我应当做我应当做完成我需要完成的任务完成我需要完成的任务站在其中一个进程的角度思考站在其中一个进程的角度思考Copyright by Li Xinliang4计算节点a基本的基本的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,datatype,dest,tag,comm,ierr)消息接收消息接收MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr)5Copyright by Li Xinliang5Copyright by Li XinliangMPI的消息发送机制的消息发送机制 两步进行两步进行MPI_Send(A,)发送发送MPI_Recv(B,)接收接收 发送发送 变量变量A接收接收 到变量到变量B配合使用配合使用6Copyright by Li XinliangMPI的消息发送机制 两步进行发送 变量A接收 到阻塞发送阻塞发送开始开始结束结束消息成功发出消息成功发出缓冲区可释放缓冲区可释放阻塞接收阻塞接收开始开始结束结束消息成功接收消息成功接收缓冲区数据可使用缓冲区数据可使用一、一、阻塞式通信与非阻塞式通信阻塞式通信与非阻塞式通信阻塞式发送与接收阻塞式发送与接收MPI_Send(A,)MPI_Recv(B,)7Copyright by Li Xinliang阻塞发送开始结束消息成功发出缓冲区可释放阻塞接收开始结束消息 MPI_Send()MPI_Send()返回后缓冲区可释放返回后缓冲区可释放 sum=sum=call MPI_Send(sum,)call MPI_Send(sum,)sum=sum=变量可重复利用变量可重复利用 MPI_Recv()MPI_Recv()返回后缓冲区数据可使用返回后缓冲区数据可使用Call MPI_Recv(sum1,)Call MPI_Recv(sum1,)Sum=sum0+sum1Sum=sum0+sum1 8Copyright by Li Xinliang MPI_Send()返回后缓冲区可释放8Copyri非阻塞发送非阻塞发送启动发送启动发送立即返回立即返回计计算算通信完成通信完成释放发送缓冲区释放发送缓冲区发发送送消消息息非阻塞接收非阻塞接收启动接收启动接收立即返回立即返回计计算算通信完成通信完成引用接收数据引用接收数据接接收收消消息息计算计算与与通信通信重叠重叠非阻塞消息发送与接收非阻塞消息发送与接收9Copyright 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,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型型)10Copyright by Li Xinliang非阻塞消息发送10Copyright by Li Xinli非阻塞通信调用后立即返回，缓冲区不能非阻塞通信调用后立即返回，缓冲区不能立即使用立即使用Sum=计算某变量计算某变量MPI_Isend(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 11Copyright by Li Xinliang非阻塞通信调用后立即返回，缓冲区不能立即使用MPI_Irec利用通信与计算重叠技术提高效率利用通信与计算重叠技术提高效率例：例：计算差分计算差分串行程序串行程序 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=N12Copyright by Li Xinliang利用通信与计算重叠技术提高效率例：计算差分0J=1,2,3并行程序并行程序 以两个进程并行为例以两个进程并行为例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_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)13Copyright by Li Xinliang并行程序 以两个进程并行为例01J=1,2 N/2AIf(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)=(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“内边界”01J=1,2 N/2特点：14Copyright by通信与计算重叠通信与计算重叠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,req2,ierr)call MPI_Isend(A(1,1),N,MPI_real,0,99,MPI_Comm_world,req1,ierr)endif01J=1,2 N/215Copyright by Li Xinliang通信与计算重叠01J=1,2 N/215CopyrighDo 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 Enddoendif01J=1,2 N/2特点：特点：传递边界信息传递边界信息 同时同时进行计算进行计算内点内点读取系统时间读取系统时间 doubleprecision time time=MPI_Wtime()16Copyright by Li XinliangDo j=2,N-101J=1,2 N/2特点：内点读二、二、如何收发非连续数据如何收发非连续数据例如：例如：发送数组的一行发送数组的一行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).17Copyright by Li Xinliang二、如何收发非连续数据A(1,1),A(1,2),A(方法方法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不足：不足：额外的内存占用额外的内存占用 额外的拷贝操作额外的拷贝操作通信不复杂的情况，内存拷贝工作量不大，该方法也可以采用。通信不复杂的情况，内存拷贝工作量不大，该方法也可以采用。效果还效果还可以可以18Copyright by Li Xinliang方法2.将发送的数据拷贝到连续的数组中不足：额外的内存占方法方法3：构建新的数据结构构建新的数据结构 Count:块的数量；块的数量；blocklength:每块的元素个数每块的元素个数Stride:跨度跨度（各块起始元素之间的距离）（各块起始元素之间的距离）Oldtype:旧数据类型，旧数据类型，Newtype:新数据类型新数据类型 （整数）（整数）例：例：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,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(3,4).A(1,:)在内存中的排列次序19Copyright by Li Xinliang方法3：构建新的数据结构 A(1,1),例：例：发送三维数组中的一个面发送三维数组中的一个面 (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,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 Xinliang例：发送三维数组中的一个面 (Fortran)20CMPI_TYPE_INDEXED(count,array_of_blocklengths,array_of_displacements,oldtype,newtype，ierr)构造数据类型更灵活的函数构造数据类型更灵活的函数 直接指定每块的元素个数及偏移量直接指定每块的元素个数及偏移量块的数量块的数量（整数）（整数）每块元素的个数每块元素的个数（整形数组）（整形数组）每块的偏移量每块的偏移量（整形数组）（整形数组）例：例：数组数组 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(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-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 XinliangNMPI_TYPE_INDEXED(count,array三、三、MPI的通信域和组的通信域和组1.预定义通讯域预定义通讯域 MPI_Comm_World:包含所有进程的组包含所有进程的组2.通讯域的分割通讯域的分割 MPI_Comm_Split(comm,color,key,New_Comm)1.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_column)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三、MPI的通信域和组02143576891011Col例：例：计算差分计算差分 三维分割三维分割 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）计算边界差分计算边界差分02143576891011MPI_Comm_World23Copyright by Li Xinliang例：计算差分 02143576891011MPParameter(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)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)定义三个方向的通信域定义三个方向的通信域24Copyright by Li XinliangParameter(M1=M/NM,N1=N/NN,P1=PCall 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)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_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定义新的数据结构25Copyright by Li Xinl 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,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内点边界点26Copyright by Li Xinlian四、分布数组的文件存储四、分布数组的文件存储 分布数组分布数组 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)”)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四、分布数组的文件存储012327Copyright by 分布数组分布数组 real A(M/m1,N/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(i,j)把把A1 装配到装配到A0 enddo Write(33)A0 else call MPI_Send(A,)endif 01230123028Copyright by Li Xinliang 分布数组 real A(M/m1,N/n1)012301存储方式存储方式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的存放格式一致。的存放格式一致。存储的文件串行程序可直接读取存储的文件串行程序可直接读取 real A(M,N)open(55,file=“A.dat”,form=“binary”)read(55)A close(55)29Copyright by Li Xinliang存储方式4.按列搜集后存储 Real Aj(M)第 1存储方式存储方式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,status,ierr)MPI_file_write_at(fileno,offset,buff,const,datatype,status,ierr)offset 偏移，偏移，buff 缓冲区，缓冲区，const 数目数目 30Copyright by Li Xinliang存储方式5 并行IO (MPI 2.0)30CopyriPart 3 实例教学实例教学 CFD程序的程序的MPI实现实现实例实例（1）用拟谱方法求解不可压用拟谱方法求解不可压N-S方程方程 实例（实例（2）用流水线方法计算紧致差分用流水线方法计算紧致差分 常用的优化方法常用的优化方法31Copyright by Li XinliangPart 3 实例教学 CFD程序的MPI实现31Cop回顾回顾 基本的基本的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,datatype,dest,tag,comm,ierr)消息接收消息接收MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr)32Copyright by Li Xinliang回顾 32Copyright 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,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型型)33Copyright by Li Xinliang非阻塞消息发送33Copyright by Li Xinli发送非连续数据发送非连续数据构建新的数据结构构建新的数据结构MPI_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)A(1,1),A(2,1),A(1,2),A(2,2).A(1,3).34Copyright by Li Xinliang发送非连续数据构建新的数据结构MPI_TYPE_VECT通讯域的分割通讯域的分割 MPI_Comm_Split(comm,color,key,New_Comm)02143576891011Color 相同的进程在同一组相同的进程在同一组根据根据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_World35Copyright by Li Xinliang02143576891011Color 相同的进程在同一组M实例实例 1.用（拟）谱方法求解二维不可压用（拟）谱方法求解二维不可压N-S方程方程2p物理模型物理模型周期性边界条件周期性边界条件按照给定能谱布置初始流动按照给定能谱布置初始流动 研究流动的演化规律研究流动的演化规律36Copyright by Li Xinliang实例 1.用（拟）谱方法求解二维不可压N-S方程2p物Fourier 变换变换 (1D)Fourier 变换变换 的特点：的特点：求导数求导数-乘积乘积困难：困难：非线性项非线性项卷积卷积计算量巨大计算量巨大在物理空间计算在物理空间计算Fourier 变换的快速算法变换的快速算法FFT37Copyright by Li XinliangFourier 变换 (1D)Fourier 变换 的特点二维 Fourier 变换变换两次一维两次一维 Fourier 变换变换38Copyright by Li Xinliang二维 Fourier 变换两次一维 Fourier 变换38求解步骤：求解步骤：1）读入初值读入初值 2）调用调用FFT 得到得到 3）计算计算 调用调用FFT 得到得到 4）计算计算 调用调用FFT 得到得到 5）计算计算 6）积分积分 求出下一时间步的值求出下一时间步的值 7）调用调用 FFT 得到得到 8）循环循环 3）-7）直到给定的时间直到给定的时间 39Copyright by Li Xinliang实际计算中，要采用抑制混淆误差的措施求解步骤：39Copyright by Li Xinlian程序的并行化：程序的并行化：二维二维 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(:)=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程序的并行化：二维 FFT二维FFT:调用两次一维数据重分布的实现数据重分布的实现A1(M/P,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数据重分布的实现A1(M/P,N)A2(M,N/P)123程序：程序：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_Real,k-1,.)A2(k-1)*M/P+1:k*M/P),:)=B2(:,:,P)Enddoend 问题：问题：全部发送，全部发送，发送成功发送成功后再启动接收。后再启动接收。容易死锁容易死锁 按行分布按行分布 -按列分布按列分布42Copyright by Li Xinliang程序：按行分布 -按列分布42Copyright bSubroutine 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 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 Xinliang43Copyright by Li Xinliang数据全交换：数据全交换：MPI_AlltoAll(sendbuf,sendcount,sendtype,recvbuf,recvcount,recvtype,comm,ierr)sendbuf 发送缓冲区（首地址）发送缓冲区（首地址）recvbuf 接收缓冲区（首地址）接收缓冲区（首地址）sendcount 发送数目发送数目 recvcount 接收数目接收数目 sendtype 发送类型发送类型 recvtype 接收类型接收类型 Comm 通信域通信域 ierr 整数，返回错误值（整数，返回错误值（0为成功）为成功）To 0To 1To 2To 3Sendbuf 的数据格式的数据格式sendcountFrom 0From 1From 2From 3Recvbuf 的数据格式的数据格式recvcount44Copyright by Li Xinliang数据全交换：To 0To 1To 2To 3Sendb程序：程序：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)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 问题：问题：无法做到计算与通信重叠无法做到计算与通信重叠 45Copyright by Li Xinliang程序：45Copyright by Li Xinliang 二维二维 并行并行FFT 的实现的实现（输入数据、输出数据均为按列分布）（输入数据、输出数据均为按列分布）1）调用一维调用一维FFT实现实现 i-方向的变换方向的变换 u-u12）重新分布数据重新分布数据（按列（按列-按行）按行）u1-u23)调用一维调用一维FFT 实现实现j-方向的变换方向的变换 u2-Fu2 重新分布数据重新分布数据（按行（按行-按列）按列）Fu2-Fu46Copyright by Li Xinliang 二维 并行FFT 的实现 1）调用一维FF实例实例（2）利用流水线利用流水线 实现紧致差分的并行化实现紧致差分的并行化紧致型差分格式：紧致型差分格式：相同网格点上引入更多信息。相同网格点上引入更多信息。性能更优化。性能更优化。是是 的差分逼近的差分逼近普通差分格式：普通差分格式：显式给出显式给出 Fi 的表达式的表达式紧致型差分格式：紧致型差分格式：隐式给出隐式给出 Fi 的表达式的表达式 6 阶中心6 阶对称紧致阶对称紧致 (Lele)5 阶迎风紧致阶迎风紧致 (Fu)j-2 j-1 j j+1 j+247Copyright by Li Xinliang实例（2）利用流水线 实现紧致差分的并行化紧致型差分格式 普通差分格式：普通差分格式：直接计算导数，并行容易直接计算导数，并行容易紧致格式的计算：紧致格式的计算：递推递推递推公式：递推公式：1)计算出计算出 （由边界条件或边界格式给出）（由边界条件或边界格式给出）2）由由 递推计算递推计算 出全部导数出全部导数 后面的数据必须等待前一步计算完成，无法并行后面的数据必须等待前一步计算完成，无法并行48Copyright 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)缺点：缺点：通信次数过多通信次数过多49Copyright 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 块块 50Copyright by Li Xinliang通信次数过于频繁解决方法：分块流水线步骤：F(kp,i对称紧致格式对称紧致格式追赶法追赶法令则代入（代入（1）得得对比（对比（2）得得边界处导数可由边界条件或边界格式给出边界处导数可由边界条件或边界格式给出：则步骤：步骤：1）2）由由（3）式递推，得到）式递推，得到 3)4)由由（2）式递推，得到）式递推，得到特点：特点：两次递推。两次递推。并行方法与前文类似并行方法与前文类似51Copyright by Li Xinliang对称紧致格式追赶法令则代入（1）得对比（2）得边界处导数3.常用的并行优化方法常用的并行优化方法 1）通信与计算重叠通信与计算重叠 采用非阻塞通信采用非阻塞通信 Isend,Irecv 2）用重复计算代替通信用重复计算代替通信 3）拆分长消息、合并短消息拆分长消息、合并短消息 4）优化通信方式优化通信方式 52Copyright by Li Xinliang常用的并行优化方法52Copyright by Li Xin 用重复计算代替通信用重复计算代替通信 例如：例如：计算差分计算差分 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)不复杂时，可提高效率不复杂时，可提高效率1,2 N N+153Copyright by Li Xinliang 用重复计算代替通信01方法 1）方法 2）方法 2 长消息切割成多个短消息发送、接收长消息切割成多个短消息发送、接收 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_Recv54Copyright by Li Xinliang 长消息切割成多个短消息发送、接收缓冲区MPI_Send合并短消息合并短消息 do 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,)55Copyright by Li Xinliang合并短消息 do m=1,10055Copy 优化通信方式优化通信方式 例：例：数据散发数据散发0号号 进程：进程：数据数据 A(100)，散发给散发给 0-99方式方式1）0 进程执行进程执行100次次 MPI_Send 其他进程执行其他进程执行 MPI_Recv MPI_Scatter()采用该算法采用该算法方式方式 2）0 进程进程 把把 A(100)切割成切割成10份份，发送给发送给10个进程个进程 10个进程接收个进程接收A1(10)后再散发后再散发 56Copyright by Li Xinliang 优化通信方式例：数据散发56Copyright by OpenMP并行编程入门并行编程入门一、一、特点特点 1.针对针对共享内存共享内存计算机结构计算机结构 全部全部CPU/线程均可访问内存线程均可访问内存 2.程序改动量小、实现方便程序改动量小、实现方便 （以编译指示符为主）（以编译指示符为主）3.适用于小规模并行或与适用于小规模并行或与MPI配合配合 进行大规模并行进行大规模并行 内存CPU（核心）CPU（核心）CPU（核心）1台PC机/1个计算节点（共享内存构架）CPU内存CPU内存CPU内存外部网络节点1节点2Cluster结构，分布内存构架OpenMP并行编程入门一、特点内存CPU（核心）CPU（print*,code 1!$OMP PARALLEL print*,code 2!$OMP END PARALLEL print*,code 3 end例1(test1.f90)：编译编译（在深腾（在深腾7000)运行结果（屏幕截图）ifort test1.f90 -openmp添加 openmp 选项运行：1.设置线程数（并行执行的数目）export OMP_NUM_THREADS=4 （例如，4个）2.执行：./a.out显示结果：code 1 code 2 code 2 code 2 code 2 code 3并行域中的代码执行了4次Test2.f90:print*,code 1!$OMP PARALLEL print*,code 2“!$OMP PARALLEL print*,“code 3”!$OMP END PARALLEL!$OMP END PARALLEL print*,code 4 end print*,code 1例1(test1.f93.DO 循环分解循环分解（openMP最常用的并行方法）最常用的并行方法）!$OMP PARALLEL!$OMP DO do k=1,12 print*,k enddo!$OMP END DO!$OMP END PARALLELend示例：线程线程0k=1,2,3线程线程1k=4,5,6线程线程2k=7,8,9线程线程2k=10,11,12!$OMP PARALLEL!$OMP DO!$OMP PARALLE