（应用数学专业论文）分布式并行环境下的krylov子空间方法.pdf

摘要 本文提出了两种适合于分布式并行环境的，求解大型稀疏非对称线性方程组的并行迭 代法即改进的双共轭残差( b i c o n j u g a t er 艄i d l l a l ) 方法和改进的双共轭残差稳定化( b i _ c o n j u g a t er e s i d u a ls t a b i l i z e d ) 方法，分别简记为i b i c r ( i m p r o v e db i c o n j l l g a t er e s i d u 以) 方法和i b i c r s t a b ( i m p r o v e db i - c o l l j u g a t er e s i d u a ls t a b i l i z e d ) 方法作为解决k r y l o v 子空间方法并行化瓶颈的方法之一，我们通过算法重构，有效地将b i c r 和b i c r s t a b 方法中所需要的两个和三个全局同步化点分别降低到了一个，从而使所有内积计算以及 矩阵向量乘积是独立的，没有数据相关性我们的方法设计组合了数值稳定性与并行算 法设计的要素，代价仅是增加了一些计算量，且相对于通讯性能的改进，这些计算量的 增加是微不足道的性能和等效率分析表明i b i c r 和i b i c r s t a b 方法分别比b i c r 和 b i c r s t a b 方法具有更好的并行性和可扩展性理论分析和数值试验表明i b i c r 方法相 比于b i c r 方法的并行性能改进趋向于5 0 ，i b i c r s t a b 方法相比于b i c r s t a b 方法的 并行性能改进趋向于6 6 7 我们比较了i b i c r 与i b i c g 方法，b i c r 与b i c g 方法， 结果显示b i c r 和i b i c r 方法具有相同的收敛速度，且分别比相应的b i c g 和i b i c g 收 敛更快，i b i c r 方法还克服了i b i c g 方法出现的残差范数的振荡现象我们还比较了 i b i c r s t a b 与b i c r s t a b ，i b i c g s t a b 方法并行计算的收敛性，得到了i b i c r s t a b 方 法比b i c r s t a b ，i b i c g s t a b 方法收敛更快的结果 关键词：稀疏非对称线性方程组；k r v l o v 子空间方法；i b i c r 方法； i b i c r s t a b 方法；全局通讯；分布式并行计算 a b s t r a c t i nt h i sp a p e r ，j m p r o v e db i - c o n j u g a t er p s i d u 扎1 ( i b i c r ，j b r i e f ) a n d1 1 【1 p r o v e db i c o i q u g a t ei 沁s i d u a ls t a h i l i z e d ( i b i c r s t a b ，i nb r i e f ) l e t h o d sf o rs o l v i n gl a r g es p a r s el i i 卜 e a rs y s t e m sw i c hu n s y m m e t r i c a lc o e m c i e i l tm a t n c e sa r ep r 叩o s e df o rd i s t r i b u t e dp a r a l k l e n v i r o n m e n t sa so n eo fw a y st os o l v et h eb o t t l e _ n e c kp r o b l e mo fk r y l o vs u b s p a c em e 恤o d 8 o fp a r a l k l ，w cr e d u c et w oa n d 吐玛eg l o b 以s y n c h r o n i z a t i o np o i n t so ft h et w om e t h o d st o o n eb yr e c o n s t r u c t i n gb i c ra n db i c r s t a bm e t h o d sr e s p e c t i v 堍s u c ht h a ta ui n n e rp r o d u c t sp e ra n dm a t r i x _ v e c t ( ) rm u l t i p l i t _ i o n si t e r a t i o na r ei n d e p e n d e n t o u rm e h t o d sc o m b i n e t h ee l e m e n t so fi l u i n e i i c a is l a b i l i t yw i t ht h ec h 盯a c t e r so fd e s i g no fp a r a u e la l g o r 扎h m s t h ec o s ti so n l yal i t t l e m c r e a s e dc o m p u t a t i o n ，w h i c hl sap i e ( _ ：eo fc a k er e l a “v i n gt oi m p r o v e m e n t0 fc o m m u n i c a 七i o np e r f o r m a n c e p e r f o r m a n c ea n di s o e m c i e n ( 1 va n a l y s i ss h o w t h a ti b i c ra n dl b i c r s r a bm e t h o d sh a 阳b e t t e rp a r a l l e l i s r na n ds c a l a b i l i t yt h 卸b i c r a n db i c r s7 r a bm e t h o d s ，r e s p e c t i 、p l ht h e o r e t i c a la n a l y s i s 蛐dn m e r i c a le x p e r i m e n t s s h o wt h a tt h cp a r a l l e lp e r f o r m a n c ee a nb ei m p r o v e db yaf a c t o ro fa b o u t2a n d3 jr e s p e c t i v l h f l l r t h e r m o r e ，w ec o m p a r e di b i c rw i t hi b i c ga n db i c rw i t hb i c gm e t h o d s t h e r e 5 u l t ss h o w 仙a tb i c ra n di b i c rm e t h o d sc o n v e r 能n ta tt h es a m en u m b e ro fi t e r a t i o 。 a n dt h e yc o n v e r j g e n tf a s t e rt h a nb j c ga n di b i c gn l e h o d s ，a n di b j c ri n e t h 。dc o n q u e r s t h ev i b r a t i o no fr e s i d u a ln o r mo fi b i c gm e t h o dw ec o m d a r e da l s oi b i c r s t a bw i 乞hb i c r s t a ba n di b l c g s t a bm e t h o d s w e8 e et h a ti b i c r s t a bm e t h o dc o n v e r g e r l t sf a s t e r t h a nb i c r s l 、a ba l 埘i b i c g s t a bm e t h o d s k e yw o r d s ：s p a r s eu n s y m m e t r i c a ll i n e a rs y s t e m s ；k r y l o vs u b s p a c em e t h o 出；i b i c r m e t h o d ；i b i c r s t a bm e t h o d ；g l o b a lc o m m u n i c a t i o n ；d i s t r i b u t e dp a r 8 】l e lc o m p u t i n g 1 1 第一章研究背景和预备知识 1 1 研究背景 求解线性方程组是数值计算的一个基本任务在科学计算中，线性方程组是经常遇到 的比如在用有限差分法或有限元法近似求解偏微分方程问题的过程中，还有在求解非线 性问题的过程中作为中间步而出现由于科学计算问题的求解要求具有更高的精度和更快 的计算速度，所以问题离散后往往得到大型稀疏线性方程组而大型稀疏线性方程组的求 解时间通常在整个计算时间中占有很大的比重，有的甚至达到8 0 所以它的高效求解成 为了科学计算的关键目前流行的求解方法是迭代法和并行迭代法| 1 1 在大型稀疏线性方程组的迭代法中，k r y l o v 子空间方法是最有效的，如对称正定问题 的共轭梯度( c g ) 方法 2 】、非对称问题的g m r e s 3 ，b i c g 4 ，q m r 5 ，b i c g s t a b 6 ， b i c r 7 和b i c r s t a b 方法等k r y l o v 子空间方法的每个迭代步仅以矩阵向量乘积的 形式涉及系数矩阵在许多情形中，特别是当系数矩阵具有好的结构时，这些运算适合于 在向量和共享存储并行机上实现但是对分布式存储计算机，矩阵和向量被分布在各处理 器上，因此，即使通过并行运算可有效地实现矩阵运算，我们仍不能避免全局通讯，即内 积计算所需要的所有处理器间的通讯 所有k r y l o v 子空间方法，包括b i c r 和b i c r s t a b 方法，在分布式存储并行机上 的并行化主要涉及三种运算4 1 = 向量校正、矩阵向量乘积和向量内积向量校正是自然 并行的对大规模稀疏矩阵，矩阵向量乘积往往只需要相邻处理器问的局部通讯，例如用 五点或九点差分法离散所得的稀疏矩阵它与向量的乘积仅需要与其相邻的四个或八个处 理机进行局部通讯由于向量被分布在各处理机上，所以向量内积的计算所需要的聚集和 广播必然需要全局通讯当处理器数目增加时，这些全局通讯的开销变得越来越大，并且 潜在地对算法的可扩展性具有严重的负面影响，因此向量内积计算引起的全局通讯成为了 k r y l o v 子空间方法并行化的主要瓶颈关于分布式存储系统上的通讯问题的详细讨论可 参见 8 解决这种引起性能退化的瓶颈问题具有三种方式，当然它们可以组合使用第一种是 消去数据相关性或减少全局同步化点个数，使得几个内积可以同时计算，同时传送；第二 种是进行算法重构，使得尽可能多的通讯可与有效的计算进行重叠；第三种是用其它不涉 1 第一章研究背景和预备知识 及全局通讯的计算来替换内积计算本文的方式属于第一种，即减少全局同步化点个数 最近，b 讧c k e r 等 9 和1 r a n g 等 1 0 基于耦合的两项递归l a n c z o s 过程提出了拟极 小残差( q m r ：q u a s i - m i n i m a lr 船i d u a l ) 方法的一种新的并行形式；s t u r l e r 等f 8 1 给出了 如何降低g m r e s ( m ) 和c g 方法中全局通讯的影响；l i u 等 1 l 】提出了并行c r 方法 1 n g 等 1 2 1 4 分别给出了改进的c g s 、b i c g s t a b 和b i c g 方法；l i u 和g u 等f 1 5 1 提出了并行q m r g g s t a b 方法所有这些都着重在前两种方式g u ，l i u 和m of 1 6 1 提 出了一种不需要整体内积的c g 类方法，即多搜索方向共轭梯度( m s d c g ) 方法，基于 区域分解，该方法将c g 方法中的内积计算用小线性方程组的求解来代替从而完全消除 了全局通讯，这属于第三种方式 本文通过对b i c r 和b i c r s t a b 方法进行算法重构，提出了改进的b i c r 和b i c r s r a b 算法对比b i c r 和b i c r s t a b 算法，i b i c r 和i b i c r s t a b 算法的数值稳定 性分别和b i c r 和b i c r s t a b 算法相同，同时考虑了在分布式存储并行机上实现时并行算 法的性能，其同步开销分别减少为原来的二分之一和三分之一对比b i c r 和b i c r s r a b 算法，从理论上证明了i b i c r 和i b i c r s 7 r a b 的可扩展性都得到了提高，性能提高比率 分别趋向于5 0 和6 6 7 通过数值实验也得出了与理论分析相吻合的结果 1 2 1 一般投影算法 1 2 预备知识 考虑求解大型稀疏线性代数方程组 a z = 6 ( 1 2 ，1 ) 其中a r 。是非奇异矩阵，。，6 r 。v 是维向量 为求解线性方程组( 1 2 1 ) ，大多数实用的迭代法都这样那样地利用了投影方法，即从 维向量空间中找出一个子空间k 。，d i m k 。= m ，在其中寻找近似解，子空间匠。称为 求解子空间或搜索子空间为在k 。中求出一个近似解，显然要有”z 个约束条件，通常 采用m 个正交性条件特别地，可采用残差向量r = 6 一a z 与m 个线性无关向量正交的 2 第二章k r y l o v 子空间方法 了g m r e s r 方法，它涉及内积和外积两类方法，内积法即g m r e s 方法，外积法即g c r 方法1 9 9 3 年s a a d 提出了灵活的g m r e s 方法，即f g m r e s 方法这种方法使得每 步所用的预条件子可以不同，更具灵活性1 9 9 6 年s t u r l e r 基于g c r 方法给出了一大类 方法，这类方法的特点是将g c r 方法的正交性与其它算法，如g m r e s 、b i c g s t a b 等方法结合起来产生新的收敛性好的方法 属于正交化方法的主要k r y l o v 子空间方法的相互关系可如下表示： l g e r l g m r e s r l i g 彳r e 三、基于双正交化的方法 取l 。= k 。( 印，) ，k 。= k 。( a ，r o ) 显然，若a 是对称矩阵时，这一类方法就 是正交投影方法下面主要考虑是非对称矩阵的情形b i c g ( b i c o n j u g a t eg r a d i e n t ) 由 l a n c z o s 基于l 锄c z o s 算法提出，1 9 9 4 年f l e t c h e r 对b i c g 算法进行了改进，现在的许 多双正交化方法都可以由b i c g 方法推导出来，算法3 给出了b i c g 算法b i c g 方法 要使用系数矩阵a 的转置，在许多情形下会出现中断f 除数为零) 现象，同时在收敛过程 中残差范数会出现增大现象c g s ( c o n j u g a t eg r a d i e n ts q u a r e d ) 方法避免了使用a 的 转置计算，计算量与b i c g 方法相同由于其一步迭代相当于b i c g 方法的两步迭代，所 以从理论上讲其收敛和发散的速度是b i c g 方法的两倍当用b i c g 方法出现中断时， c g s 方法也出现中断 算法4 b i c g 方法 1 )c o m p u t er o = 6 一a z o ，c h o o s er js u c ht h a t ( r ；，r o ) o s e tp o = 7 o ，瑞= 喵 2 ) f b rj = o ，1 ，2 ，一，u n t i lc o n v e r g e n c e ，d o 3 ) a ，= 器 4 )+ 1 = z ，+ n j 功 8 r s s 日 e 兄 m m 一 哪 日 + 一 r r s c ，) s e e e 冗 r r m m m g g g b 0 f 第二章k r y l o v 子空间方法 取l 。= k 。= 匠。( 4 丁a ，a t r o ) ，这类方法是将c g 方法应用到法方程a t a z = a 丁6 和印a o = 6 上对大多数问题，印a 条斋甘甘锱茳警，薹蓄奇善品堪臻晰瑚判晤扮 期辩捱扯；薹誓藜摹l 釜磊掣! l m 要誊l = t i | | 斗羹耋誊e ! 型i 5 豢差! j ；皇 挺军j ；? 删萎攀；哥璧 带：样蠹量e 棼霪萎茧堑型型雾量喜享一霹鞴l 蒌鞲蓖的影响下面考虑利用一次全局通讯 完成所有内积计算来减少这种影响 i b i c r 算法设计的主要思想就是通过数学推导消除6 ) 和1 0 ) 步的内积计算具有的数 据相关性，使它们可以同时计算，并利用一次全局同步来完成所有内积计算 令k = a ，鲒= a 和：，骱= a p 。，o 。= ( r ：，k ) ，6 n = ( 簖，) ，考虑b i c r 算法中的 6 ) 和1 0 ) 步 第6) 步变为 驴揣黠= 黜= 毒 c 。z ， 第1o ) 步变为 风： 错= 等= 等 饵。， 下面考虑( 簖，) = ( a t 成却。) ，步3 ) 和步4 ) 的公式两边分别乘4 和4 r 并做内积 有 (瞻，) = ( a 丁r ：+ 风一1 a 丁璐a r n + 风一1 a 砌1 ) = ( 4 丁r 三，a r 。) + 风1 ( a t r ：，却。一1 ) ( 3 2 3 ) + 风一1 ( a 丁p ：一。，a h ) + 藤一1 ( a 一1 ，却。一1 ) 令 坛= a 丁惦，c 。= ( a t r ：，a ) = ( 坛，t 。) ，厶= ( a 7 蠊巾a ) = ( 酩一。，t 。) ， = ( 4 了r ：，4 一1 ) = ( ，q 忆一1 ) 第二章k r y l o v 子空间方法 = ( 一1 ，a t r ：一1 ) + 风一2 ( r 。一1 ，a 丁p ：一2 ) 一0 _ 。一1 ( 却。一l ，a t p ：一1 ) 利用( 224 ) 可得( ，a ：一1 ) 一( 7 。- 1 ) a t 璩一。) = o ，所以 同样可利用r i 一，和a m 做内积来得到 ( r ：_ 1 ) 却。) = ( r ：_ 1 ) a i 囊m 一_ ! 銎g 鬟一i ! i囊多t!皇倒!匿膨一i；孽二；委i威；饕铬；拍一li|g=a1h= ( 1 z 。1 + 岛。 。一1 一如q 。一1 ) = e 。一1 + 风q 。一1 一矗c 。一1 ， r 手d n = 曙a k= 曙a ( u 。一1 一a 。1 ) = r 手e 。一1 一o 。r 孑c 。= 如一1 一o 。， = r 手= r手( 一1 + 风一1 一如1 ) = 如一1 + 风h l 一如7 r 。一1 ， p 。= 曙t h1= r 彳a r 。l = r 孑4 ( s 。一1 一u 。一1 k 一1 ) = r 手“一1 一u 。一1 r 吾d 。一l ， = 咖。一l un 一1 ( 6 n 一2 一n n 一1 丌n1 ) = 一1 一u n 1 6 n 一2 一u n 一1 n 一1 7 r 一1 ， 民= r 手e n =r 手a “。= 曙a a = r 手a a ( s n 一亡n ) = 霄a ( 一如) 改进的b i c rs t a b 方法用算法的形式写出为 算法8 i b ic r s i a b 方法 1 )c o m p ut er o = 6 一a z o ，c h o o s er ；s u c ht h a t ( r ；，r o ) o ， “o = a r 0，e o = a “o ，0 = a t 喵，g o = t b = p o = c 0 = 0 ， 口一1 = 丌o = 0 = 0 ，妒o = r j t 札o ，印= 喵t e o ，伽= o o = 岫= 1 f ( i r 他= 0 ，1 ，2 ，u n t i lc o n v e r g e n c e ，d o p n = 咖n一1 一u n 一1 d k 一2 + u n 一1 0 n 一1 丌n 一1 矗= 急扎- l，风= 告 p n2r n一1 + 尻鼽。1 一矗 。一1 第二章k r y l o v 子空间方法 组合上述公式，便可得到b i c r 方法 算法5 b i c r 方法 c o m p u t en ) = 6 4 z o ，c h o o s er ；s u c ht h a t ( r j ，4 r o ) o s e tp 二l = p 一1 = 0 ，p 一1 = o f b rn = o ，1 ，2 ，一，u 1 1 t i lc o n v e r g e n c e ，d o p 。一+ 岛一1 p 。一l 蜣= r ：+ 风一1 峨一l a p 。= a + 风一1 a p h 一1 。= 若锚 ( 全局同步点) z n + 1 = z 十n p “ r n + 1 = r n o n a p n r 0 1 = r ：一q 。4 t 戚 风：譬# ；掣( 全局同步点) 一“一 ( r ：，a r n ) 2 匕h j lh j 二y7 、 e n d d o 从推导过程可以看出，当系数矩阵a 是对称时，且选择，；= r o ，那么b i c r 方法就 退化为c r 方法 表2 1 ：每次迭代的计算量对比 方法向量校正矩阵向量乘内积全局同步点 b i c r5 22 2 b i c g6222 在表2 1 中，比较了b i c r 算法与b i c g 算法在一个迭代步中的基本操作从表2 1 可以看出b i c r 算法与b i c g 算法的计算量几乎相同 2 3b i c r s t a b 方法 为了克服在b i c r 算法中使用a 7 的缺点，并增加数值稳定性，我们提出了b i c r s t a b 1 3 ” 劫勘回砷踟加n 第二章k r y l o v 子空间方法 算法其设计的主要思想来自于对b i c r 算法的分析在b i c r 算法中，第n 个迭代步的 残差向量可以记为 r 。= ( a ) r o ，( 2 3 1 ) 其中是次数为n 的特定多项式，满足关系式 西。( o ) = o ；西。+ l ( ) = ，。0 ) 一o 。t 7 h ( f ) ( 2 32 ) 类似地，共轭方向多项式记为 m = ( ) r o ，( 2 3 ，3 ) 其中丌n 是次数为n 的多项式，满足关系式 丌。+ 1 ( ) = 毋。+ l ( ) + 魔7 r 。0 ) ( 2 34 ) 从b i c r 算法可知，r i ，蜣与，m 具有相同的递推关系式，所以将( 2 3 1 ) 、( 2 3 3 ) 中 的4 用a r 替换就可以得到 壤= 。( a 丁) a 丁7 ；，p ：7 r 。( 4 t ) 4 t r ； ( 2 3 5 ) 另外，从b i c r 算法中可知参数。，。可记为 ( ( a ) r o ，咖。( a t ) a 7 喵)( 4 候( a ) r 0 ，r j ) ( 23 6 ) 从式子( 23 6 ) 我们分析出只要能找到两个向量来分别近似织( 4 ) r 0 和”：( a ) r o ，就可以避 免使用a 7 因此，我们把b i c r s t a b 方法中的残差向量定义为如下形式 r 。= 妒。( a ) 砂。( 4 ) r o ，( 23 7 ) 其中，( ) 是b i c r 算法的残差多项式，( ) 是为了保证b i c r s t a b 算法相对于b i c r 算法具有更好的稳定性和收敛光滑性而引入的一个满足短递推的多项式 妒。+ l ( t ) = ( 1 一u 。t ) 掣( t ) ，( 23 8 ) 其中是常数 1 4 第二章k r y l o v 子空间方法 为了便于推导我们先不考虑上述式子中的常数，由( 2 3 2 ) 和( 2 3 8 ) 可得 西。+ l + 1 = ( 1 一w 。) ( t ) 咖。+ 1 = ( 1 一u 。t ) ( 妒。九一a 。7 r 。) ( 2 3 9 ) 由( 2 3 9 ) 可知，如果知道妒。丌n ，则可得到+ + ，的递推关系式 而 丌n = ( + 风一1 7 r 。1 ) = 讥+ 阮一1 ( 1 一u 。一l ) 妒竹一1 _ l( 23 1 0 ) 我们定义： 乜。= a = a ( a ) 蛾a ) r o ，阮一呶( a ) ( a ) r o ，= 4 p n ，= a 由( 2 3 9 ) 、( 2 3 1 0 ) 可得 r 。+ 1 = ( j 一“。a ) ( f k q 。a f k ) ，( 2 3 1 1 ) 肼1 = r 。+ 1 + 风( 1 一u 。a ) p 。( 2 3 1 2 ) 下面我们来计算迭代过程中的常数a 。与风在算法b i c r 中风= 忍+ 瓦，其中 西= ( ( a ) r 0 ，九( a t ) a t r ；) = ( a ( 4 ) 2 r 0 ，r ；) 我们令 尸。= ( 砂。( a ) r o ，砂。( a 7 ) a 7 r ；) = ( a ( 4 ) 砂。( a ) ，r ；) = ( u 。，r ；) 为了推导出陬与磊之间的关系，我们将讥( a 丁) a t r j 展成幂级数可得 p 。：( 。( a ) m ，q i ”( a 7 ) “a 丁r j + q 5 “( a 丁) ”1 a t r j + ) 又因为当 t 。比较( 4 ，2 4 ) 和( 4 2 5 ) 式，我 们得到i b i c r 方法具有比b i c r 方法更好的并行性 由( 4 24 ) 和4 2 5 ) 式关于处理器台数p 极小化7 售只和z 。篙e r ，我们得两种方法具 有最小并行时间的处理器台数分别为 r = 譬掣= 罢铲z 吻煳5 1 i 再了2 可i 酝丁一 0 j 和 c r = 群= 堕掣( 4 。7 ) 所“c r2 可丁西丽一2 百砑- 卜6 。j 2 7 第四章并行化方法的实现和性能分析 由于对大规模分布式并行计算机， t 。 屯，从而甓等z2 + j i 岛 2 ( 对任意 n ： 0 ) 因此，i b i c r 方法比b i c r 方法具有更好的并行可扩展性 当固定，j p 充分大时，t 。 气，此时i b i c r 方法相对于b i c r 方法的改进比率 为 ”= 警。高筹藩筹杀一， 。s ，。 码：c r4 t 。p l o g p + ( 4 札：+ 1 4 ) 2 n 。 r 这说明对比b i c r 算法，i b i c r 算法的性能改进比率趋向于5 0 4 2 2i b i c r s t a b 与b i c r s t a b 的性能比较 同样从i b i c r s t a b 算法的推导过程也可以看出，i b i c r s t a b 算法和b i c r s t a b 算法数学上是等价的表4 2 给出b i c r s t a b 和i b i c r s t a b 算法每次迭代的计算量和 所需要的全局通讯次数比较情况 表4 2 ：每次迭代的计算量对比 从表4 2 可看出：为了消除2 个全局同步点，与b i c r s t a b 算法相比，i b i c r s t a b 算法增加了三次内积计算，四次向量校正由于k r y l o v 子空间方法的主要计算时间集中 在矩阵向量乘积和全局通信，相对于全局通讯的减少，这点代价是完全可以接受的，这也 可从数值试验中得到验证 由表4 2 可得b i c r s t a b 每次迭代所需要的时间为 赐脚a b 2 9 删+ 3 p r o d ( 1 ) + 2 一 触9 1 = ( 4 扎。+ 2 2 ) o ，i p + 6 l o g p ( t 。+ k ) + 4 n m 如+ 4 ( 2 礼5 + 咒m ) t 而i b i c r s t a b 每次迭代所需要的时间为 ，j 翟g r s t a b= 1 3 。“州+ t m 。母，础( 6 ) + 2 t m 。t 。 。 = ( 4 n ：+ 3 6 ) o 爿p + 2l o gp ( 。+ 6 t w ) + 4 扎m 。+ 4 ( 2 礼古+ 讫。；) t w f 4 2 1 0 ) 第四章并行化方法的实现和性能分析 并行算法的效率e 定义为加速比与处理器台数p 之比，所以 e = ； 瓦。( ) 一 e 。( ) p 。，( ，尸)正。口+ 瓦( ，p ) 1 耳i 羹嚣黼1 8 戮 麓戆豁巡跫甄；瑚瞒嚣赜澍嚣髓烈疆羿藩；鍪嚣饕羟鹌靛泛墨鬣驺甓i 矮嚣攫凌囊蔫 羹i ，。鬣篱露嚣冀羹。弱潜； ￥i l 挞，裂翡攀冀冀嚣嚣豁錾嚣塑浚秀襄疆净搿茎，囊 餮夔嚣 1 0 8 0 1 0 8 0 1 2 0 0 “1 2 0 0 1 3 2 0 1 3 2 0 1 4 4 0 1 44 0 1 5 0 0 x1 50 0 2 9 20 34 1 1 5 8 3 4 3o 681 9 7 7 4 2 81 593 7 2 3 5 3 52 354 3 9 6 6 0 933 35 4 7 1 6 6 84 036 0 3 3 76 64 8 56 3 3 5 1 0 7 07 6 77 1 6 6 1 2 4 09 4 27 5 9 7 1 29 988 46 8 0 7 1 6 0 21 23 47 7 0 3 1 5 5 2 1 20 47 7 5 6 2 9 5o 2o 6 7 7 3 2 50 ，561 7 1 5 4 1 41 172 8 2 1 4 6 5 1 51 3 2 3 7 5 4 92 ，093 8 1 2 6 1 42 393 8 ，8 5 第四章并行化方法的实现和性能分析 图2 ：i b i c r 方法的等效率曲线 同样我们利用等效率分析法来比较b i c r s r a b 与i b i c r s t a b 方法的可扩展性 从表4 2 我们可以得到在每一个迭代步中串行求解b i c r s t a b 所需要的时间为 嚆各雕丁a b = 9 篓墨d + 3 景嚣，。d ( 1 ) + 2 篇恶。= ( 4 几z + 2 2 ) o ，l ( 4 3 1 3 ) 每个迭代步b i c r s t a b 方法所需要的通讯时间为 丁茹淼t b = 3 t 嚣器。d ( 1 ) + 2 嚣羔。= 6 l o g p ( 如+ ) + 4 s + 4 ( 2 + ) ，( 4 3 1 4 ) 每个迭代步b i c r s t a b 方法所需要的通讯开销时间为 7 苫琵s t a 日= p 7 譬苫r s t a b z 鬣昌r s t b = p 瑶子昌t a b ( 4 3 1 5 ) 把方程( 4 3 1 3 ) ，( 4 3 1 5 ) 代入( 4 3 4 ) 可得 茗昌脚旷高；。t 。b 一口= 蒜肌心而蒜引。g p ( 4 s - 1 6 ) 从表4 2 我们也可以得到在每一个迭代步中串行求解i b i c r s t a b 所需要的时间为 丁篙c 咒s t a 日= 1 3 t i ：嚣p d + 6 t 嚣2 。d ( 1 ) + 2 t i 翟己。一( 4 n ：+ 3 6 ) o ，i ( 4 3 1 7 ) 3 2 第四章并行化方法的实现和性能分析 每个迭代步i b i c r s t a b 方法所需要的通讯时间为 丁器：品s 丁 口= 6 i ：嚣甜( 1 ) + 2 f 篙：。= 2 i o g p ( 如+ 6 ) + 4 n m t 。+ 4 ( 2 礼b + 几m ) 埘，( 4 3 1 8 ) 每个迭代步i b i c r s t a b 方法所需要的通讯开销时间为 瑶弓r s 丁a 日= p 7 菘r s t a b 一：e 翟c r s 丁n 日= p t j 苫；苫盖s t 且日 ( 4 3 1 9 ) 把方程( 4 3 1 7 ) ，( 4 3 1 9 ) 代入( 4 3 4 ) 可得 瑶g 脚旷高磁赫， a 。b 。c r s t a b = i i t i ：；豢p l 。g p “i i i 丢i ：j p l 。g 尸( a 。z 。) 其中了。篇c r s t a b 和嚆昌脚t a b 分别是i b i c r s t a b 和b i c r s t a b 方法的墙上时 间，而巧暑努r s 刚b 和赡嚣裔r a 。，巧鬣弓r s 7 4 b 和昭i 孙s r 。b 分别是它们的通讯时间和冗 余时间在图3 ，4 中我们取不同的效率e ，利用( 4 3 1 7 ) 和( 4 3 2 0 ) 两个函数式分别画出 了i b i c r s t a b 和b i c r s t a b 方法的等效率曲线我们由方程( 4 3 2 0 ) ，( 4 3 1 6 ) 和图3 ，4 都可以看出i b i c r s t a b 比b i c r s t a b 方法的可扩展性好同时我们得出两种方法都满 足关系式一0 ( 尸l o g 尸) ， 图3 ：b i c r s t a b 方法的等效率曲线 第凹章并行化方法的实现和性能分析 图4 ： i b i c r s7 r a b 方法的等效率曲线 第五章数值实验 另外，我们也比较了b i c r 与b i c g 、i b i c r 与i b i c g 1 4 方法并行计算的收敛性， 得到了与 7 中串行计算相类似的结果，见图7 ，8 图中可以见到， b i c r 和i b i c r 方法 具有相同的收敛速度，且分别比相应的b i c g 和i b i c g 收敛更快，i b i c r 方法还克服了 i b i c g 方法出现的残差范数的振荡现象 图7 ：b i c g 与b i c r 方法在1 6 台处理器上的收敛结果 图8 ：i b i c g 与i b i c r 方法在1 6 台处理器上的收敛结果 第五章数值实验 表5 2 ：性能比较 c p u# h n 妇o w d sb i c r s t a bi b i c r s t a b i m p r o v e r r 忙n t s t bt b l cc bt i bt f b cc i b 叩靶 2 21 2 0 1 2 01 9 2 0 5 4 2 8 1 1 1 7 9 0 1 7 9 3 6 6 7 7 6 8 5 2 4 42 4 0 2 4 04 4 21 5 43 4 8 3 2 4 0 0 5 6 2 3 2 4 4 5 7 0 6 3 6 4 8 84 8 0 4 8 0 4 4 428 3 6 3 8 32 8 51 0 13 5 3 13 5 8 16 4 3 1 1 0 l o6 0 0 6 0 0 4 9 0 3 1 66 4 3 53 1 21 2 64 0 4 43 6 3 36 0 1 3 1 2 1 27 2 0 x7 2 05 4 83 ，7 06 7 6 5 3 3 8 1 5 44 5 5 63 8 3 2 5 8 3 8 1 4 1 48 4 0 8 4 06 9 85 3 07 5 9 83 6 51 8 04 9 3 24 7 7 l6 6 0 4 1 6 1 69 6 0 9 6 06 2 044 17 1 1 53 9 62 1 15 3 4 13 6 1 35 2 1 5 5 2i b i c r s t a b 与b i c r s t a b 我们用五点差分对( 5 1 ) 进行离散，固定每台处理器求解3 6 0 0 个未知量，进行2 0 0 0 步 迭代，我们比较了i b i c r s t a b 和b i c r s t a b 方法的总墙上时间性能和通讯时间性能， 所得结果列在表5 2 之中，其中码和乃b 分别为b i c r s7 r a b 和i b i c r s t a b 方法的 并行求解时间，而r 和乃b g 分别为两种方法的通讯时间，c 名= ( 珏r ) 1 0 0 和 d b = ( 乃b _ g d b ) 1 0 0 分别是两种方法的通讯时间占墙上时间的百分比改进( ) 是 用分式q = ( 1 一乃口乃) 1 0 0 和q g = ( 1 一乃日一g 一e ) 1 0 0 计算得到的表5 2 表 明：随着处理器台数的增加，改进是明显的随着处理器台数的增加，内积计算的并行瓶 颈越来越明显，因此降低全局通讯次数所获得的改进就越来越显著 由表52 我们也看到，在通讯时间方面， i b i c r s7 r a b 方法比b i c r s7 r a b 方法的改 进趋向于6 6 7 ，这与我们的理论分析是相吻合的并且b i c r s t a b 和i b i c r s t a b 分 别在6 4 和2 5 6 台处理器时通讯占优性能改进结果也见图9 可见i b i c r s t a b 方法比 b i c r s t a b 方法具有更好的可扩展性 第五章数值实验 图9 ：i b i c r s t a b 方法相对于b i c r s t a b 方法的性能改进 同时，我们固定问题规模= 9 6 0 9 6 0 ，进行3 0 0 0 次迭代分别检测了b i c r s t a b 和i b i c r s t a b 方法的加速比，我们是利用在一台处理器的串行时间比上在p 台处理器 上的并行时间来得到加速比。在图1 0 中，l i n e a r 代表线性加速比， b i c r s r a b 代表 b i c r s t a b 方法的加速比，i b i c r s t a b 代表i b i c r s t a b 方法的加速比由图1 0 我们 也可以看出i b i c r s t a b 方法比b i c r s t a b 方法具有更好的可扩展性同时，也可以看 到，两种并行方法都出现了超线性加速比，其原因是我们的并行方法在单处理机时并非最 优的 图1 0 ：i b i c r s t a b 与b i c r s t a b 方法的加速比对照 第五章数值实验 对( 51 ) 的九点差分离散所得的线性方程组的试验结果与表5 2 类似 另外，我们也比较了i b i c r s r a b 与b i c r s t a b 、i b i c g s t a b 【1 3 方法并行计算 的收敛性，得到了与串行计算相类似的结果，见图1 1 ，1 2 图中可以见到，i b i c r s t a b 方 法比b i c r s t a b 、i b i c g s t a b 方法收敛的更快 图1 1 ：i b i c r s t a b 与b i c r s t a b 方法在1 6 台处理器上的收敛结果 图1 2 ：i b i c r s t a b 与i b i c g s t a b 方法在1 6 台处理器上的收敛结果 4 1 第六章总结与展望 本文提出了两种适合于分布式并行环境的，求解大型稀疏非对称线性方程组的并行迭 代法即改进的双共轭残差( b i c o n j u g a t er e s i d u a i ) 方法和改进的双共轭残差稳定化( b i c o n j u g a t er e s i d u a ls t a b i l i z e d ) 方法，分男9 简记为i b i c r ( i m p r o v e db i - c o l l j u g a t er e s i d u a l ) 方和i b i c r s t a b ( i m p r o v e db i c 。n j u g a t er e s i d u a ls t a b n i z e d ) 方法作为解决k r y l o v 子 空间方法并行化瓶颈的方法之一，我们通过算法重构，有效地将b i c r 和b i c r s t a b 方法 中所需要的两