LecNote-13-MPI并行程序开发基础(示例程序)

第十三讲并行算法设计(一)并行算法设计的评判典型的并行应用问题矩阵乘FOX算法二维FFT问题Laplace问题、Jacobi迭代法多体问题N-Body的Barnes-Hut算法一维FFT问题Gauss-Seidel迭代法求解线性方程组常用的三种并行算法设计方法理想并行(embarrassinglyparallelcomputation)分治并行(partitioninganddivide-and-conquerstrategies)流水并行(pipelinedcomputation):本讲理想并行：矩阵乘FOX算法、二维FFT问题、Laplace问题、Jacobi迭代法下一讲分治并行：一维FFT问题、多体问题N-Body的Barnes-Hut算法流水并行：Gauss-Seidel迭代法求解线性方程组并行算法设计的评判:什么样的并行算法是好的什么样的并行程序是好的?效率高:使用N个处理器时,加速比能够接近N伸缩性好当N很大时,也能充分利用N个处理器的计算能力无论问题规模W多大,通过增加处理器的数量N,总是不需要修改程序就能让问题可解、计算需要的时间可以控制如何提高效率?按照BSP模型每个超级计算步上,各个处理器承担的负载要均衡要交换的数据量少超级计算步上不能有临界区如果要有临界区,就拆成多个超级计算步(例如SUM运算)如果不能用多个超级计算步来避免临界区,应该考虑修改计算算法了(N-body计算)临界区:概念上说,多个子任务要并发读写某个公共的数据结构,且每次只能有一个任务占有该数据结构的访问权;从实现上看,访问权的管理机制包括pthread模型:互斥锁操作,lock/unlockCELL/BE模型:MAILBOX操作,由PPE维护公共数据结构的值MPI模型:SEND-RECV操作,由master维护公共数据结构的值BSP模型:临界区与任务分配方法临界区:子任务间的相关读写相关:按照BSP模型,这两个子任务不在一个超级计算步上写写相关:这意味着子任务间写的次序不影响结果,也就是规约操作(reduction,如SUM/MAX/MIN).如果不是规约操作引起的写写相关,计算结果就不确定了.引入临时变量,保存每个子任务的局部结果多个超级计算步:第一个超级计算步计算各个子任务，后面N个超级计算步规约它们的结果,N是处理器数P的对数任务分配方法:并行程序实现阶段考虑的事情,如何优化超级计算步的实现效率peer-to-peercollaborating(静态任务分配):master-slave(动态任务分配):超级计算步上子任务的数量不能预先知道子任务的负载不一致使用的处理器性能不一致(以网络环境为计算平台时,会出现这种情况)求素数问题的任务分配实现:优化BSP的任务分配pthread使用临界区:MPI使用master-slaveBSP:并行程序的效率高超级计算步上没有临界区REDUCTION操作引起的临界区超级计算步上,各个处理器的负载要均衡:有足够的子任务子任务负载不均衡、或者各个处理器的性能不一致时,通过master-slave分配任务,可以达到基本负载的基本平衡超级计算步上子任务复杂性远大于它的数据通信开销需要动态任务分配时,子任务的复杂性远大于一次任务分配的开销

pthread模型:一次临界区操作CELL/BE模型:一次MAILBOX操作MPI模型:一次SEND-RECV操作BSP与计算机的体系结构基本没有关系分布存储计算机:消息传递模型实现数据交换、超级计算步的同步共享存储计算机基于cache的多核处理器、SMP：锁、信号量实现超级计算步的同步CELL处理器：DMA实现数据交换、mailbox实现同步BSP：并行程序的伸缩性好当N很大时,也能充分利用N个处理器的计算能力超级计算步上，要有足够多的子任务无论问题规模W多大,通过增加处理器的数量N,总是不需要修改程序就能让问题可解、计算需要的时间可以控制(Gustafson定律)超级计算步上，W增加了，子任务的数量也能增加，才能利用新增的处理器子任务的规模与W无关，否则处理器的寻址空间是个问题子任务的数据通信量与W无关如果有临界区操作,子任务的复杂性远超过临界区操作的开销规约操作引起的临界区操作动态任务分配引起的临界区操作并行程序要有好的伸缩性，需要面向分布存储的计算机共享存储计算机中CPU的数量很有限、总的存储空间有限BSP：好的并行算法超级计算步上，有足够多的子任务超级计算步上，或者没有临界区、或者只REDUCTION操作引起的临界区超级计算步上子任务的数量随着数据规模增加而增加子任务的数据规模、通信规模与总的数据规模没有关系子任务的计算开销远大于数据传输开销和临界区操作开销之和面向分布存储结构的计算机并不排斥使用共享存储结构的并行计算机：牺牲对问题规模的伸缩性好的并行算法是好的并行程序的基础，但不意味着并行程序一定好实现代码面向共享存储结构的并行计算机消息交换中会在问题规模大时出现死锁问题规模大时，消息交换出现死锁示例if(myRank<comm_size-1)MPI_Send(buf1,buf_size,MPI_INT,myRand+1,tag,comm);elseMPI_Send(buf1,buf_size,MPI_INT,0,tag,comm);if(myRank>0)MPI_Recv(buf2,buf_size,MPI_INT,myRand-1,tag,comm,&status);elseMPI_Recv(buf2,buf_size,MPI_INT,comm_size-1,tag,comm,&status);每个处理器都要将buf1中的数据交给MPIRuntime后，才执行MPI_Recv语句MPIRuntime需要缓存了comm_size-1个处理器中的buf1数据后，才能开始执行MPI_Recv语句处理器数量越大，意味着可并行执行的子任务数量越多，需要缓存的数据越多，容易出现缓存空间不够在SMP服务器上运行时，这种情况更容易出现：原来一台4路的SMP没有问题，新买一台8路的SMP后就运行不了了解决办法用非阻塞通信修改算法矩阵乘参考《AUser’sGuidetoMPI》第5.1节矩阵乘AB=CA是nm的矩阵B是mk的矩阵C是nk的矩阵并行性c(i,j)可以并行计算a(i,l)b(l,j)可以并行计算数据局部性：在计算c(i,j)的处理器上，需要A的第i行和B的第j列如何获得高性能？在数据划分与通信之间寻找trade-offFOX算法把处理器组织成一个二维的网格结构：

ppA、B、C都按照(block,block)划分FOX算法：循环从每一行处理器中，选择一个其中一个处理器，向所在行广播A的本地局部片段curA用广播得到的A片段curA与本地的B片段locB相乘，locC=locC+curAlocB把本地的B片段locB发送给“上邻居”，把从“下邻居”接收的数据作为本地的B片段locBA00A01A02A03A10A11A12A13A20A21A22A23A30A31A32A33B00B01B02B10B11B12B20B21B22B30B31B32C00C01C02C10C11C12C20C21C22C30C31C32=B03B13B23B33C03C13C23C33A00A01A02A03A10A11A12A13A20A21A22A23A30A31A32A33B00B01B02B10B11B12B20B21B32B30B31B32C00C01C02C10C11C12C20C21C22C30C31C32=B03B13B23B33C03C13C23C33A00A01A02A03A10A11A12A13A20A21A23A30A31A32A33B10B11B12B20B21B22B30B31B00B01B02C00C01C02C10C11C12C20C21C22C30C31C32B13B23B33B03C03C13C23C33A22B22=A00A01A02A03A10A11A12A13A20A21A22A23A30A31A32A33B20B21B22B30B31B32B00B01B12B10B11B12C00C01C02C10C11C12C20C21C22C30C31C32=B23B33B03B13C03C13C23C33A00A01A02A03A10A11A12A13A20A21A23A30A31A32A33B30B31B32B00B01B02B10B11B20B21B22C00C01C02C10C11C12C20C21C22C30C31C32B33B03B13B23C03C13C23C33A22B02=FOX算法分析FOX矩阵乘：代表了一种计算类型。有两个对象集合A和BA的每个成员与B的每个成员都有联系对A和B中的对象分组，每个作为一个子集：任取B中的一个子集，该子集分别与A的每个子集联系特点：可伸缩性(scalability)增加处理器的数量，可以解决更大规模的问题：数据划分单个处理器上，数组的大小有限制：内存空间、寻址空间增加处理器的数量，可以提高计算的效率计算开销：nm2kCmultiple+n(m-1)kCadd通信开销：p2CBStartup+nmCBmove+p2CPStartup+mkCPmovep2个处理器CBStartup：Broadcast通信的启动开销CBmove

：Broadcast通信中单个元素的开销CPStartup：“点到点”通信的启动开销CBmove

：“点到点”通信中单个元素的开销如何实现FOX算法？分布存储结构的并行计算机处理器组需要组织成一个二维的拓扑结构需要把处理器划分成子组，而且要有两种划分办法二维拓扑结构的每一行作为一个子组：广播A的数组片段二维拓扑结构的每一列作为一个子组：循环移动B的数组片段用“点到点”通信完全能够满足上述要求，只是编程比较麻烦MPI提供了一组通信子管理函数，这样编程起来比较方便进程组拓扑结构通信子的创建进一步理解MPI的进程组、通信域、通信子进程组：参与计算的进程，无拓扑结构。在代码设计阶段，度量伸缩性通信域：参与计算的处理器、何时参与、及其拓扑结构。代码实现阶段，表达对运行平台的要求通信子：进程组+通信子。代码运行阶段，满足程序的要求不同阶段，要不同的通信子多个并行程序的并发执行二维FFT什么是一维DiscreteFourierTransform：给定一个序列(a0,a1,…,an-1)，希望得到另一个序列(b0,b1,…,bn-1)，其中什么是二维DiscreteFourierTransform：给定一个nm的阵列(as,t)，希望得到另一个nm的阵列(cs,t)，其中

其中ei=cos+isine-i=cos-isin行变换列变换jkaj,kcj,kbj,k二维FFT分析计算特征：对nm的阵列(as,t)中每一行进行FourierTransform，得到一个新的nm阵列(bs,t)对nm的阵列(bs,t)中每一列进行FourierTransform，得到一个最终的nm阵列(cs,t)并行性从(as,t)nm变换到(bs,t)nm，每一行都是可以并行的数据局部性：计算bs,t时，涉及(as,t)nm的第s行全部元素从(bs,t)nm变换到(cs,t)nm，每一列都是可以并行的数据局部性：计算cs,t时，涉及(bs,t)nm的第t列全部元素同步：在开始计算(cs,t)nm前，必须完成(bs,t)nm的计算二维FFT的并行显然，(按照BSP)二维FFT问题的并行策略进程的拓扑结构是一个线性序列计算任务划分成两个步骤步骤1：(as,t)nm按照(block,*)的方式划分，每个进程负责局部(as,t)nm的行变换步骤2：(cs,t)nm按照(*,block)的方式划分，每个进程负责局部(cs,t)nm的列计算通信：步骤1结束后，(bs,t)nm按照(block,*)的方式划分步骤2要求(bs,t)nm按照(*,block)的方式划分需要对(bs,t)nm进行一次转置矩阵转置A0B0C0D0E0F0A1B1C1D1E1F1A2B2C2D2E2F2A3B3C3D3E3F3A4B4C4D4E4F4A5B5C5D5E5F5A0A1A2A3A4A5B0B1B2B3B4B5C0C1C2C3C4C5D0D1D2D3D4D5E0E1E2E3E4E5F0F1F2F3F4F5进程的拓扑结构是一个线性序列原分布：A[N,M]按照(block,*)的方式划分；各个进程中，局部的A片段按“列优先”方式存储转置在每个进程上，对自己的局部的A片段执行MPI_Scatter对每次MPI_Scatter所接收到的数组(Ai、Bi、Ci、Di、Ei、Fi)执行转置二维FFT的特征总结符合BSP模型整个计算问题从时间轴上可以划分成一组super-step，总的super-step数量是有限的在每个super-step上，有一组并行执行的子任务在相邻的两个super-step之间，数据访问的局部性不同，需要通信实际上，任何一个应用问题都具备下列特征从时间轴上，可以把整个计算问题的执行划分成有限个顺序执行的阶段(super-step)。例如在Laplace计算中，可把整个问题划分成三个步骤分发初始数据执行迭代计算收集计算结果每个阶段分别完成自己的任务，相对独立一个阶段的计算任务由一个进程完成：Laplace计算中的第一个阶段计算任务(读文件数据、组织成数组结构)和第二个阶段计算任务(写文件数据)都是由一个进程完成一个阶段的计算任务由多个进程并行完成：Laplace计算中的第二个阶段在相邻的两个super-step之间，数据访问的局部性不同，需要通信理想并行什么是理想并行？(embarrassinglyparallelcomputation)计算任务可以划分成一组子任务子任务之间是完全独立的——不需要交换数据每个子任务处理一份数据——每个子任务处理的数据可以相同，也可以不同子任务之间不需要交换计算的结果/中间结果对子任务的结果不需要额外的处理(如分发distribute、收集collect、合并combine等)，就可以得到整个任务的结果举例：二维FFT整个计算被划分成两个super-step，每个super-step上的计算任务都是理想并行的一个super-step可以划分成一组子任务，每个子任务执行一行(列)数据的一维FFT在单个super-step上，每个一维FFT的子任务处理完全不同的数据，没有任何的数据交换FOX算法也是理想并行FOX算法：循环从每一行处理器中，选择一个其中一个处理器，向所在行广播A的本地局部片段curA用广播得到的A片段curA与本地的B片段locB相乘，locC=locC+curAlocB把本地的B片段locB发送给“上邻居”，把从“下邻居”接收的数据作为本地的B片段locB在每个super-step上,要改变一次数据的划分输入数据计算结果进程Laplace问题典型的“场”问题：连续的数据区域（datadomain），每一点都有若干感兴趣的、随时间变化的物理量（速度，温度，压力等），且通常它们在时间和空间上都是连续的。“标准的”处理方法：在空间和时间上做离散化处理。

Laplace问题特征并行性：在时间轴上，每个时间步串行执行在空间轴上，数据区域上的每一点都是可以并行的数据局部性：在局域范围内，只与“邻居”有数据交换关系按照BSP，我们将每个时间步看作一个super-step每个super-step上的通信模式相同每个super-step上的计算划分相同每个super-step上的数据划分相同Jacobi迭代法一种线性方程组的数值解法什么是Jacobi迭代法：AX=B表示一个方程组A是mn的矩阵，且

A(i,i)0X是长度n的未知向量B是长度为m的已知向量令X(0)=(c0,c1,…,cn-1)当max(|X[i](t)-X[i](t-1)|)<时，X=X(t)Jacobi迭代法的理解在工程计算和科学分析领域，经常需要对系统的稳定性进行分析：在一定的约束条件下，我们关系的物理量在整个问题空间中的分布态势。例如：在考虑三峡大坝时，设计者必须关心的一件事情是三峡坝区的应力场分布，即这个区域中各个地质断层板块之间的稳定性大坝和蓄水都会增加有关地质断层板块上的压力在哪里建坝、建多大的坝、蓄水位最大可到多高时，这种新的压力不至于引起原有地质断层板块的稳定性我们不知道地质断层板块之间的作用力，但我们知道有哪些地质断层板块地质断层板块之间的拓扑结构山体的高度、各个地质断层板块的刚度、……人类已有的认识可以让我们确定一组约束条件：为了保持地质断层板块的稳定性，它们之间的各种作用力(压力、摩擦力等)与地质断层板块的特征(大小、刚度、位置、密度等)一定满足某个关系我们需要根据这些约束关系推算作用力在整个坝区的分布，从数学模型上考虑，可以把每个作用力看作是该系统中的一个自由度X的每个分量表示一个自由度A和B表示约束条件Jacobi迭代法的特征在每个时间步上，执行的计算任务完全相同并行性X[0](t)、X[1](t)、…、X[m-1](t)的计算可以并行进行数据局部性：计算X[i](t)

，涉及A的第i行和B[i]、X(t-1)数据相关性：每个X[i](t)都需要X(t-1)并行算法进程的拓扑结构是一个线性序列A按照(block,*)方式划分、B按照“block”方式划分对X(t-1)进行数据复制X(t)按照“block”方式划分通信：维护X(t-1)按照BSP模型分析把每个时间步上的计算看作一个super-step在单个super-step上时间的scalability：每个处理器有一次通信，执行n/p个X[i](t)的计算规模的scalability：A和B都被划分，只有X被复制运行的速度取决于问题的收敛速度——循环迭代的次数近似理想并行什么是近似理想并行？(nearlyembarrassinglyparallelcomputation)计算任务可以划分成一组子任务子任务之间是完全独立的——不需要交换数据每个子任务处理一份数据——每个子任务处理的数据可以相同，也可以不同子任务之间不需要交换计算的结果/中间结果经过对子任务结果额外的处理(如分发distribute、收集collect、合并combine等)，得到整个任务的结果举例：Jacobi迭代(求解方程组AX=B)，在每个时间步上的计算任务都是近似理想并行的每个X[i](t)的计算只涉及A的第