LecNote-14-并行算法设计(二)

第十四讲并行算法设计(二)分治并行(partitioninganddivide-and-conquerstrategies)一维FFT问题多体问题N-Body的Barnes-Hut算法流水并行(pipelinedcomputation)Gauss-Seidel迭代法求解线性方程组一维FFT在前面的二维FFT算法中，我们没有考虑矩阵的每一行(列)是如何实现，只是假定已经有了一个串行的一维FFT/DFT函数，并没有考虑每一行(列)内部是否有并行性但是，如果计算任务是对一个一维的向量进行DFT，且n非常大，比如n=1M/G并行性：每个bk都是可以并行计算的局部性：每个bk的计算都需要整个的向量(a0,a1,…,an-1)从上一讲中，我们知道，为了提高并行的性能，必须对(a0,a1,…,an-1)进行划分提高程序对问题规模的scalability：最大可解问题的规模可以随处理器数量增加提高程序中的数据访问效率：通常，涉及的数据规模越小，访问的命中率就越高通信问题？bk=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15)(0,2,4,6,8,10,12,14)(0,4,8,12)(0,8)wk(4,12)wk(2,6,10,14)(2,10)wk(6,14)wk(1,3,5,7,9,11,13,15)(1,5,9,13)(1,9)wk(5,13)wk(3,7,11,15)(3,11)wk(7,15)a0a1a2a3a4a5a6a7a8a9a10a11a12a13a14a15b0b1b2b3b4b5b6b7b8b9b10b11b12b13b14b15a0a1a2a3a4a5a6a7a8a9a10a11a12a13a14a15b0b1b2b3b4b5b6b7b8b9b10b11b12b13b14b15bk=(0,1,2,3,4,5,6,7,8,9,10,11,12,13,14)(0,2,4,6,8,10,12,14)(0,4,8,12)(0,8)wk(4,12)wk(2,6,10,14)(2,10)wk(6,14)wk(1,3,5,7,9,11,13)(1,5,9,13)(1,9)wk(5,13)wk(3,7,11)(3,11)wk(7)a0a1a2a3a4a5a6a7a8a9a10a11a12a13a14b0b1b2b3b4b5b6b7b8b9b10b11b12b13b14一维FFT并行特征如果在存在关系的两个数据之间连接一条边，所有数据之间形成一个全连同的图每个处理器要执行log(p)次通信，其中p是处理器的数量每次通信交换数据量为n/p每一次通信时，都要对处理器进行重新分组每个处理器上的循环空间nlog(n)p每个循环的计算量：一维并行FFT的scalability每个处理器上的负载n/p个bk计算每个处理器上的数据量于p成反比通信的次数log(p)消息大小n/p分治并行partitioninganddivide-and-conquerstrategies划分策略(partitioning)：把一个问题分解成若干个组成部分通常，需要对每个组成部分的结果进行合成(combine)后，才能够获得整个问题的结果例如：理想并行采用的是一种划分策略，这种策略下，在把各个部分的结果合成为整个问题的结果时，几乎不需要什么额外的运算。划分策略的分类数据划分(datapartitioning)/域分解(domaindecomposition)：把对计算的数据分解成一组数据子集，在不同的子集上并行的执行处理。要求每个数据子集上执行的运算没有依赖关系不同的数据子集可以采用数据复制(datareplication)的策略，使它们有数据“重叠”例如在二维的FFT计算中，每一个super-step上，都执行的是数据分解，把整个二维阵列按照“行”/“列”划分，在每一“行”/“列”上并发执行一维FFT功能分解(functionaldecomposition)：把计算任务划分成一组独立的功能模块，并发执行每一个功能模块。不同的功能模块所需要的初始数据可以有“重叠”例如Jacobi迭代(求解方程组AX=B)，每个X[i](t)的计算作为一个功能模块，每个功能模块出了涉及A的第i行、B[i]为。，都需要X(t-1)(“重叠”的数据部分

)分治策略(divide-and-conquer)：把一个复杂的问题分解成一组子问题，每个子问题与原问题的形式相同，但规模比原问题小，而且对子问题可以采用这一策略进一步细分这是另一种特殊情形的partitioning策略：子问题与原问题除了问题规模外，完全相同例如：一维FFT(A)，其中A是一个长度为n的向量。我们采用的并行策略是令FFT(A)=FFT([FFT(A0dd),FFT(Aeven)])把FFT(A)划分为三个一维FFT子问题：FFT(A0dd)、FFT(Aeven)和FFT([FFT(A0dd),FFT(Aeven)]FFT(A0dd)执行对A的奇数位置的元素组成的向量的一维FFT，问题规模是原问题规模的一半FFT(Aeven)执行对A的偶数位置的元素组成的向量的一维FFT，问题规模是原问题规模的一半FFT([FFT(A0dd),FFT(Aeven)])执行对一个由两个元素所组成向量的一维FFTM-ary的分治策略把一个复杂的问题分解成M个子问题，M2，每个子问题与原问题的形式相同当子问题的规模还是太大时，进一步把每个子问题分解成M个粒度更细的子问题如此递归，直到每个子问题的规模合适为止一个4-ary分治的例子用分治策略解决N_Body问题有N个粒子(body)

在天体物理学里代表星体，例如地球、太阳、火星等在分子动力学里代表构成一个分子的各个原子……每个粒子有一定的状态：位置、速度、加速度、能量、温度等，这些状态是随时间变化的影响一个粒子状态变化的因素是粒子的初始状态其它粒子对它的作用力任意两个粒子之间都有牛顿万有引力，与两个粒子的距离有关对不同物理问题，粒子之间还存在其它与双方状态有关的作用力。例如在两个原子间还存在势能力(由原子的自转速度、原子间的距离等有关)。问题：对一个N_Body系统，给定其中各个粒子的初始状态最终状态是什么样的，即其中每个粒子的状态不再变化时，各个粒子的状态是什么样的？从初态到终态的变化过程中，变化的轨迹是什么样的？即对于我们关心的某(几)个状态量(例如电磁场势能的分布)，随时间变化的规律是什么样的？参考《DesigningandBuildingParallelPrograms》第1.4.2节把N个粒子按照“block”方式，划分成M组，M是处理器的数量每组作为一个粒子的子集，分配给一个处理器在每个处理器上，除了存储本地粒子的子集local_body_set外，开辟一个buffer，其容量大小是能够存储一个粒子的子集。在每个处理器上，用它的local_body_set对buffer进行初始化在每个时间步上，执行下列循环从1到M，执行循环对local_body_set中的每个粒子，计算buffer中各粒子对它的作用力把buffer中的内容发送给“左”邻居(0#处理器的“左”邻居是M-1号处理器)从“右”邻居接收消息，把消息数据存储在buffer中(M-1#处理器的“右”邻居是0#处理器)根据得到的作用力，更新local_body_set中的各个粒子的状态以100各粒子、4个处理器为例每个时间步上粒子0~24粒子0~24粒子25~49粒子25~49粒子50~74粒子50~74粒子75~99粒子75~99粒子0~24粒子25~49粒子25~49粒子50~74粒子50~74粒子75~99粒子75~99粒子0~24粒子0~24粒子50~74粒子25~49粒子75~99粒子50~74粒子0~24粒子75~99粒子25~49粒子0~24粒子75~99粒子25~49粒子0~24粒子50~74粒子25~49粒子75~99粒子50~740#处理器0#处理器0#处理器0#处理器N_Body问题的一种并行算法：性能分析对问题规模的scalability数据的scalability：计算数据以“block”方式划分在各个处理器上，数据存储的能力与处理器数量成线性关系速度的scalability：每个处理器上的主要运算开销是计算local_body_set各个粒子所受作用力，每个时间步上的复杂度为O(N2M)=O(N2)，其中N是粒子数量、M是处理器数量通信复杂度：在每个时间步上通信启动O(M)每次通信的数据传递开销O(NM)总的数据传递开销O(N)N是粒子数量、M是处理器数量因此，采用这种算法解决N_Body问题：随着问题规模的上升，无论是否增加处理器的数量，计算的性能都按照O(N2)的比例下降在大的分子模拟或天体物理学计算中，N的数量会很大，104或者更多怎样才能够提高实际应用问题的计算性能？N_Body问题一种优化并行算法Barnes-Hut算法：采用分治策略、进行近似计算基本思想：对任何一个粒子A，当一群粒子与A的距离“足够远”，这群粒子对A的影响可被一个“个体”来(近似)代替。在牛顿力学中，任何粒子B对A的作用与r2成反比，r是A与B间的距离对于一组粒子B1、B2、…、Bk，如果它们在空间位置上位于一个边长为d的立方体范围内，该立方体中心与A的距离为r。只要r“足够”大、d“足够”小，那么在计算B1、B2、…、Bk对A的影响时，可以把它们看作单个“个体”B来，用B对A的影响近似代替B1、B2、…、Bk对A的影响关键问题是：什么叫“足够远”、如何确定“一群”、用哪“一个体”来取代？N_Body问题一种优化并行算法：实现思路建树：确定一个包含所涉及空间的立方体；用一棵8叉树来表示对该空间的如下划分，把空间中的粒子分布的区域性和层次性同时表现出来。该立方体为根，如果其中只有一个粒子，停止；否则按照自然的方式，将它划分为8个相等的子立方体，如果其中有粒子，则用一个子节点代表。分别以子节点为根，重复上述过程。结果：每个叶节点代表一个粒子，每个粒子由唯一一个叶节点代表。内节点代表一个空间单元。可以考虑为节点实现为一个数据结构节点所代表立方体的质心的属性，包括空间位移、质量、速度、加速度等节点所代表立方体的边长d如果粒子分布在三维空间中，得到一棵8叉树如果粒子分布在平面空间中，得到一棵4叉树N_Body问题一种优化并行算法：实现思路计算一个粒子(叶节点)所受的力：从根开始做树的遍历，如果一个空间单元的质心距离本粒子足够远，则用在该质心的一个等价体来计算，不再考察空间单元下面的子树。足够远的标准：设一个立方体空间的边长为d，本粒子到它的质心的距离为r若则可以用该立方体包含星体的总质量和质心来计算，不需要再考虑下面的个别星体。其中，选在0.5—1.2之间，它越小，意味着近似的精度越高。这个式子也表达了“距离越远，能被近似的空间越大”的含义。

注意不同粒子的计算量不同对一个粒子而言，在遍历树的时候，“邻近”的分枝会遍历较深。由于树的平均高度为logn，计算一个星体的受力平均也就是logn。于是整个一个时间步的计算就是大约nlogn。N_Body问题一种优化并行算法：并行性分析采用Barnes-Hut算法解决N_Body问题时，在每个时间步上，都要依次执行下列四个super-step建树计算数上各个内节点的参数(质心、质量等)遍历树，计算各个粒子所受到的力更新粒子的属性

实际上，每个super-step都是可以并行的建树：不同的子空间可由不同的处理器负责向下分解，独立的操作计算数上各个内节点的参数：每个内节点的参数计算都完全独立遍历树，计算各个粒子所受到的力：以粒子为单位，完全独立的更新粒子的属性：以粒子为单位，完全独立的从各super-step的时间开销比例来看，计算各个粒子所受到的力是主要的时间开销N_Body问题一种优化并行算法：实现的难点涉及两个数据结构粒子，包括每个粒子的空间位置、质量、速度、加速度等属性通过建树过程得到的8(或4)叉树：在每个时间步上，得到的树不同如何划分这两个数据结构、如何分配计算任务建树过程中，自然的想法是：相邻的粒子位于同一个处理器上。但粒子是运动的，每个时间步上，粒子间的邻接关系会改变计算任何一个粒子所受的力时，都需要从数的根开始遍历，如何划分树？每个处理器上都保存一棵完整的树吗？计算各个粒子所受到的力时，自然的想法是每个处理器负责一组粒子的作用力计算，但是每个粒子作用力的计算复杂性不同在不同的时间步上，同一个粒子作用力的计算复杂性不同更新粒子的属性时，最好的分配办法是：把粒子以“block”方式均匀划分给各个处理器请考虑考虑数据划分的最简单办法是：每个处理器都保留这两个完整的数据结构但是，粒子的数量越大，树中节点的数量也越多，这样损伤了对问题规模的scalability比较二维FFT、一维FFT和FOX算法三者都可以用BSP模型刻画，而且super-step数量与被计算的数据值无关二维FFT可以表示成两个super-step一维FFT可以划分成log(p)个super-stepFOX算法中super-step取决于并行程序的进程拓扑结构Gauss-Seidel迭代法一种线性方程组的数值解法，比Jacobi迭代法的收敛速度快什么是Gauss-Seidel迭代法：AX=B表示一个方程组A是mn的矩阵，且

A(i,i)0X是长度n的未知向量B是长度为m的已知向量令X(0)=(c0,c1,…,cm-1)当max(|X[i](t)-X[i](t-1)|)<时，X=X(t)Gauss-Seidel迭代法的特征Gauss-Seidel迭代法的并行性分析与Jacobi迭代法不同的是，在每个时间步上，X[0](t)、X[1](t)、…、X[m-1](t)的计算必须是顺序执行的下列计算可以并行执行X[i](t)=(b[i]-left(i)(t)-right(i)(t-1))a[i,i]对于j>i，left(j)(t)=left(j)(t)+a[j,i-1]X[i-1](t)对于j<i，right(j)(t)=right(j)(t)+a[j,i-1]X[i-1](t)并行算法进程的拓扑结构是一个线性序列A按照(block,*)方式划分，left、right、B按照“block”方式划分X(t)按照“block”方式划分通信：每个时间步上，每个处理器执行p次广播MPI_Bcast其中一次是把局部的X(t)片段广播给其他进程p次用于接收其它处理器广播的X(t)片段流水并行(pipelinedcomputation)什么是流水并行(pipelinedcomputation)把一个问题分解成一组功能相同的子任务task，每个子任务处理的数据不同，这些子任务需要依次执行把每个子任务的执行划分成多个阶段stage每个进程负责实现一个stageP0stage0P1stage1P2stage2P3stage3P4stage4data0data1data2data3data4data5data6data0data1data2data3data4data5data0data1data2data3data4data0data1data2data3data0data1data2时间Gauss-Seidel迭代的流水并行什么是Gauss-Seidel迭代法：AX=B表示一个方程组A是mn的矩阵，且

A(i,i)0X是长度n的未知向量B是长度为m的已知向量令X(0)=(c0,c1,…,cm-1)当max(|X[i](t)-X[i](t-1)|)<时，X=X(t)X[0](t)、X[1](t)、…、X[m-1](t)的计算顺序执行下列计算并行执行X[i](t)=(b[i]-left(i)(t)-right(i)(t-1))a[i,i]对于j>i，left(j)(t)=left(j)(t)+a[j,i-1]X[i-1](t)对于j<i，right(j)(t)=right(j)(t)+a[j,i-1]X[i-1](t)并行算法进程的拓扑结构是一个线性序列A按照(block,*)方式划分，left、right、B按照“block”方式划分X(t)按照“block”方式划分第一步第二步第三步第四步第五步第六步第七步第八步规则计算与非规则计算对并行计算而言，我们可以这样理解(注意：我没有试图给一个精确的定义)规则计算：访问异地数据时，不需要使用间接地址划分子任务的数据时，不需要使用数据划分的索引非规则计算：不是规则计算的其他计算非规则计算的例子：线性稀疏方程求解一个大的稀疏方程，为了提高性能，常使用压缩存储的办法避免“0”元素的乘法：减少Jacobi迭代法或Gauss-Seidel迭代法中乘法运算量避免“0”元素占用的存储空间：提高数据访问效率例如用一个一维的数组A存储系数矩阵中的非“0”值，用另一个与A相同长度的数组B作为地址索引，B[i]记录A[i]在系数矩阵中的行号和列号采用前面的Jacobi迭代法或Gauss-Seidel迭代法并行算法计算，在划分方程组的系数矩阵时，需要根据B的值划分A的数据，B的作用是访问A的索引为什么会出现非规则计算的问题原因一：问题空间本身具有不规则性，通常表现为：有一组数据并行的子任务它们分别对不同的数据集合执行相同的运算但每个数据集合的规模不同例如：非连续弹簧体问题(multi-elasticbodysystem)每个弹簧体被离散化成一组节点，由于弹簧体的大小不同，节点的数量不同根据各个弹簧体边界上的受力，计算弹簧体的各个节点的移动情况12345678原因二：为了性能优化例如：稀疏线性方程的求解性能优化时常需要用到领域知识例如：对于分块矩阵形式的稀疏方程，我们可以把它划分成若干个非稀疏矩阵的小方程对于一般的稀疏矩阵，则必须采用索引的方法非规则计算的表现用两个不同的数组(数据文件)描述原始数据索引存储原始数据的数组/文件：一组定长的记录稀疏矩阵的一个非零元素：同一行/列的非零元素连续存放多弹簧体中，弹簧体上一个“刚点”：同一弹簧体的“刚点”连续存放存储索引的数组/文件：一组定长的记录稀疏矩阵：每一行/列用一个记录，这一行在原始数据(的数组/文件)中起始记录号有多少个记录多弹簧体：每个弹簧体用一个记录，这个弹簧体原始数据(的数组/文件)中起始记录号有多少个记录