petsc求解二维bratu方程.pdf

1、PETSc 求解二维 Bratu 方程 * 程强 1 张海斌 2 曹建文 3 中国科学院计算机网络信息中心超级计算中心， 100080 北京工业大学应用数理学院， 1000 22 中国科学院软件研究所并行软件研究与开发中心， 100080 摘要 本文以求解二维 Bratu 方程为例，详细介绍了使用 PETSc 求解非线性方程组的方法和实现。在固定问题计算规模下，对程序的并行可扩展能力和浮点性能进行了测试和分析，最后从 KSP 迭代方法、预条件子选择、内层迭代次数和雅可比矩阵求解方法等方面讨论了最优运行参数的选择方法。 关键 词 非线性方程求解； 自动微分；浮点性能 SNES 是可移植可扩展科学

2、计算工具箱 1（ PETSc）的三个核心解法器之一，用于求解大规模非线性方程组。它采用牛顿迭代方法，而每个牛顿迭代步则使用 Krylov子空间迭代方法 2和不同预条件子的组合实现，这非常适合于高效求解一大类 PDE 方程离散化后形成的大型稀疏非线性方程组。 用户使用 SNES 求解一个非线性方程组时，仅仅需要提供非线性函数本身和可选雅可比矩阵的计算程序。通常，雅可比矩阵可采用有限差分由 PETSc 提供的函数自动实现；但是，当网格规模过大时，雅可比 矩阵的计算代价将非常昂贵，并且计算精度受截断误差的影响。 自动微分方法 3（ AD ）在无截断误差意义上分析求解函数导数，具有计算代价小、精度高和

3、可靠性好等优点。特别地，运用 AD 方法计算函数梯度的计算时间代价一般不超过其本身代价的五倍 4 ，并且与独立变元的数目无关。目前， PETSc 使用自动微分软件（ ADIC/ADIFOR）来提供求解雅可比矩阵的源代码。 对于 Bratu 方程的求解，除了雅可比矩阵的计算精度和代价对 SNES 的整个运行代价影响较大之外，在每个牛顿迭代步内如何选择合适的 Krylov 子空间迭代方法和不 同预条件子的组合对整个计算性能的影响也不可能被忽略。此外，为了追求尽可能低的计算代价，在固定问题规模和计算精度约束下，我们往往需要还在外层牛顿迭代步数和内层 KSP 迭代精度（或步数）之间寻求一个合适的折衷。

4、幸运的是， PETSc 为实现和调试最优性能提供了一个理想的试验平台。 本文受到国家自然科学基金项目“自动微分方法研究及其实现”（ 60503031）和十五信息化建设专项“超级计算环境建设与应用” (INF105-SCE)的资助。 E-mail: , Tel: 010-5881-2132 1 问题的描述 我们考虑区域 )1,0()1,0( = 上的二维 Buratu方程 =),(,0),(,02yxuyxeu ul(1) 采用均匀网格及五点中心差分格式将（ 1）离散化后，就得到一个标准形式的非线性方程组 0)( =uF 。 PETSc 提供了两种优化方法可供选择，即一维线性搜索和信赖域方法。而

5、在每个牛顿迭代步内， SNES 都要近似求解一个线性方程组来得到了一个下降或搜索方向，因此用户需要设置不同 Krylov 子空间迭代方法和不同预条件子的类型，以及内层 KSP 迭代的最大迭代步数或计算精度。 2 固定计算规模的并行可扩展性能 为了考察 PETSc 程序的并行可扩展能力，我们固定总迭代步数，并使用一维线性搜 索方法求解 Bratu 方程 。其中内层 KSP 迭代方法采用共轭梯度法（ CG）和块雅可比预条件子组合，固定牛顿迭代步数为 8 或 32，而每个牛顿迭代步中线性方程组求解的迭代步数为 32。迭代误差取 NORM2 范数。 表 1给出了在两种不同的网格规模下对 SNES 范例

6、 EX5 在深腾 6800 机群系统上的部分测试结果，包括函数范数、计算时间（秒）、单机浮点性能（ Mflops ）、总体浮点性能（ Mflops）和通信次数。 表 1 线性搜索方法在深腾 6800 系统上固定迭代步数的测试结果 处理器数 1 2 4 8 16 32 64 函数范数 5.0e-5 3.6e-4 1.2e-6 6.2e-5 7.90e-5 8.5e-4 6.0e-4 计算时间 201.7 99.32 44.56 19.27 7.044 2.956 1.568 平均性能 42.53 43.20 48.01 55.42 75.90 90.40 84.96 浮点性能 42.53 86.

7、39 192.0 443.4 1214 2893 5438 512512，3232通信次数 0 2.302e3 9.176e3 2.288e4 5.510e4 1.198e5 2576e5 函数范数 2.9e-3 3.9e-3 4.9e-3 4.0e-3 3.9e-3 5.2e-3 4.4e-3 计算时间 838.90 444.3 210.3 98.92 46.80 19.88 7.620 平均性能 41.04 38.79 40.11 43.56 46.00 53.26 70.80 浮点性能 41.04 77.59 160.4 348.5 736.0 1704 4531 10241024，83

8、2通信次数 0 2.318e3 8.992e3 2.320e4 5.563e4 1.169e5 2.616e5 1 2 4 8 16 32 6432641282565121024204840968192Total Performance (Mflops)Number of ProcessorsLinear Growth512x512,32x32Linear Growth1024x1024,8x32图 1 线性搜索方法在深腾 6800 上固定迭代步数的整体浮点性能测试果 图 1 显示了两种不同规模在深腾 6800 上固定迭代步数后的整体浮点计算性能随处理器数变化的结果。不难看出，二者都获得了超线

9、性的加速性能，并且随着处理器数目的增加越来越明显，这是局部内存需求降低后引起单机平均浮点性能的提高所致；同时，较小计算网格规模则获得了更为明显的超线性加速现象，但随着处理器数目的增加其平均性能逐渐变差 。 3 最小时间代价驱动下的最优参数选择 在给定计算精度约束下，我们来讨论如何应用 PETSc 算法资源获得尽可能小的计算时间代价。 就求解二维 Bratu 方程而言， 一维线性搜索方法比信赖域方法在计算稳定性方面要优越得多。如果内层 KSP 迭代步数太小，不仅将降低 SLES 环境设置和启动带来的额外开销的复用率，而且使得每个牛顿迭代的计算误差过大，信赖域方法因而表现出极慢的收敛性能，许多计算

10、甚至中途终止。在单机浮点性能方面，两种方法大致相当。此外，两种最优化方法在数据存储开销方面的差别并不大。 在以下讨论中，我们选择 一维线性搜索方法，处理器数目取 8（默认值）。除非特别指出，计算规模取 512128（默认值），最大相对误差精度为 1.0e-6（默认值），内层迭代方法默认为共轭梯度方法和块雅可比预条件子的组合。另外，固定最大牛顿迭代步数为 104。 3.1 KSP 内层迭代步数的选择 表 2 给出了一维线性搜索方法取不同 KSP 迭代步数在曙光 2000 系统上固定计算精度的测试结果，包括内层迭代步数、外层迭代步数、计算时间（秒）、函数范数（ 10-6）和单机浮点性能（ Mflo

11、ps）。不难看出，整个计算代价对于不同 KSP 迭代步数是非常 敏感的。当 KSP 迭代步数比较小时，由于每个牛顿迭代步求解线性方程组的精度太低，导致整个求解过程收敛很慢。随着内层 KSP 迭代步数的增加，线性方程组的求解精度逐渐提高，加速了牛顿迭代的收敛过程。但是，当牛顿迭代步数减小到一定程度时，仅靠增大 KSP迭代步数毫无意义。 表 2 一维线性搜索方法在不同内层迭代步数下固定精度的测试结果 内迭代步数 2 4 8 16 32 64 128 256 外迭代次数 58 158 153 109 15 7 3 3 计算时间 11.97 42.06 59.67 68.04 17.31 14.74

12、12.29 24.20 函数范数 9994 3232 1.5 1.5 0.55 0.91 0.005 0.0002 平均性能 5.168 5.926 6.505 6.925 7.176 7.645 7.757 7.823 3.2 雅可比矩阵的求解 PETSc 为雅可比矩阵求解提供了三种方式：有限差分、自动微分和“无矩阵”方法。在下面的测试中，我们将看到不同方法对求解性能的影响程度。 这里内层迭代方法均选共轭梯度方法和雅可比预条件子的组合，最大牛顿迭代步数为1024，内层 KSP 迭 代步数为 32。对于“无矩阵”方法， PETSc 没有提供任何预条件子，用户需要自己提供雅可比预条件子。表 3

13、给出了三种雅可比求解方法在曙光 2000 系统上的对比测试结果 ， 包括函数范数 ( 10-6 )、计算时间、迭代步数、单机浮点性能和存储开销。 从图 2 不难看出，三种计算方法随处理器数目增长呈现出并不规则的计算时间代价。当处理器数不超过 8 时，有限差分方法比其它两种方法具有更高的时间代价，除了处理器数为 2时自动微分方法最为突出外，其它两种方法大致相当。当处理器数目超过 8 后，三种方法的时间代价曲线就变得更加不规则了， 总的趋势是“无矩阵”方法在时间代价上的优势十分明显。 在计算精度方面，自动微分方法具有非常明显的有优势，处理器数目为 4 和 32 时超线性收敛速度甚为明显。这同自动微

14、分方法具有高精度方面的理论优势是一致的。另外，在存储需求方面，自动微分方法也是最小的。 表 3 三种雅可比矩阵求解方法在曙光 2000 系统上的对比测试结果 处理器数 1 2 4 8 16 32 函数范数 1.48 1.43 1.45 1.39 1.49 0.07 计算时间 1695 828.1 447.5 220.0 127.3 88.66 牛顿迭代 340 334 366 244 226 268 平均性能 10.82 10.90 11.05 7.483 6.003 4.698 有限差分方法存储需求 27.3 20.6 10.3 5.32 2.74 1.33 函数范数 1.42 0.19 0

15、.012 1.40 1.47 0.07 计算时间 1724 556.7 346.3 190.9 108.8 88.85 牛顿迭代 328 268 330 249 224 268 平均性能 9.455 11.96 11.86 8.124 6.406 4.689 自动微分方法存储需求 23.3 17.2 8.61 4.40 2.38 1.31 函数范数 5.60 2.82 3.65 1.30 12.4 0.724 计算时间 1596 707.8 353.5 172.3 154.6 65.34 牛顿迭代 282 282 282 227 281 141 平均性能 13.38 15.10 15.11 1

16、2.48 8.596 5.107 “无矩阵”方法存储需求 23.8 17.4 8.87 4.60 2.42 1.32 1 2 4 8 16 321001000WallTime (second)Number of ProcessorsFinite DifferenceAutomatic DifferentiationMatrix Free Method图 2 三种计算微分方法求解雅可比矩阵的时间代价 3.3 Krylov 子空间迭代方法和预条件子的选择 在每个牛顿迭代步内， SNES 都要近似求解一个线性方程组，因此 KSP 求解在整个计算中占有很大比重，如何选择 Krylov 子空间迭代方法和

17、预条件子对于提高性能具有重要影响。这里选择几种有效的迭代方法和并行预条件子，在曙光 2000 系统的 8 个处理器上做如下组合测试：固定内层 KSP迭代步数为 32，雅可比矩阵求解采用默认的有限差分方法。 表 4 给出了 不同 Krylov 子空间迭代方法和不同预条件子组合测试结果，包括函数范数（ 10-6 ）、计算时间、牛顿迭代步数和单机浮点性能（ Mflops）。不难看出，在固定的计算精度下，双共轭梯度平方法（ BCGS）和加性 Schwarz 预条件子（ ASM）的组合在二维 Bratu方程求解中具有最小的计算时间代价和最少的收敛次数，而双共轭梯度平方法与块雅可比预条件子的组合则稍差一些

18、，但计算精度最好。在平均性能方面，它们均比较低，仅有广义最小残差法（ GMRES）和雅可比预条件子组合的 50%左右。 我们注意到，经典 Gram-Schmidt 方法（ CGS）、无转置拟最小残量法（ TFQMR）以及共轭梯度方法（ CG）与加性 Schwarz 预条件子（ ASM）组合在中途终止了收敛过程。这是因为内层 KSP 线性方程求解的精度较低所致，只要适当增大内层 KSP 迭代次数就可以保证计算收敛到固定计算精度。在很多情形下，牛顿迭代次数随处理器数目的不同而存在着较大的差异。另外，在迭代方法和预条件子的选择方面，用户还须特别考察几种较小时间代价数值方法的并行可扩展能力。 表 4

19、不同 Krylov子空间迭代方法和不同预条件子组合测试结果 Krylov子空间方法 函数范数 计算时间 牛顿迭 代 平均性能 cg bicg gmres bcgs cgs tfqmr cr 1.40 1.46 1.50 1.40 49677 17442 1.51 186.9 219.5 469.6 129.9 15.28 898.0 307.6 249 222 408 117 14 793 471 jacobi 8.297 10.23 17.11 10.24 10.00 10.73 10.17 0.55 0.55 1.49 0.07 367 502 1.14 18.08 23.74 19.17

20、 12.76 29.47 234.4 15.17 15 15 11 7 17 106 16 bjacobi 6.828 9.045 12.62 8.545 8.749 7.467 9.347 181.1 0.32 0.97 0.35 319 4329 285 20.41 46.13 22.11 10.22 20.42 6.676 16.37 17 22 11 5 10 3 13 预 条 件 子 asm 7.113 6.909 10.97 7.691 7.525 7.826 7.265 3.4 最优参数测试结果 在以上讨论的基础上，我们继续讨论求解二维 Bratu 方程 在固定计算规模和计算精度

21、下的最小时间代价，并与 PETSc 默认参数选择的结果进行对比。这里网格规模取 512512，计算精度取相对误差 1.0e-8（ snes_rtol）。在深腾 6800 机群系统上，我们选择一维线性搜索方法进行求解，其中内层线性方程求解取双共轭梯度平方法（ BCGS）和加性 Schwarz 预条件子（ ASM）组合，而最优 KSP 内层迭代步数随处理器数目的不同而存在着较大差异（见表4）。另外，这里采用自动微分方法来求解雅可比矩阵。对于 选择默认参数的情形，为了提高牛顿迭代的效率，我们固定最大 KSP 迭代次数为 512。 表 5 给出了两种不同参数选择下的测试结果，包括函数范数（ 10-8

22、）、计算时间（秒）、牛顿迭代步数、内层 KSP 迭代步数、单机平均浮点性能（ Mflops）、时间加速比和并行效率。不难看出，选择默认参数的计算时间大致为选择最优参数的 7-9 倍，而计算精度与后者大致相同。但在单机平均浮点性能方面，选择默认参数明显比选择最优参数得到的结果高出近一倍。另外，在时间加速比和并行效率方面，当处理器数不超过 32 时选择最优参数一般比选择默认参数略 低一些，而达到 32 后情况则完全相反，并且二者随着处理器数目的增加均达到了超线性加速比（见图 3）。 表 5 深腾 6800 机群系统上选择最优参数后得到的对比测试结果 处理器数 1 2 4 8 16 32 函数范数

23、0.66 0.12 0.40 0.59 1.01 0.91 计算时间 1130 496.1 228.7 111.4 51.22 22.24 牛顿迭代 8 7 7 7 7 8 KSP 迭代 512 512 512 512 512 512 平均性能 81.37 81.12 87.97 90.25 97.99 129.8 加 速 比 1 2.28 4.94 10.14 22.06 50.81 默认参数选择效率 (%) 100 114 124 127 138 159 函数范数 0.28 0.51 0.54 0.71 0.19 0.086 计算时间 142.1 72.62 32.90 14.47 6.8

24、0 2.50 牛顿迭代 4 4 4 4 5 5 KSP 迭代 95 95 100 95 85 90 平均性能 41.00 40.15 46.62 50.51 60.25 86.93 加 速 比 1 1.96 4.32 9.82 20.90 56.84 最优参数选择效率 (%) 100 97.8 108 123 131 178 相对加速比 7.95 6.83 6.95 7.70 7.53 8.90 1 2 4 8 16 321248163264SpeedupNumber of ProcessorsReserved running parametersOptmized running parame

25、tersLinear Speedup图 3 在深腾 6800 上选取默认参数和最优参数的并行加速比曲线 4 结论 最优化方法的选择、 KSP 迭代步数选取、雅可比矩阵的求解方法和 Krylov 子空间迭代方法和预条子组合的选择都是影响二维 Bratu 方程 求解性能的重要方面。 PETSc 是一个优秀的软件包，它为用户在这些方面的选择和权衡提供了丰富的资源和试验平台，以追求最小的计算时间代价。在很多情形下，用户还必须在计算时间、存储开销和计算精度等方面做合理权衡。 对于复杂的大规模科学与工程计算而言，自动微分方法相比起有限差分在计算时间方面具有无比的优越性。事实上， SNES 求解效率的关键因

26、素之一就是如何高效地求解大型稀疏雅可比矩阵。 参考文献 1 Satish Balay, et. al PETSc 2.0 users manual. Technical Report ANL-95/11-Revision 2.1.3, Argonne National Laboratory, May 2003 2 戴华，求解大规模矩阵问题的 Krylov 子空间迭代方法，南京航空航天大学学报，第 33 卷第 2 期， 2001 3 Andreas Griewank, Evaluating Derivatives, SIAM, 2000 4 Cheng Qiang, Zhang Linbo and Wang Bin. Model Adjointisation and Its Costs, Science in China (F), 2004 作者简介： 程强， 1971年，男，四川邻水人，副研，博士，主要研究领域为自动微分方法、并行计算。