结构分析的并行预处理共轭梯度法

结构分析的并行预处理共轭梯度法 张朝阳 孙炳楠 唐锦春 ( )浙江大学土木系, 杭州, 310027 提 要 () () 共轭梯度法 法和预处理共轭梯度法 法由于其固有的并行性, 在有限元并行计算C G PC G 中得到很大的重视. 本文对 法和 法求解线性静力问题的并行算法给出了严密、清晰的数 C G PC G , 并提出了单元预处理矩阵的概念, 学推导过程, 给出了适合于工程应用的简单有效的预处理矩阵 有利于发展不同预处理矩阵的研究. 作者在局部内存式并行机 上实现了该算法.T ran sp u te r 关键词: 有限元; 并行算法; 共轭梯度法 中图法分类号: 242. 21O 0 引 言 随着科学研究与工程技术的发展, 大规模科学计算变得越来越普遍, 计算机硬件的发展也 是日新月异, 计算机体系结构打破了传统的冯r诺依曼体系结构, 走向了并行化的道路. 相应 地, 适合于各种结构的并行算法和并行化软件也不断出现. 在有限元分析中, 经常遇到有限元方程组 K x = f 的求解问题. 在使用迭代法求解时, 并 行化的共轭梯度法及预处理共轭梯度法将传统迭代步中每一步的计算分解到各单元上进行计 算, 避免了形成总刚阵, 使得小型机求解大规模问题成为可能, 所以这种方法得到了广泛的重 视. 文献1 仅给出了 C G 法的推导, 本文给出 PC G 法的推导, 并提出单元预处理阵的概念, 不 需形成总体预处理阵, 为并行预处理阵的研究提供了一个有效的方法, 并解决了文献4 中的 并行算法仍需形成总体预处理阵的问题. 文献2 、3 论述了 C G 法与 PC G 法在共享内存式 并行机上的并行算法, 本文给出在局部内存式并行机上的并行算法及简明的数学推导. 另外, () 很多文献 如文献1 、4 、5 中都未提到并行 C G 法和并行 PC G 法迭代初始值的选取, 本 文对这一问题作了论述. 1 有限元法 公式 小学单位换算公式大全免费下载公式下载行测公式大全下载excel公式下载逻辑回归公式下载 首先给出一些基本约定, 设: () e V — 相应于每个单元的单元向量, 即单元结点向量. () e e V — 各单元的V 按单元编号顺序直接组合而未经装配的向量. () e V — 表示各单元的V 经装配后的整体向量. E (e) V 表示V 对各单元分配而形成的向量, 即V 装配后再向各单元分配的量.(E ) E V — 表示V 相应于每个单元的单元向量, 即单元结点向量. 本文所说的“单元”的概念也包括“子结构”的含义. 注意本文约定的这些向量只是一种向 量定义, 不一定有物理意义, 但有限元方法中的结点力向量和结点位移向量也适用于这种定 义. 所以有限元原理的一些基本公式可表示如下: () () ()eeE () f = K x 1 E T e() 2x = A x f = A f e eE() 3f = K X e eE()4 f = K x T e ()5 f = A KA x = K x 其中A 在有限元方法中称为局部化矩阵, 当取单元坐标系与总体坐标系一致时, A 是布尔阵, T 其元素由 0 和 1 组成, 表达了单元间的连接信息; A 称为整体化矩阵. 故有以下关系式成立: E () V = AV 6 Te() 7V = A V E Te()8 V = AA V T Te 根据局部化矩阵A 和整体化矩阵A 的功能, AA V 显然等价于在整体上有关联的局部量求和 之后再将总体量向各单元分配, 即: , (e)(e)E ()) + V () 9 ( V V = ie6i?ad j. , (e) (e) 其中V ( )表示单元 e 的相邻单元的结点向量, V ()表示单元 e 的结点向量, 也不是简单的向i e 6 i?ad j. 量相加, 而是按“对号入座”的原则相加, 即整体上相关联的局部量相加, 以装配成整体向量, i? a d j. 表示单元 e 的所有相邻单元. 有了以上的约定, 我们就可推导C G 法与 PC G 法的并行算 法. 2 C G 法与并行 C G 法公式 ( 众 所周知, 解线性方程组 Kx = f 的共轭梯度法迭代过程如下 其中要求矩阵 K 是正定 ) 的: u = Kp () 10T T Α= rr p u ƒ() 11 x = x + Αp new () 12 r= r - Αu () new 13 TT () 14Κ= rrrrƒew new n ()15 p = r+ Κp new new () 式中 r 表示势能梯度或残余力向量 f - K x ; p 表示位移向量的共轭梯度方向. 开始迭代时, 取 初始值: () r = f - K x 160 () 17 = r p 一般取 x = 0.0 ()()() () 应用 6、7、8式对 C G 法进行并行化处理, 则10式成为: T eT eEu = A K A p = A Kp e e E T e () 若记 u= Kp , 则有 u = A u. 11式成为:P P ()e () () T T e T e Te T E eT E (( ( ( ) ) ) ) Θr r = A r r = rA r = rr = r r = 66e = 1e = 1 , ()()()() ) (() eeeeTE TE ) ( 式中 r = + r () , Θ = r r , P 为单元数, 所以内积 r r 的计算即可通过并行求解 ( ) r ei6 i?ad j.(e) 各单元的 Θ, 然后将其相加求得内积. 另外: T TT e T e E T e () () p u = p A u = A p u= p u ()() 将12、15等式的两边同乘局部化矩阵A , 得: E E E x = x + Αp w A x = A x + ΑA pnew ne 即: E E E A p = A r + ΚA pnew new p w = rw + Κp ne ne T e T e T e T e e e e e() ( ) ()() 而13式成为: A r= A r- ΑA u = A r- Αu, 可取 r = r- Αu . 14式算法与 11ew ew n n () 式相同. 初始迭代值若取 x = 0, 则16式为 r = f , 即:0 T e T e A r= A f e 这里须说明, f 是各单元的结点力向量, 但对于各单元相互关联的结点, 由于单元间的相 , e e (e) (e) (E ) (E ) 互作用, f 是未知的, 故 r 的相应分量可任选, 但须保证r )+ r ()= f , 其中f 的值为实( i e 6 i?ad j. e 际外载荷. r的取法在很多文献上被忽略. 这样可记为: ,() ()ee r = f , , () () () ()eeeE () )()( 其中 f 向量须满足f + f = f . 若将 17式方程两边同乘局部化矩阵A , 即有:i e 6i?ad j. E E p = r E E 此时 r 即取实际外载荷的值 f . 此时建立起以下并行 C G 法: e = 0x 0 ()18 ,e e()初始值: 19 r = f E E() 20p = r e eEu = K p ()21 e T E E T e) () ( Α= rrp u ƒ()22 EEE = x + Αp x new () 23eee = r ()- Αu r new 24 e T E e T E () ) ( ) 25( Κ= rr ƒrr ew w n ne E E E()26 p = r + Κpn ew new E e 迭代过程中 r 的值须由相关单元间的通讯实现 r 的迭加而得. 3 PC G 法与并行 PC G 法公式 () CG 法理理论上经过 n 步n 为矩阵 K 的阶数迭代一定能收敛于真实解, 但实际上, 由于它的计算过程对舍入误差十分敏感, 导致某些情况下求解有限元方程组收敛很慢, 甚至迭代次 数会超过 n. 为加速收敛, 人们提出对方程组系数矩阵进行预处理, 引入预处理矩阵, 即有以下 () 预处理共轭梯度法 PC G 法: u = K p () 27 T T Α= r d p u ƒ() 28 x = x + Αp new () 29 () r= r - Αu30new () 31d = M rnew new TT () 32Κ= rd ƒrd ew new n () 33p = d + Κp new new () 34r = f - K x 0 初始迭代值取: ()35 d = M r ()36 p = d- 1 一般取 x = 0, 其中M 称为预处理阵, 一般M 也须是正定矩阵, 且最好接近于 K , 可以看到,0 () PC G 法相对于 C G 法实际上只是多了 31式的运算. ()()()() 下面对 PC G 法进行并行化. 其中 27、29、30、34式的并行化与 CG 法并行化相同. ()()() () 应用 6、7、8式, 28式成为: T T e T e Te T E() ) ) ( ( r d = A rd = rA d = rd () 31式成为: T e T eT e Er A d = A M A r = A M newnew new e e E e e Te可记 d ew = M rw , 这里提出了单元预处理矩阵M 的概念, 且令M 满足A M A = M . 关于计n ne e e 算过程中子结构预处理阵M 的取法, 由其定义可知, 须使M 能装配成整体阵M , 本文中M 取 e J aco b i 预处理阵, 即取整体刚度矩阵对角元的倒数, 对应于各子结构相关联结点的M 的元素 应满足: , () ()ee)()( m = i e M + M 6j ?ad j. , () () ee ee 其中m 表示M 的元素,M 表示本子结构的M 的元素, M 表示相邻子结构的M 的元素 (e) ( i) 6i?ad j. e () () 按“对号入座”的原则求和. 所以M 可在满足上述条件下任意取值. 32式的并行算法同 28 () 式. 对33式等式两边同乘矩阵A , 即有: E E Ep = d + Κp new new () () () 35式并行算法同31式. 对36式等式两边同乘矩阵A , 即有: E Ep = d 综合上述, 有以下并行 PC G 法: ,e e()r= f 37 e e E () 38初始值:d = M r E e() p = d 39e eEu = K p ()40 E E ee T ) () ( Α= rd ƒp Tu () 41EEE= x + Αp new x () 42eee = r- Αu r()new 43 e Ee()= M 44 new r d e T E e T E () 45) ( ) ( Κ= rd ƒr d ew w n ne E E E()46 p = d + Κpn ew new EE ee 其中 d 、r分别通过 d 、r的通信求和获得. 本文推导的并行 PC G 法解决了文献4 中求向量 d 不能并行化的问题, 将迭代过程各步 e 骤完全分解到各单元上进行, 并且提出了预处理矩阵M , 的单元化矩阵M 的概念, 使得各单元 的预处理阵的选取更直观化, 有利于发展不同的预处理阵. 4 C G 法和 PC G 法并行算法的实现 从 P CG 法公式可以看出, P CG 法的每一次迭代过程包括下列运算: 一次矩阵向量乘法 T T () () Kp , 二次内积运算p u 和 rd , 三次比例向量运算 x = x + Αp、r = r - Αu 及 d = dnew new new ) ()+ Κp , 二次标量除法Α和 Κ的计算. 对于式 d = M r , 由于预处理阵M 的取法不同, 其计算方 法也不同. C G 法由于其固有的并行性, 已引起人们的很大兴趣, 但由于目前并行机的结构差别很 大, 而并行算法总是与并行机的结构紧密相连, 所以对于不同类型的并行机, 其并行 C G 法也 不同. 基地共享内存式并行机, 文献2 给出了椭圆型偏微分方程边值问题的有限元 PC G 解 法, 文献3 给出了基于 PC G 方法的结构优化问题的解法. 根据共享内存式并行机的内存特 点, 其并行算法中各向量以整体方式存储于共享式内存中, 但不需形成整体刚阵 K , 而以其单 (e) 元形式 K 存储. PC G 迭代步骤中最耗机时的计算为矩阵向量积算法. 文献2 、3 根据共享内存并行机 () () ee的特点, 对各单元循环计算 K p , 然后求和计算 p , 但由于存在内存访问冲突现象, 单元向量 对总体向量的叠加以及内积的求法不易并行化. 由于这个问题, 文献3 所获得的加速比不够 理想. 文献2 将结构划分为结构化网格, 并采用着色的方法以避免以存取冲突, 但这种方法 显然有局限性并增加计算量和复杂性. 本文所使用的是 Inm o s T 800 T ran sp u te r 并行机, 属于M IM D 局部内存式并行机, 各处理 器间的通讯由信息传递来实现, 不存在存取冲突的问题. 基于局部内存式并行机的特点, 矩阵 () e向量乘积可按下列方法实现: 即发送向量 r 到相邻单元所在的处理器, 并从该处理器接收向 ,( ) ()ie( ) () iE r + r . 内积的计算也可按第二节所述并行化. 量 r , 然后计算 r = 6i?ad j. 另外一个重要问题是预处理矩阵的选取及其并行化. 构造预处理矩阵的方法很多, 文献 ( )3 比较了不完全 分解法 、法、多项式法, 并对多项式法实现了并行算 C ho le sk y ICC G EB E 法, 但其方法较适合于共享内存式并行机, 并且由于其计算量较大, 耗费机时反而多于 法. C G 由于各种预处理阵各有其优缺点, 所以关于预处理矩阵的构造问题仍然有等于进一步探讨. 2 的预处理矩阵由于其形式简单并且易于并行化, 所以仍是一种有吸引力的方法. 文献J aco b i 并行算法即是采用 预处理矩阵. 本文亦采用 预处理矩阵, 在局部内存式并行 J aco b i J aco b i 机上实现并行 法.PC G 5 数值算例 为检验并行 法与并行 法的有效性, 作者对江苏某电视塔桁架结构作了风载荷静 C G PC G 2 ( ) 力分析 图1, 塔高134. 5, 其中天线杆39, 塔架95. 5, 采用基本风压0. 4ƒ, 塔架为四 m m m kN m 边形钢管结构, 共分16层, 详细数据参见文献6 . 所使用的 并行机具有4个 . T ran sp u te r C PU 将塔架结构划分为4个子结构, 对应4个 每个子结构有4层塔架. 连接的 , T ran sp u te r C PU C PU 拓扑结构如图2, 数值计算结果如表1: 图2 连接拓扑图T ran sp te r C PU 表1 数值计算结果 串行 法 并行 法 串行 法并行 法C G C G PC G PC G ()运行时间 秒 42. 611 12. 286 30. 690 9. 314 98 97 69 69 迭代次数 加速比 3. 47 3. 47 86. 8% 82. 5% 效率 从表1可以看到, 并行 法的效率略低于并行 法, 其原因 PC G C G 应是并行 法增加了处理器间的通讯量及同步点. 另外可以看到, PC G 使用 预处理阵的 法比 法节约时间30% 左右, 而且使 J aco b i PC G C G 图1 电视塔桁架结构 用 预处理阵编程简单, 所以不失为一种简单实用的方法. J aco b i 参 考 文 献 () 1 . , . . , 1986, 23 6: 845, 858L aw K HA P a ra lle l F in ite E lem en t So lu t io n M e tho dCom p u tS t ruc t 2 , . - - , . . . . . , 1988, 26:B a r ragy E e t a lA P a ra lle l E lem en tby E lem en t So lu t io n Sch em eIn tJN umM echE ng 2367, 2382 3 , . , Schm it L A e t a lS t ruc tu ra l O p t im iza t io n B a sed o n P reco nd it io ned Co n juga te G rad ien t A na ly sisM e tho d s . . . . . , 1994, 37, 943, 964 In tJN umM echE ng 4 姚坚. 工程结构分析中的并行计算方法研究: 学位 论文 政研论文下载论文大学下载论文大学下载关于长拳的论文浙大论文封面下载 . 北京: 中国建筑科学研究总院, 1993, 5 () 5 邓绍忠, 周树荃. 一类有限元结构分析的 计算方法. 计算结构力学及其应用, 1994, 11 3: 8EB E 6 王肇民, . 塔桅结构. 上海: 同济大学出版社, 1989P e il U Co n ju ga te g rad ien t m e tho d fo r s t ru c tu re p a ra lle l a lgo r ithm Zh an g C h ao yan g Su n B in gn an T an g J in ch u n (). , , , 310027D ep to f C iv il E ng inee r ingZh e jiang U n ive r sityH angzho u A bstrac t () () Co n ju ga te G rad ien t m e tho d C G an d P reco n d it io n ed Co n ju ga te G rad ien t m e tho d PC G a re ve ry im po r tan t m e tho d s in f in ite e lem en t p a ra lle l p ro ce ssin g b ecau se o f th e ir in h e ren t . p a ra lle lismT h e p a ra lle l a lgo r ithm s o f C G m e tho d an d PC G m e tho d fo r ca lcu la t in g th e lin ea r . ,sta t ic p ro b lem a re p re sen ted in th is p ap e rB e side s p re sen t in g a sim p le an d eff ic ien t m a t r ix w e a lso p ropo sed th e co n cep t o f e lem en t p reco n d it io n ed m a t r ix w h ich is h e lp fu l fo r th e p re2 . co n d it io n ed m a t r ix re sea rch in gA N um e r ica l ex am p le is ca lcu a ted o n a T ran sp u te r p a ra lle l .sy stem w ith d ist r ib u ted m em o ry : ; ; Key word sF EM p a ra lle l a lgo r ithm co n ju ga te g rad ien t m e tho d