中国科技大学并行计算课件10线性方程组的求解

并 行 计 算 中国科学技术大学计算机科学与技术系国家高性能计算中心(合肥)2004年12月Date1现代密码学理论与实践之五第三篇 并行数值算法 第八章 基本通讯操作 第九章 稠密矩阵运算 第十章 线性方程组的求解 第十一章 快速傅里叶变换 Date2现代密码学理论与实践之五第十章 线性方程组的求解 10.1 三角形方程组的求解 10.2 三对角方程组的求解 10.3 稠密线性方程组的求解 10.4 稀疏线性方程组的求解 Date3现代密码学理论与实践之五10.1 三角形方程组的求解 10.1.1 基本术语 10.1.2 上三角方程组的求解 Date4现代密码学理论与实践之五 基本术语 线性方程组的定义和符号 a1,1x1 + a1,2x2 + + a1,nxn = b1 a2,1x1 + a2,1x2 + + a2,nxn = b2 an,1x1 + an,1x2 + + an,nxn = bn 记为 AX=b Date5现代密码学理论与实践之五10.1 三角形方程组的求解 10.1.1 基本术语 10.1.2 上三角方程组的求解 Date6现代密码学理论与实践之五 上三角方程组的求解 上三角方程组的回代解法并行化 (1)SISD上的回代算法 Begin (1)for i=n downto 1 do (1.1)xi=bi/aii (1.2)for j=1 to i-1 do bj=bj-ajixi aji=0 endfor endfor End可并行化Date7现代密码学理论与实践之五 上三角方程组的求解 上三角方程组的回代解法并行化 (2)SIMD-CREW上的并行回代算法 - 划分: p个处理器行循环带状划分 - 算法 Begin for i=n downto 1 do xi=bi/aii for all Pj, where 1jp do for k=j to i-1 step p do bk=bk-akixi aki=0 endfor endfor endfor End / p(n)=n, t(n)=nDate8现代密码学理论与实践之五第十章 线性方程组的求解 10.1 三角形方程组的求解 10.2 三对角方程组的求解 10.3 稠密线性方程组的求解 10.4 稀疏线性方程组的求解 Date9现代密码学理论与实践之五10.2 三对角方程组的求解 10.2.1 直接求解法 10.2.2 奇偶规约法 Date10现代密码学理论与实践之五 三对角方程组的直接求解法 Gauss消去法(难以并行化) 消元 回代 注：由于三对角的 方程组的特殊性， 一次消元或一次 回代，只涉及邻 近一个方程，故 难以并行化。Date11现代密码学理论与实践之五10.2 三对角方程组的求解 10.2.1 直接求解法 10.2.2 奇偶规约法 Date12现代密码学理论与实践之五 三对角方程组的直接求解法 奇偶规约求解法(可并行化) 三对角方程可以写成如下形式 fixi-1+gixi+hixi+1=bi i=1n f1=hn=0 串行算法描述 利用上下相邻方程消去偶序号方程中的奇下标变量: f2i-1x2i-2+g2i-1x2i-1+h2i-1x2i =b2i-1 f2ix2i-1 + g2ix2i + h2ix2i+1 =b2i f2i+1x2i +g2i+1x2i+1+h2i+1x2i+2 = b2i+1 2i-1方程乘上某个数消去2i方程中的f2ix2i-1项, 2i+1方程乘上某个数 消去2i方程中的h2ix2i+1项, 使2i方程变为 ix2i-2+ix2i+ix2i+2=i i=1,2,n/2Date13现代密码学理论与实践之五 三对角方程组的直接求解法重复最终可得: case 1: case 2: g1x1+h1x2 =b1 . f2x1+g2x2+h2x3 =b2 . f3x2 +g3x3+h3x4=b3 . f4x3 +g4x4 =b4 . 可以分别得到 g1x1+h1x2 =b1 或 g1x1+h1x2 =b1 f2x1+g2x2 =b2 f2x1+g2x2+h2x3 =b2 f3x2+g3x3 =b3 解得 x1,x2 或 x1, x2, x3 回代求解x 并行化分析: 、消去奇下标可以并行化； 回代求解可以并行化Date14现代密码学理论与实践之五第十章 线性方程组的求解 10.1 三角形方程组的求解 10.2 三对角方程组的求解 10.3 稠密线性方程组的求解 10.4 稀疏线性方程组的求解 Date15现代密码学理论与实践之五10.3 稠密线性方程组的求解 10.3.1 有回代的高斯消去法 10.3.2 无回代的高斯-约旦法 10.3.3 迭代求解的高斯-赛德尔法Date16现代密码学理论与实践之五 有回代的高斯消去法 算法基本原理 求解过程分为消元和回代两个阶段，消元是将系数矩阵A化为上三角阵T，然后对TX=c进行回代求解。 消元过程中可以应用选主元方法，增加算法的数值稳定性。 下面是消元过程图：Date17现代密码学理论与实践之五 有回代的高斯消去法 并行化分析 消元和回代均可以并行化； 选主元也可以并行化； 消元过程的并行化图示：处理器按行划分Date18现代密码学理论与实践之五10.3 稠密线性方程组的求解 10.3.1 有回代的高斯消去法 10.3.2 无回代的高斯-约旦法 10.3.3 迭代求解的高斯-赛德尔法Date19现代密码学理论与实践之五 无回代的高斯-约旦法 串行算法原理 消元: 通过初等行变换,将(A,b)化为主对角线矩阵, (方便起见, 记b为A的第n+1列) 求解: xj=aj,n+1/ajj Date20现代密码学理论与实践之五 无回代的高斯-约旦法 SIMD-CREW上的并行算法 (1)处理器: n2+n个处理器, 这些处理器排成n(n+1)的矩阵, 处理器编号为Pik, i=1n, k=1n+1 (2)并行化分析 消元的并行化: / O(n) for j=1 to n-1, each Pik Par-do /第j次消元 Pij(ij): aij 0 Pik(ij, k=j+1n+1): aik aik-ajk(aij/ajj) end for 求解: for each Pjj(j=1n) Par-do: xj aj,n+1/ajj /O(1) (3)时间分析: t(n)=O(n), p(n)=O(n2), c(n)=O(n3) 成本最优？Date21现代密码学理论与实践之五 无回代的高斯-约旦法 成本最优？ 串行算法的最优时间：由于 x=A-1b A-1b(假设已有A-1): O(n2) 求A-1: 求A-1需要: 2次n/2n/2矩阵的逆 i(n/2) 6次n/2n/2矩阵的乘 m(n/2) 2次n/2n/2矩阵的加 a(n/2) i(n)=i(n/2)+6m(n/2)+2a(n/2) a(n/2)=n2/2, m(n/2)=O(n/2)x) 2x i(n)=O(nx) 综上，串行算法的最优时间为O(nx) 2x2.5 Date22现代密码学理论与实践之五10.3 稠密线性方程组的求解 10.3.1 有回代的高斯消去法 10.3.2 无回代的高斯-约旦法 10.3.3 迭代求解的高斯-赛德尔法Date23现代密码学理论与实践之五迭代求解的高斯-赛德尔法 串行算法原理 如果对某个k, 给定的误差允许值c有 则认为迭代是收敛的。 并行化分析 由于每次迭代需要使用本次迭代的前面部分值，因而难以 得到同步的并行算法，下面给出一个异步的并行算法。Date24现代密码学理论与实践之五迭代求解的高斯-赛德尔法 MIMD异步并行算法 N个处理器(Nn)生成n个进程, 每个进程计算x的一个分量 算法 Begin (1)oldi xi0, newi xi0 (2)生成进程i (3)进程i repeat (i) oldi newi (ii)newi (bi-kiaikoldk)/aii until i=1n| oldi - newi |c xi newi End Date25现代密码学理论与实践之五第十章 线性方程组的求解 10.1 三角形方程组的求解 10.2 三对角方程组的求解 10.3 稠密线性方程组的求解 10.4 稀疏线性方程组的求解 Date26现代密码学理论与实践之五10.4 稀疏线性方程组的求解 10.4.1 线性方程组的并行化方法 10.4.2 稀疏线性方程组的迭代解法 10.4.3 高斯-赛德尔迭代法的并行化 Date27现代密码学理论与实践之五线性方程方程的并行化方法 线性方程组选择算法的考虑因素 系数矩阵A的结构 dense Gaussian elimination, etc Sparse iterative method triangular substitution, odd-even reduction certain PDEs multigrid, etc 计算精度要求 Gaussian elimination: more accurate, more expensive Conjugate gradients: less accurate, less expensive 计算环境要求 architecture, available languages, compiler quality libraries?Date28现代密码学理论与实践之五线性方程方程的并行化方法 求解方法的并行化 (1)直接解法的并行化(用于稠密线性方程组) - Gauss消去法(包括选主元的Gauss消去法) - Gauss-Jordan消去法 - LU分解法 (2)迭代法的并行化(用于稠密和稀疏线性方程组) - Jacobi - Gauss-Seidel(可异步并行化) - Jacobi OverRelaxation(JOR) - Gauss-Seidel OverRelaxation(SOR) - Conjugate GradientDate29现代密码学理论与实践之五10.4 稀疏线性方程组的求解 10.4.1 线性方程组的并行化方法 10.4.2 稀疏线性方程组的迭代解法 10.4.3 高斯-赛德尔迭代法的并行化 Date30现代密码学理论与实践之五稀疏线性方程方程的迭代解法 迭代解法 Date31现代密码学理论与实践之五10.4 稀疏线性方程组的求解 10.4.1 线性方程组的并行化方法 10.4.2 稀疏线性方程组的迭代解法 10.4.3 高斯-赛德尔迭代法的并行化 Date32现代密码学理论与实践之五高斯-赛德尔迭代法的并行化 由PDE离散产生的稀疏线性方程组 (1)Laplace方程Date33现代密码学理论与实践之五高斯-赛德尔迭代法的并行化(2)由五点格式的离散化(假设g(x,y)=0)Date34现代密码学理论与实践之五高斯-赛德尔迭代法的并行化A为稀疏的块三对角矩阵 Date35现代密码学理论与实践之五高斯-赛德尔迭代法的并行化 Gauss-Seidel迭代解法的并行化 (1)两种串行算法的计算顺序及其并行化 顺序1 顺序2 注: 顺序1难以并行化；顺序2可以小规模并行化Date36现代密码学理论与实践之五高斯-赛德尔迭代法的并行化(2) 红黑着色并行算法 并行计算所有的红点 并行计算所有的黑点 重复、直至满足精度要求 Date37现代密码学理论与实践之五