HPC-基础并行数值算法.ppt

工程与科学中的高性能计算—并行数值算法基础 张健飞 河海大学力学与材料学院 Outline 概述 矩阵向量相乘 矩阵相乘 线性方程组求解 基本代数库 概述 算法是解题方法的精确描述，是一组有穷的规则，它们规定了解决某一特定类型问题的一系列运算。 第1节 矩阵向量相乘 预备——几个记号 第2节 矩阵相乘 串行矩阵乘法 A00B00 A00B01 A00B02 A11B10 A22B20 A11B11 A11B12 A22B21 A22B22 A01B10 A01B11 A01B12 A12B20 A20B00 A12B21 A12B22 A20B01 A20B02 A02B20 A02B21 A02B22 A10B00 A21B10 A10B01 A10B02 A21B11 A21B12 + + 第一步 第二步 第三步 ＝ A00B00＋ A01B10＋ A02B20＝C00 A00B01＋ A01B11＋ A02B21＝C01 A00B02＋ A01B12＋ A02B22＝C02 A11B10＋ A12B20＋ A10B00＝C10 A11B11＋ A12B21＋ A10B01＝C11 A11B12＋ A12B22＋ A10B02＝C12 A22B20＋ A20B00＋ A21B10＝C20 A22B21＋ A20B01＋ A21B11＝C21 A22B22＋ A20B02＋ A21B12＝C22 Cannon算法具有很好的负载平衡，其特点是在同一行中滚动，在同列中滚动，每一步计算出的部分值。从所给出的这五个并行计算矩阵乘积的算法可以看到，由于计算量对每个处理机来说是相同的，因此在选择算法时只需考虑通信量，对于方阵的乘积，当 时，Cannon 算法具有优越性。 第三节 线性代数方程组并行求解方法 串行算法 寻找置换矩阵Q，使得 线性方程组求解任务可以分为两方面：一方面是计算矩阵的LU分解，其中L和U分别是下三角矩阵和上三角矩阵；另一方面是求解三角形方程组。 Q是置换矩阵，L为单位下三角矩阵，U为上三角矩阵，采用in-place技术，将上、下分解因子分别存储在矩阵A得上、下三角部分中，替换掉A原来相应位置上的元素。 这里考虑的线性方程组是： 串行LU分解算法（列向选主元） 在上述算法中，主要计算工作量是修正矩阵A，即做 因此并行计算的主要任务就是在多处理机上同时对矩阵的不同部分做修正。对于每个j，对应的j越大，相应位置对应的计算量越小，所以对j连续地按块划分的方式进行并行计算是不合适的，将导致有的进程计算量大，有的进程计算量小，即出现负载不平衡。在多处理机上LU 分解的重要工作是使载荷尽可能的平衡，我们采用卷帘(wrap) 存储方式在各处理机上分配矩阵。把矩阵的第i列存放在Pi mod p中。 负载平衡 假设，在下面算法中n=p×m，A’的第i列为原来A的第(i-1)×p＋myid+1列，9×9矩阵在3个处理机中的存放方式为： 分布式系统的并行LU分解算法 以9列，3个处理器为例，即n＝9，p＝3，m＝3 以下为了便于演示，示意图中未考虑选主元; n=9,p=3 P0 j=0 P1 P2 P0 j=1 P1 P2 P0 j=2 P1 P2 P0 j=3 P1 P2 j=4 j=5 j=6 j=7 j=8 分布式系统上的三角方程组的并行解法 串行算法 在这个算法中每次对b进行修正时用到L的一列，如果按这种方式并行修正b，则称为列扫描方法。 讨论下三角方程组： 矩阵乘法中三个求和循环可任意排序，从而一共3！＝6种形式，如jki形式（标量描述） 这6种可能性（ijk,jik,ikj,jki,kij,kji）的任何一个都对应一内层运算，而且具有自己数据流动形式。例如ijk形式，内层运算是点积，数据用到的是A的行和B的列。 矩阵乘法：点积形式（ijk） 记 则上述算法可表示为 进一步可表示为 其中 或表示为 ijk形式的内部两层循环实质上是Gaxpy运算 矩阵乘法：Saxpy形式（jki） 假定矩阵A和C的列划分如下 比较两边第j列，可知 相当于一系列Saxpy运算 进一步，k循环实质上是一Gaxpy运算 得到如下算法 矩阵乘法：外积形式 考虑kji形式 记 则内部的两层循环是一个外积修正 可得到以下算法 可见其中的C=AB是p个外积之和。 并行矩阵乘法 假设m,k,n均能被p整除，即： 行列划分算法 将A和B矩阵分别划分为如下的行块子矩阵和列块子矩阵 算法： 算法中 , 在处理器中每次循环向前移动一个处理器。 初始状态，Ai、Bi和Cij(j=0,1,…,p-1)存放在Pi中。 为 矩阵 将Cij的计算放在Pi上进行，这