计算流体力学课件：MPI并行程序设计初步1

MPI并行程序设计（1）知识点：

MPI程序的运行机制——拷贝N份，给每个进程一份MPI的基本库函数——6个库函数

“对等式”编程方法——站在单个进程的角度思考12Part1:基本概念及MPI并行编程入门

并行计算基本概念MPI并行编程入门1）MPI简介2）MPI的基本函数及消息传递3）MPI的常用函数4）对等式编程思想的举例说明——全收集、矩阵相乘一、基本概念并行计算机简介大规模并行计算——超级计算（Supercomputing)/高性能计算(HPC)为何需要超级计算？

应用领域的巨大计算需求单CPU的计算能力有限

应用对计算资源的需求3

●CFD的计算资源依赖性

计算量大——流动的多尺度性（湍流）大飞机全部流动细节完全分辨：

最小尺度：mm~mm量级；

计算网格：1012-1016;需计算量：~1020-30

；工程需求：8个小时之内完成计算预计：LES：2045年；DNS:2080年最大尺度~mmm几种我们常见的并行计算机CPUCPUCPU…总线或交叉开关Memory(a)SMP,物理上单一地址空间CPUCPUCPU…定制网络LMLMLM虚拟分布共享存储(DSM)(b)DSM,逻辑上单一地址空间P/CP/CP/C…定制/标准网络LMLMLM(c)Cluster/COW,物理/逻辑上多地址空间SMPMPPMPP…WANLMDSMSM(d)Grid(ClusterofClusters)并行计算机体系结构5内存带宽瓶颈访存冲突机制控制复杂虚拟共享存储“NUMA”访存冲突机制控制复杂克服了访存冲突及内存瓶颈访存的局部性——对网络要求不严各系统的性能210021002100210021002100210021002100单处理器共享存储局域并行机群广域并行机群G

F

L

O

P

S6低价格可扩展自行搭建的简易机群7并行机群：搭建简单——简单的局域网并行机群=局域网早期作者搭建的简易机群机群软件：Linux/Windows;套件OSCAR;MPICH.NT，……我国最早搭建的机群：LSEC张林波搭建的32节点机8美洲虎/1700万亿次天河1号千万亿次CPU+GPU混合系统单精度千万亿次的GPU系统Mole-xx顶级的超级计算机天河2号：CPU+加速卡（XeonPhi)3.39亿亿次9Top5超级计算机（2013-6）排名计算机核心数Rmax(Tflops)Power(KW)1Tianhe-2(MilkyWay-2)-TH-IVB-FEPCluster,IntelXeonE5-269212C2.200GHz,THExpress-2,IntelXeonPhi31S1P

NUDT，China312000033862.7;178082Titan-CrayXK7,Opteron627416C2.200GHz,CrayGeminiinterconnect,NVIDIAK20x

CrayInc,USA56064017590.0;82093Sequoia-BlueGene/Q,PowerBQC16C1.60GHz,Custom

IBM,USA157286417173.2;78904Kcomputer,SPARC64VIIIfx2.0GHz,Tofuinterconnect

Fujitsu,Japan70502410510.0;126605Mira-BlueGene/Q,PowerBQC16C1.60GHz,Custom

IBM,USA7864328586.6;39452.并行程序设计工具1）共享存储式

自动并行（并行编译器）

IntelFortran/C编译器ifcaa.for-parallel

编译目标：多线程程序OpenMP内存CPU编译指示符：！ompparallel……10局域网络计算机Cluster系统2）分布存储式HPF(High-PerformanceFortran)

基于数据并行，程序改动较小

效率较低PVM(ParallelVirtualMachine)

MPI(MessagePassingInterface)

基于消息传递

效率较高11

MPI的编译、运行环境1）并行计算机（力学所机群、深腾7000，曙光5000A)

编译：mpif90/mpif77/mpiccf90/f77-I/usr/local/mpi/include-L/usr/local/mpi/lib-lm-lmpi

运行：mpirun/bsub2）MPIforNT(Windows2000,XP)

编译环境：VisualFortran/MSDevelopStudio

设置：头文件路径、连接路径

运行：mpirun二、MPI并行编程入门1.简介1213设置Windows下的MPI环境Step1:下载并安装mpich.nt.1.2.5安装包；Step2:更改VisualFortran的环境设置，添加MPICH的include及lib路径

1）Tools->options->Directories；

在“showdirectoriesfor:”

栏目选择“includefiles”;在“Directories:”下的对话框里面添加MPICHinclude的路径，例如“C:/Porgramfiles/mpich/SDK/include”(根据安装的具体位置而定）

在“showdirectoriesfor:”的栏目选择“Libraryfiles”,在“Directories:”下的对话框里面添加MPICHLib的路径，例如“C:/Porgramfiles/mpich/SDK/lib”

2)程序编译时，请把mpich.lib添加在链接库里。project->settings->link;在objcet/Librarymodules下的对话框里面添加mpich.lib（例如“kernel32.lib”变更为“kernel32.lib；mpich.lib”)Step3:编译程序，生成可执行文件14Step4:利用mpirun运行程序。（该工具在Windows桌面的“开始->程序->mpich->mpd->MPIRun”)

在“Application：”对话框里面选择刚编译生成的可执行文件。在“NumberofProcesses”对话框里面选择运行运行的进程数（即所谓“CPU个数”）。

在“Advancedoptions”对话框里面选择“Alwayspromptforpassword”。MPIRun运行时会要求用户输入计算机的用户名和密码。

点击“Run”即可运行(需要输入计算机的用户名和密码）。注意：如果程序为f90程序，请修改mpif.h。将行首的注释符“C”替换为“!”，否则编译会出错。（mpif.h在安装路径的include目录下，通常在C:/Porgramfiles/mpich/SDK/include里面）通常MPIRun需要以计算机管理员的身份运行，需要对计算机设置用户名和密码。如果计算机没有设置密码，则需要在控制面板中设置。

些防火墙及杀毒软件会阻止MPIRun的运行，如出现问题需要关闭这些防火墙及杀毒软件。服务器/前端机计算节点a.exea.exea.exeMPI程序的运行原理：

服务器（前端机）编译

可执行代码复制

N份，每个节点运行一份

调用MPI库函数

得到每个节点号my_id

根据my_id不同，程序执行情况不同

调用MPI库函数进行通讯MPI编程的基本思想：主从式，对等式重要！15重点：对等式程序设计16计算节点a.exea.exea.exea.exe对等式设计“对等式”程序设计思想如果我是其中一个进程；我应当做……完成我需要完成的任务站在其中一个进程的角度思考一个简单的MPI程序hello.f90include‘mpif.h’integermyid,numprocs,ierr

callMPI_Init(ierr)callMPI_Comm_rank(MPI_COMM_WORLD,myid,ierr)callMPI_Comm_size(MPI_COMM_WORLD,numprocs,ierr)Print*,"HelloWorld!myidis:",myid!添加自己的程序……!callMPI_Finalize(ierr)end17运行结果：18基本MPI函数（MPI子集）1)MPI初始化callMPI_Init(ierr)(out)Integer::ierr2)MPI结束callMPI_Finalize(ierr)(out)Integer::ierr19(in)：输入参数；(out)：输出参数；

整数，返回值非0表示出错3)得到当前进程标识callMPI_Comm_rank(MPI_COMM_WORLD,myid,ierr)(In)Integer::MPI_COMM_WORLD为进程所在的通信域(Out)Integer::myid,ierr

4)得到通信域包含的进程数CallMPI_Comm_size(MPI_COMM_WORLD,numprocs,ierr)(In)Integer::MPI_COMM_WORLD(Out)Integer::numprocs,ierr20进程的ID号（从0开始）最重要的参数！MPI_COMM_WORLD:MPI预定义的通信域；可换成创建的通信域Comm基本概念：通信域（通信子）01234567891011012345012345MPI_COMM_WORLDMPI_Comm_1MPI_Comm_221把全班分成几个组，执行任务更方便“班名”，包含全班同学MPI预定义my_id“学号”组的名字（编号）组内编号MPI消息传递函数

消息发送

MPI_Send(buf,count,datatype,dest,tag,comm,ierr)MPI消息：数据描述+信封数据描述：＜起始地址，数据个数，数据类型＞信封：＜源/目，标签，通信域＞22

buf:数据起始地址(Fortran:变量名，C：

变量地址/指针）count:数据数目（以datatype为单位，必须连续）

MPI_Send(buf,count,datatype,dest,tag,comm,ierr)23Datatype:数据类型

MPI_INTEGER,MPI_REAL,MPI_DOUBLE_PRECISION,MPI_COMPLEX,MPI_LOGICAL,MPI_CHARACTER,MPI_BYTE,MPI_PACKEDReal*8x(10)……（给x赋值）CallMPI_send(x(1),10,MPI_double_precision,…..)数据的首地址（不是变量的值）10个数据（不是10个字节）Fortran:按地址传送：x(1)或xC:按值传送：&(x[0])或xdest:发送目标的ID（integer）Tag:消息标签(integer)Comm:通信域(integer),例：MPI_COMM_WORLDierr:整数，如果成功返回0例：realA…….if(my_id.eq.0)CallMPI_Send(A,1,MPI_REAL,1,27,MPI_COMM_WORLD,ierr)

标签tag的作用：区别不同的消息MPI_Send(buf,count,datatype,dest,tag,comm,ierr)24标签。举例：0号同学向1号同学发出3封信，内容分别为3科考试成绩，为了防止弄混，必须用约定的数做标签。目的ID消息接收MPI_Recv(buf,count,datatype,source,tag,comm,status，ierr)

参数：数据起始地址，数据，类型，源地址，标签，

通信域，返回状态

integerstatus(MPI_STATUS_SIZE)MPI接收匹配：数据类型、源地址、标签要一致；否则不接收

数据发送缓冲区消息装配消息传递消息拆卸数据接收缓冲区25返回状态和调试信息MPI_Recv(buf,count,datatype,source,tag,comm,status，ierr)源地址，标签

26integerstatus(MPI_STATUS_SIZE)任意源地址:

MPI_ANY_SOURCE任意标签：MPI_ANY_TAG包含必要信息MPI预定义的常量可匹配任意源、任意标签“无论谁来的信都接收”status(MPI_SOURCE)：消息的源地址status(MPI_TAT)：消息的标签status(MPI_ERROR)：错误码MPI的消息发送机制——两步进行MPI_Send(A,…)发送MPI_Recv(B,…)接收

发送变量A接收到变量B配合使用发/收两步机制；避免直接读写对方内存；保证安全性!sumfrom1to100!Runonlyfornp=2!!!!include'mpif.h'integermyid,sum_local,sum_local1,sum_global,ierr,status(MPI_STATUS_SIZE)callMPI_Init(ierr)callMPI_Comm_Rank(MPI_COMM_WORLD,myid,ierr)sum_local=0

doi=myid+1,100,2sum_local=sum_local+ienddoprint*,'myid=',myid,'sum_local=',sum_local

if(myid.eq.1)callMPI_SEND(sum_local,1,MPI_INTEGER,0,99,MPI_COMM_WORLD,ierr)if(myid.eq.0)thencallMPI_RECV(sum_local1,1,MPI_INTEGER,1,99,MPI_COMM_WORLD,status,ierr)sum_global=sum_local+sum_local1print*,'sum_global=',sum_globalendifcallMPI_Finalize(ierr)end例：计算1+2+3…+100,采用2个CPU并行CPU0:1+3+5+7…+99CPU1:2+4+6+8…+100“对等式”编程思路：

站在每个进程的角度思考28结果显示29消息传递时要注意防止死锁缓冲区（“信箱”）MPI_SendMPI_RecvMPI_Send,MPI_Recv属于阻塞式发送/接收

发送成功

接收成功发送和接收是两个独立过程

子程序返回Send与Recv一定要配合好发送到“信箱”即为成功30重要！缓冲区发生死锁的情况只发送，不接收只接收，不发送缓冲区？31例，任务：进程0发送变量A给进程1

进程1发送变量B给进程0

if(myid.eq.0)thencallMPI_send(A,1,MPI_real,1,99,MPI_Comm_World,ierr)callMPI_recv(B,1,MPI_real,1,99,MPI_Comm_World,ierr)Elseif(myid.eq.1)thencallMPI_recv(A,1,MPI_real,0,99,MPI_Comm_World,ierr)callMPI_send(B,1,MPI_real,0,99,MPI_Comm_World,ierr)endif

Step1Step2Step1Step2不会死锁32死锁的例子

if(myid.eq.0)thencallMPI_recv(B,1,MPI_real,1,99,MPI_Comm_World,ierr)callMPI_send(A,1,MPI_real,1,99,MPI_Comm_World,ierr)Elseif(myid.eq.1)thencallMPI_recv(A,1,MPI_real,0,99,MPI_Comm_World,ierr)callMPI_send(B,1,MPI_real,0,99,MPI_Comm_World,ierr)endif

Step1Step1Step2不会死锁??Step233有可能死锁的例子

if(myid.eq.0)thencallMPI_send(A,1,MPI_real,1,99,MPI_Comm_World,ierr)callMPI_recv(B,1,MPI_real,1,99,MPI_Comm_World,ierr)Elseif(myid.eq.1)thencallMPI_send(B,1,MPI_real,0,99,MPI_Comm_World,ierr)callMPI_recv(A,1,MPI_real,0,99,MPI_Comm_World,ierr)endif

Step2Step2Step134使用MPI_Sendrecv()函数来避免死锁MPI_SendRecv(buf1,count1,datatype1,dest,tag1,

buf2,count2,datatype2,source,tag2,comm,status,ierr)=MPI_Send(buf1,count1,datatype1,dest,tag1,comm,ierr)+MPI_Recv(buf2,count2,datatype2,source,tag2,comm,status,ierr)次序由系统决定35基本的MPI函数（6个）——MPI的子集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)

MPI只能点到点通信，其他函数是用这6个函数开发的；使用这6个函数，可以实现MPI的全部功能。37系统时间（墙钟）函数:MPI_Wtime()

real*8::timetime=MPI_Wtime()返回当前墙钟时间（单位：秒）（Walltime与CPUtime不同，Walltime更可靠）CPU分时为所有进程服务Real*8::Time_begin,Time_end…(初始化)Time_begin=MPI_Wtime()….….（计算任务）….CallMPI_Barrier(comm,ierr)Time_end=MPI_Wtime（）Print*,“计算时间为：”,Time_end-Time_begin…可用来测量程序的执行速度测量加速比及并行效率加速比=N个进程的执行速度/单个进程的执行速度并行效率=加速比/N作者的Hoam-OpenCFD软件加速比测试CPUCorenumber1024-16384并行效率89.6%通常要进行同步，然后测量时间，否则各进程报出的时间不同。3.常用的MPI函数1）广播MPI_Bcast(buff,count,datatype,root,comm,ierr)

参数：数据缓冲区，数目，数据类型，根进程，通讯域

例：realAif(myid.eq.0)thenopen(55,file=‘data.dat’)read(55,*)Aclose(55)endifcallMPI_Bcast(A,1,MPI_REAL,0,MPI_COMM_WORLD,ierr)广播：树状传播，效率较高38……广播的逻辑图广播的实际实现方式——树状传播MPI的消息：只能点到点传递；不能真正“广播”2）规约（求和，求最大值，求最小值，……)

MPI_Reduce(sendbuf,recvbuf,count,datatype,op,root,comm,ierr)发送缓冲区，接收缓冲区，数目，数据类型，规约操作，通讯域!sumfrom1to100!Runonlyfornp=2!!!!include'mpif.h'integermyid,sum_local,sum_global,ierrcallMPI_Init(ierr)callMPI_Comm_Rank(MPI_COMM_WORLD,myid,ierr)sum_local=0doi=myid+1,100,2sum_local=sum_local+ienddo

callMPI_Reduce(sum_local,sum_global,1,MPI_INTEGER,&MPI_SUM,0,MPI_COMM_WORLD,ierr)print*,'sum_global=',sum_globalcallMPI_Finalize(ierr)end39预定义的规约操作：MPI_MAX最大值MPI_LXOR逻辑异或MPI_MIN最小值MPI_BXOR按位异或MPI_SUM求和MPI_MAXLOC最大值及位置MPI_PROD求积MPI_MINLOC最小值及位置MPI_LAND逻辑与MPI_BAKD按位与MPI_LOR逻辑或MPI_BOR安位或40同步MPI_Barrier(comm,ierr)Comm:通讯域，ierr:返回值等待所有进程都调用MPI_Barrier(),函数才能返回…………CallMPI_barrier(MPI_COMM_WORLD,ierr)……41424.“对等式”编程示例例1：全收集的实现MPI_Allgather()题目：N个进程，每个进程有一个数A；把所有进程的数收集起来，按进程号的次序形成数组A0(1:N),存放到所有进程中。

把每个同学的电话号码收集起来，形成通讯录，发给全班同学也可以是数组，各进程的数A可以不同

A0(1)

A0(2)

A0(3)

A0(4)方式1：

根进程收集所有数据;根进程发送到所有进程

if(myid.eq.0)thenA0(0)=Adoid=1,Nproc-1

callMPI_Recv(A0(id),1,MPI_Real,id,99,MPI_Comm_World,status,ierr)enddoelsecallMPI_Send(A,1,MPI_Real,0,…...)endifif(myid.eq.0)thendoid=1,Nproc-1callMPI_Send(A0,Nproc,MPI_Real,id,……)enddoelsecallMPI_Recv(A0,Nproc,MPI_Real,0,……)endif43“班长”依次与所有同学通信，收集信息；收集后依次通信，发放信息负载不均衡效率最低可能会死锁方式2：

根进程收集所有数据;根进程广播到所有进程

if(myid.eq.0)thenA0(0)=Adoid=1,Nproc-1

callMPI_Recv(A0(id),1,MPI_Real,id,99,MPI_Comm_World,status,ierr)enddoelsecallMPI_Send(A,1,MPI_Real,0,…...)endifcallMPI_Bcast(A0,Nproc,MPI_Real,0,MPI_Comm_world,ierr)

效率高于（1）是MPI_Allgather()的原有的标准方式44广播的实现方式“班长”依次收集信息后，“广播”给全班45……1AsendtoID00123N-1……1B0123N-1Step1:“我”（my_id进程）向my_id+1进程发数据;

“我”收my_id-1进程发来的数据该步完成后：

“我”（my_id进程）得到了my_id-1的数据全收集的实现图解方式3：循环通信46……2AsendtoID10123N-1……2B0123N-1Step2:

“我”向my_id+2进程发数据;

“我”收my_id-2进程发来的数据；该步完成后：“我”得到了my_id-1,my_id-2进程的数据Step3：我向my_id+3发数据，我收my_id-3发来的数据……StepN-1完成后，我得到了全部数据；

全体进程也得到了全部数据对等式编程思想：每个人做好自己的工作，全部工作就做好了不设班长，所有人工作量相同循环通讯:由张林波研究员首次提出

dostep=1,Nproc-1id_send=myid+step;if(id_send.ge.Nproc)id_send=id_send-Nprocid_recv=myid-step;if(id_recv.lt.0)id_recv=id_recv+NproccallMPI_Send(A,1,MPI_Real,id_send,99,&MPI_Comm_World,ierr)callMPI_Recv(A0(id_recv),1,MPI_Real,id_recv,99,&MPI_Comm_World,status,ierr)enddo

效率高于（1）(2)是MPI_Allgather()的现有的标准方式47

计算矩阵A*B=C，A,B,C：N*N矩阵

采用P个进程计算(N能被P整除）

存储方式：分布存储，A,C按行分割，B按列分割矩阵ABC48例2.计算矩阵乘积“对等式”程序设计思想：站在每个进程角度思考

“我”的数据：dimensionA1(N/P,N),B1(N,N/P),C1(N/P,N)

“我”的任务：计算C10,……K…P-1ACBA1B1C149需要得到整个矩阵B——自己只有Bk——向他人索取doid=0,P-1callMPI_Recv(B_tmp,N*N/P,MPI_REAL,id,……..)……计算出C1(id)……(C1(id)=A1*B_tmp)enddo只索取数据，何人提供？

“命令”他人提供数据？不符合“对等式”程序设计思想

0,……K…P-1ACB如何完成任务？B_tmp(N,N/P,P)C1(0)C1(1)C150B顺次发送数据dostep=0,P-1if(myid.eq.step)thendoid=0,P-1callMPI_Send(B1,N*N/P,MPI_REAL,id,……)enddoendifcallMPI_Recv(B_tmp,N*N/P,MPI_REAL,step,…)……计算出C1(id)……(C1(id)=A1*B_tmp)enddo

问题：负载不均衡实际上是串行

51一窝蜂地涌向同一个进程，负载不均衡，容易阻塞B按节拍循环发送数据（同“全收集”）

dostep