稀疏线性方程组的求解共轭梯度算法步骤.ppt

2020/6/5,ParallelAlgorithmsChapter13NumericalAlgorithms,Spring,2018,2020/6/5,Solvinglinearsystemsofequations:Example:solvingordinaryorpartialdifferentialequationsnumericallyLeastsquaresproblems:Example:fittingcurvesorsurfacestoexperimentaldataEigenvalueproblems:Example:analyzingvibrations,StandardLinearAlgebraProblems,2020/6/5,StructureofA:denseGaussianelimination,etcSparseiterativemethodtriangularsubstitution,odd-evenreductioncertainPDEsmultigridDesiredAccuracy:Example:Gaussianelimination:moreaccurate,moreexpensiveConjugategradients:lessaccurate,lessexpensiveComputerSystemarchitecture,availablelanguages,compilerqualitylibraries?,ChoiceofAlgorithm,2020/6/5,讲授内容,13.1稠密线性方程组求解13.2稀疏线性方程组的求解13.3非线性方程的求根,2020/6/5,13.1稠密线性方程组求解,13.1.1一般求解方法1.线性方程组的定义和符号a1,1x1+a1,2x2+a1,nxn=b1a2,1x1+a2,1x2+a2,nxn=b2an,1x1+an,1x2+an,nxn=bn记为AX=b,2020/6/5,13.1稠密线性方程组求解,2.可并行化的求解方法(1)直接解法的并行化-Gauss消去法(包括选主元的Gauss消去法)-Gauss-Jordan消去法-LU分解法(2)迭代法的并行化(也可用于稀疏线性方程组)-Jacobi-Gauss-Seidel(可异步并行化)-JacobiOverRelaxation(JOR)-Gauss-SeidelOverRelaxation(SOR)-ConjugateGradient,2020/6/5,3.示例:上三角方程组的回代解法并行化(1)SISD上的回代算法Begin(1)fori=ndownto1do(1.1)xi=bi/aii(1.2)forj=1toi-1dobj=bj-ajixiaji=0endforendforEnd,13.1稠密线性方程组求解,可并行化,2020/6/5,13.1稠密线性方程组求解,(2)SIMD-CREW上的并行回代算法-划分:p个处理器行循环带状划分-算法Beginfori=ndownto1doxi=bi/aiiforallPj,where1jpdofork=jtoi-1steppdobk=bk-akixiaki=0endforendforendforEnd/p(n)=n,t(n)=n,2020/6/5,13.1稠密线性方程组求解,13.1.2SIMD-CREW上的Gauss-Jordan算法1.串行算法原理(一种直接解法)消元:通过初等行变换,将(A,b)化为主对角线矩阵,(方便起见,记b为A的第n+1列)求解:xj=aj,n+1/ajj注:选主元的Gauss-Jordan消去法,2020/6/5,13.1稠密线性方程组求解,2.SIMD-CREW上的并行算法(1)处理器:n2+n个处理器,这些处理器排成n(n+1)的矩阵,处理器编号为Pik,i=1n,k=1n+1(2)并行化分析消元的并行化:/O(n)forj=1ton-1,eachPikPar-do/第j次消元Pij(ij):aijj,k=j+1n+1):aikaik-ajk(aij/ajj)endfor求解:foreachPjj(j=1n)Par-do:xjaj,n+1/ajj/O(1)(3)时间分析:t(n)=O(n),p(n)=O(n2),c(n)=O(n3)成本最优？,2020/6/5,13.1稠密线性方程组求解,3.串行算法的最优时间由于x=A-1bA-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)=2i(n/2)+6m(n/2)+2a(n/2)a(n/2)=n2/2,m(n/2)=O(n/2)x)2i(n)=O(nx)综上,串行算法的最优时间为O(nx)2x2.5,2020/6/5,讲授内容,13.1稠密线性方程组求解13.2稀疏线性方程组的求解13.3非线性方程的求根,2020/6/5,Newtonnonlinearsolverasymptoticallyquadratic,Krylovacceleratorspectrallyadaptive,Schwarzpreconditionerparallelizable,ThreeVIPsforIterativeMethods,2020/6/5,13.1稠密线性方程组求解,13.1.3MIMD-CREW上的Gauss-Seidel算法1.串行算法原理(一种迭代解法)如果对某个k,给定的误差允许值c有则认为迭代是收敛的。2.并行化分析由于每次迭代需要使用本次迭代的前面部分值，因而难以得到同步的并行算法，下面给出一个异步的并行算法。,2020/6/5,13.1稠密线性方程组求解,3.MIMD异步并行算法-N个处理器(Nn)生成n个进程,每个进程计算x的一个分量-算法Begin(1)foralli:oldixi0,newixi0(2)生成进程i(3)进程irepeat(i)oldinewi(ii)newi(bi-kiaikoldk)/aiiuntili=1n|oldi-newi|cxinewiEnd,2020/6/5,13.1稠密线性方程组求解,-示例:P454例13.2-算法分析(1)从算法运行过程知,oldk需要同时读,即CREW;(2)读取/修改参数oldk时,需要一个信号灯Sk进行控制;(3)异步算法难以进行时间分析.,2020/6/5,13.2稀疏线性方程组的求解,13.2.1三对角方程组的求解1.Gauss消去法(难以并行化)消元回代,2020/6/5,13.2稀疏线性方程组的求解,2.奇偶规约求解法(可并行化)三对角方程可以写成如下形式fixi-1+gixi+hixi+1=bii=1nf1=hn=0-串行算法描述利用上下相邻方程消去偶序号方程中的奇下标变量:f2i-1x2i-2+g2i-1x2i-1+h2i-1x2i=b2i-1f2ix2i-1+g2ix2i+h2ix2i+1=b2if2i+1x2i+g2i+1x2i+1+h2i+1x2i+2=b2i+12i-1方程乘上某个数消去2i方程中的f2ix2i-1项,2i+1方程乘上某个数消去2i方程中的h2ix2i+1项,使2i方程变为ix2i-2+ix2i+ix2i+2=ii=1,2,n/2,2020/6/5,13.2稀疏线性方程组的求解,重复最终可得:case1:case2:g1x1+h1x2=b1.f2x1+g2x2+h2x3=b2.f3x2+g3x3+h3x4=b3.f4x3+g4x4=b4.可以分别得到g1x1+h1x2=b1或g1x1+h1x2=b1f2x1+g2x2=b2f2x1+g2x2+h2x3=b2f3x2+g3x3=b3直接解得x1,x2或x1,x2,x3回代求解x-并行化分析:、消去奇下标可以并行化；回代求解可以并行化,2020/6/5,13.2稀疏线性方程组的求解,13.2.2稀疏线性方程组的迭代解法1.迭代解法,2020/6/5,13.2稀疏线性方程组的求解,2.由PDE离散产生的稀疏线性方程组(1)Laplace方程,2020/6/5,13.2稀疏线性方程组的求解,(2)由五点格式的离散化(假设g(x,y)=0),2020/6/5,13.2稀疏线性方程组的求解,A为稀疏的块三对角矩阵,2020/6/5,13.2稀疏线性方程组的求解,3.Gauss-Seidel迭代解法的并行化(1)两种串行算法的计算顺序及其并行化顺序1顺序2注:顺序1难以并行化；顺序2可以小规模并行化,2020/6/5,13.2稀疏线性方程组的求解,(2)红黑着色并行算法并行计算所有的红点并行计算所有的黑点重复、直至满足精度要求,2020/6/5,13.2稀疏线性方程组的求解,4.共轭梯度法的并行化SOR等迭代方法依赖于很难选取的参数，而Hestenes-Stiefel共轭梯度法可避免此困难。下面考虑A为对称正定矩阵(1)共轭梯度法的基本思想将线性方程组AX=b的求解问题化为一个二次齐次函数的极小优化问题。,2020/6/5,13.2稀疏线性方程组的求解,AX=b与q(X)的关系,2020/6/5,13.2稀疏线性方程组的求解,重要定理,2020/6/5,13.2稀疏线性方程组的求解,(2)最优斜量法(最速下降法)求解最小优化问题，通常使用最速下降法，一般描述如下：,2020/6/5,13.2稀疏线性方程组的求解,最速下降法的迭代步骤：,2020/6/5,13.2稀疏线性方程组的求解,最速下降法的迭代图示：优点：.程序简单，每步迭代计算量少，存储省。.对于不太好的初始点X(0)，往往也能收敛。缺点：最速下降法是名不符实的，一般来说，它只有线性的收敛速度。,椭球面簇图：q(x)=Ciq(x1,x2)=Ci,2020/6/5,13.2稀疏线性方程组的求解,(3)共轭梯度法1952年，Hestenes和Stiefel发现了共轭梯度法，对于n阶的线性方程组，至多只要n步就能得到解。以n2为例从椭圆q(X)=q(X(0)上的X(0)点出发，用最优斜量法得到X(1)，d(0)=-r(0)=b-AX(0)(负梯度方向)r(1)=AX(1)-b,可以验证(r(0),r(1)=0,说明-r(1)是q(X)=q(X(1)的内法向(负梯度).第二步不走-r(1)方向，而走共轭方向p(1),可以验证p(1),过椭圆中心X*.p(1)方向的极小点就是X*.,2020/6/5,13.2稀疏线性方程组的求解,共轭梯度算法步骤：,2020/6/5,13.2稀疏线性方程组的求解,共轭梯度法的并行化分析(思考题)：考虑PRAM模型-数据划分每个处理器存有A的若干行，以及b、d、r的部分分量-计算和通讯计算：矩阵和向量乘积、向量内积时间？通讯：矩阵和向量乘积、向量内积计算时的数据交换时间？-总时间？,2020/6/5,13.2稀疏线性方程组的求解,稀疏系数矩阵的存储方法,2020/6/5,13.2稀疏线性方程组的求解,稀疏系数矩阵的存储方法,2020/6/5,讲授内容,13.1稠密线性方程组求解13.2稀疏线性方程组的求解13.3非线性方程的求根,2020/6/5,13.3非线性方程的求根,13.3.1SIMD-CREW上的等分求根算法1.二分法求根原理若f(x)在a,b上是连续的,且有f(a)f(b)0,则存在z(a,b)使f(z)=02.SISD上的二分求根算法:P458算法13.63.SIMD-CREW上的等分求根算法repeat将a,b分为N等分,等分点a=x0x1x2xN=b,Pi计算f(xi)存在某个f(xi)=0,输出xi若不成立,找xi,xi+1,使f(xi)f(xi+1)0,axi,bxi+1until|a-b|c并行时间O(logN+1|a-b|),2020/6/5,13.3非线性方程的求根,13.3.2MIMD-CREW上的Newton求根算法1.Newton求根法设f(x)在a,b上是连续且f(x)0,f”(x)0,f(a)f(b)0,取x0=b第n+1次迭代:xn+1=xn-f(xn)/f