并行可扩展科学计算工具箱PETSc简介与应用

莫 则 尧

（北京应用物理与计算数学研究所）

一、 PETSc的起源与现状 二、 PETSc的成功应用典范 三、 PETSc的体系结构 四、 PETSc的核心组件 五、 PETSc程序示例 六、 PETSc的具体应用与比较 七、 PETSc的优点与缺陷

1

一、PETSc的起源与现状

1． 全称：

 并行可扩展科学计算工具箱

 Parallel Extensible Toolkits for Scientific Computing  http://www.mcs.anl.gov/petsc

2．起源：

 美国能源部ODE2000支持的20多个ACTS（Advanved Computing Test & Simulation program, 美国能源部超级计算中心http://www.nersc.com/）工具箱之一，其中包括：  能提供算法的工具：

 Aztec ：分布式存储并行机求解大规模稀疏线性代数系统库；

 Hypre ：线性系统求解预条件库；  Opt++ ：串行非线性优化问题数值库；

 PETSc ：并行可扩展科学计算工具箱，提供大量面向对象的并行代数数据结构、解法器和相关辅助部件，适合并行可扩展求解PDE方程（有限差分、有限元、有限体离散的隐式和显示格式）；  PVODE ：并行常微分方程求解库；

2

 ScaLAPACK ：

 SuperLU ：代数系统直接求解库；  算法开发辅助工具：

 Global Arrays ：以共享存储并行程序设计风格简化分布存储并行机上程序设计；  Overture ：网格生成辅助工具；

 POET ：并行面向对象环境与工具箱；  POOMA ：并行面向对象方法与应用，提供大量适合有限差分和粒子类模拟方法的数据并行程序设计（HPF）的C++类；

 运行调试与支持工具：CUMULVS，Globus，PAWS，SILOON，TAU，Tulip；  软件开发工具：

 ATLAS & PHiPAC ：针对当代计算机体系结构的高性能特征，自动产生优化的数值软件，可与手工调试的BLAS库相比较（？）；  Nexus , PADRE, PETE； 3．现状

 时间 ：1995年—现在；

 目前版本：PETSc-2.0.28 + patch，源代码公开（不包

3

含用户自己加入的核心计算子程序）；

 核心人员：数学与计算机部，Argonne国家重点实验室，Satish Balay, William Gropp, Lois C.McInnes, Barry Smith；

 参研人员：相关访问学者（几十人次，不同组件实现）；  可移植性：CRAY T3D，T3E，Origin 2000, IBM SP, HP UX, ASCI Red, Blue Mountain, NOWs，LINUX，ALPHA等；

 目前，已下载上百套；

4

二、PETSc的成功应用典范

1．PETSc-FUN3D：

 参考：W.K.Anderson etc., Achieving high sustained performance in an unstructured mesh CFD applications, SC’99. (http://www.sc99.org)

 FUN3D：四面体三维无结构网格离散、求解可压或不可压Euler和Navire-Stokes方程、串行程序、百万量级的非结构网格点， NASA Langley 研究中心W.K.Anderson开发，适合飞行器、汽车和水下工具的设计优化；  核心算法：非线性方程拟时间步逼近定常解、隐格式离散、Newton线性化、Krylov子空间迭代算法、加性Schwarz预条件（每个子区域近似精确求解），具有很好的数值可扩展性（即非线性迭代次数不随处理机个数的增加而显著增加）；

 移植周期：五个月（初步1996.10—1997.3），包括熟悉FUN3D与网格预处理器、学习ParMetis无结构网格剖分工具并集成到PETSc中、加入和测试PETSc的新功能、优化FUN3D面向向量机的代码段到面向cache的代码段、PETSc移植（非常少的时间，小于20天），并行I/O与后处理阶段还没完成；

5

 并行性能：

 代码行从14400减少为3300行（77%），包含I/O；  优化后，串行程序发挥各微处理器峰值性能的78%-26%；（附页1）

 ONERA M6 Wing, 2.8百万个网格单元（11百万个未知量），512—3072个 ASCI Red 节点（双Pentium Pro 333MHz，每节点一个进程），保持95%以上的并行效率，发挥峰值性能的22.48%；（附页2）  其他并行机：CRAY T3E、Origin 2000、IBM SP；  奖励：SC’99 Gordon Bell最佳应用奖；

2．石油：21世纪新一代油藏数值模拟框架； （USA Texas 大学油藏数值模拟中心）

3．空气动力学数值模拟中多模型多区域耦合流场问题：（USA自然科学交叉学科重点项目）；

4．天体物理中恒星热核爆炸问题数值模拟；（USA Chicago大学）

6

三、PETSc的体系结构

PETSc层次 BLAS LAPACK MPI KSP (Krylov子空间方法) PC (预条件) DRAW SNES (无约束优化、非线性解法器) SLES (线性方程解法器) PDE解法器 TS时间步Matrices Vectors Index Sets 7

四、PETSc的核心组件

1．程序设计技术

 面向对象程序设计风格 + 标准C语言实现；  标准C语言（C++）和FORTRAN语言接口；  强调以面向对象的数据结构为中心设计数值库软件，并组织科学数值计算程序的开发；  PETSc应用：

a) 根据应用需求，通过调用PETSc子程序建立数据结构（如向量、规则网格阵列、矩阵等）；

b) 调用PETSc各功能部件（如索引、排序、基于规则网格的拟边界数据分布与收集、线性解法器、非线性解法器、常微分方程解法器、简单的图形输出等）的子程序来对数据对象进行处理，从而获取PETSc提供的科学计算功能；

2．核心组件（component ，class）：

 向量（Vectors）：

 创建、复制和释放指定长度和自由度的串行或MPI并行向量（数据段被自动分配到不同的进程）；  向量操作：元素的赋值与索引、类似于BLAS的向量运算、向量的可视化显示等；

8

 索引与排序（Index Sets，Ordering）：

 向量和矩阵元素的局部与全局、自然与多色序号的对应关系；

 建立和释放任意两个集合的元素之间的对应和映射关系；

 适合无结构网格在进程间的的任意网格剖分，以及通过数据映射操作，完成相应无规则网格拟边界数据交换的消息传递；

 分布阵列（DA：Distributed Array）：

 建立在规则网格之上，一维、二维和三维（i=1,...,L, j=1,...,N, k=1,...,N），自动或指定阵列在进程间的区域划分，并沿拟边界设置宽度任意（离散格式需求）的影象（ghost）数组，存储邻近进程在相应位置包含的网格点的数值；

元素索引点

9

 阵列元素可包含多个自由度，且可以任意索引和访问，属于向量的一种特殊情形，多有向量操作均适合于它；

 数值计算中应用最为广泛的数据结构；  局部序和全局序可以索引和映射；  矩阵（Matrices）：

 指定维数和自由度大小的串行和MPI并行矩阵的生成、复制和释放；  矩阵元素的索引与访问；

 稀疏矩阵压缩存储格式：AIJ稀疏行、BAIJ块稀疏行、Bdiag块对角；  稠密矩阵；

 类似BLAS的基本矩阵操作，以及矩阵元素的标准或可视化输出；

 矩阵的隐式形成与使用；

 基于无结构网格划分工具（ParMetis）的并行矩阵的形成与使用；

 线性代数方程解法器（SLES）：

 基于稀疏矩阵与向量数据结构；  SLES的建立、访问、设置和释放；

 目前实现的解法器：Krylov子空间方法（GMRES、

10

CG、CGS、BiCGSTAB、TFQMR、Richardson、Chebychev），预条件（Additive Schwarz、Block Jacobi（ILU）、Jacobi、serial ILU、serial ICC、serial LU）；  收敛性测试与监控（实时图形显示迭代误差下降趋势）；

 非线性代数方程与无约束优化方程解法器（SNES）：

 基于稀疏矩阵、向量和SLES数据结构；  SNES的建立、访问、设置和释放；  Newton线性化：line serach、trust region；  收敛性测试与监控（实时图形显示迭代误差下降趋势）；

 PDE或ODE时间依赖方程解法器（TS）：

 基于稀疏矩阵、向量、SLES和SNES数据结构；  TS的建立、访问、设置和释放；

 方法：Euler、Backward Euler、拟时间步逼近定常解等；

 对象的打印、图形和可视化输出：  选项数据库支持：

 对所有用户指定的算法和功能部件的性能监控，可在MPI程序运行时由命令行参数输入，非常方便；

mpirun –np 4 example -ksp_type bcgs –ksp_xmonitor

11

 并行性能自动统计、输出（— log_summary）；  用户自定义选项（如网格规模、进程个数、图形可视化输出等）；

3．与其他库软件的功能互用与接口：

 BlockSolve95（并行ICC（0）、ILU（0）预条件）；  ESSL（IBM 快速稀疏矩阵LU分解）；  Matlab（数据的图形和数值后处理）；  ParMeTis（并行无结构网格图剖分）；  PVODE（并行常微分积分）；  SPAI（并行稀疏近似逆预条件）；

 SAMRAI，Overture（并行网格管理软件包）；

12

五、PETSc示例

1． 例一：（petsc-2.0.28/src/sles/examples/tutorial/ex2f.F）

! 求解二维规则区域上Dirichlet问题，其中调用PETSc的SLES部件求解 ！ 有限叉分离散所得的稀疏线性代数方程组。 !

! 程序运行: mpirun -np ex2f [-help] [all PETSc options] !

program main implicit none ! 头文件

#include \"include/finclude/petsc.h\" ！base PETSc routines #include \"include/finclude/vec.h\" ！vectors #include \"include/finclude/mat.h\" ！matrices

#include \"include/finclude/pc.h\" ！preconditioners

#include \"include/finclude/ksp.h\" ！Krylov subspace methods #include \"include/finclude/sles.h\" ！SLES interface #include \"include/finclude/sys.h\" ！system routines #include \"include/finclude/viewer.h\" ！viewer routines ！变量说明

double precision norm

integer i, j, II, JJ, ierr, m, n

integer rank, size, its, Istart, Iend, flg Scalar v, one, neg_one ！标量

Vec x, b, u ！近似解、右端项、精确解 Mat A ！线性系统系数矩阵对象 SLES sles ！线性解法器对象 KSP ksp ！迭代KSP方法 PetscRandom rctx ！随机数对象

! Beginning of program ！ 进入PETSc环境

call PetscInitialize(PETSC_NULL_CHARACTER,ierr) ！

m = 3 ！ 省缺X方向网格点个数 n = 3 ！ 省缺Y方向网格点个数 one = 1.0 neg_one = -1.0 ！

！从命令行输入X、Y方向网格点个数

13

call OptionsGetInt(PETSC_NULL_CHARACTER,'-m',m,flg,ierr) call OptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,ierr) !

call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) ！ 进程序号 call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr) ！ 进程个数 !

！创建系数矩阵A（m*n,m*n），按省缺AIJ系数格式存储，元素在运行时由PETSc ! 按连续行分块分配到各进程

call MatCreate(PETSC_COMM_WORLD,m*n,m*n,A,ierr) !

! 各进程获取自己包含的系数矩阵的行的范围

call MatGetOwnershipRange(A,Istart,Iend,ierr) !

! 5点格式有限叉分等网格步长离散二维规则区域上Dirichlet问题, ! 网格变量按自然序排列,系数矩阵元素赋值 do 10 II=Istart,Iend-1 v = -1.0

i = II/n ! 网格点所在的行坐标 j = II - i*n ! 网格点所在的列坐标 if ( i.gt.0 ) then JJ = II - n

call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr) !A(II,JJ)=A(II,JJ)+v endif

if ( i.lt.m-1 ) then JJ = II + n

call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr) !A(II,JJ)=A(II,JJ)+v endif

if ( j.gt.0 ) then JJ = II - 1

call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr) !A(II,JJ)=A(II,JJ)+v endif

if ( j.lt.n-1 ) then JJ = II + 1

call MatSetValues(A,1,II,1,JJ,v,ADD_VALUES,ierr) !A(II,JJ)=A(II,JJ)+v endif v = 4.0

call MatSetValues(A,1,II,1,II,v,ADD_VALUES,ierr) !A(II,JJ)=A(II,JJ)+v 10 continue !

! 并行系数矩阵形成

call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)

14

call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) !

! 创建长度为m*n的PETSc向量,并由PETSc省缺平均分配到各进程中, ! 也可以由用户在命令行参数中指定

call VecCreateMPI(PETSC_COMM_WORLD,PETSC_DECIDE,m*n,u,ierr) !

call VecSetFromOptions(u,ierr) ! 接受用户的特殊指定,若无,则为省缺 call VecDuplicate(u,b,ierr) ! 向量复制 call VecDuplicate(b,x,ierr) ! 向量复制 !

! 设置精确解,对应Drichlet边界条件

call OptionsHasName(PETSC_NULL_CHARACTER, & \"-random_exact_sol\

! 命令行输入是否采用随机数产生器生成精确解 if (flg .eq. 1) then

call PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,

& rctx,ierr) ! 创建PETSc随机数产生器

call VecSetRandom(rctx,u,ierr) ! 由随机数产生器创建精确解向量 call PetscRandomDestroy(rctx,ierr) ! 释放随机数产生器 else

call VecSet(one,u,ierr) ! 精确解赋值为1 endif

call MatMult(A,u,b,ierr) ! 矩阵向量乘,计算右端项 !

! 可视化观察精确解

call OptionsHasName(PETSC_NULL_CHARACTER, & & \"-view_exact_sol\ if (flg .eq. 1) then

call VecView(u,VIEWER_STDOUT_WORLD,ierr) endif !

! 创建线性解法器对象

call SLESCreate(PETSC_COMM_WORLD,sles,ierr) ! 设置预条件矩阵为自身

call SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN, & ierr) !

! 命令行输入确定迭代算法、预条件技术、收敛条件（残差下降比例）、 ！ 是否可视化观察残差迭代下降历史

-6

！ 省缺为GMRES（30）、块JACOBI（ILU（0））、10、否。

! -ksp_type -pc_type -ksp_monitor -ksp_rtol call SLESSetFromOptions(sles,ierr)

15

！

! 求解线性系统

call SLESSolve(sles,b,x,its,ierr) ！

！误差计算

call VecAXPY(neg_one,u,x,ierr) call VecNorm(x,NORM_2,norm,ierr) if (rank .eq. 0) then

if (norm .gt. 1.e-12) then write(6,100) norm, its else

write(6,110) its endif endif

100 format('Norm of error ',e10.4,' iterations ',i5) 110 format('Norm of error < 1.e-12, iterations ',i5) ！

! 释放内存空间

call SLESDestroy(sles,ierr) call VecDestroy(u,ierr) call VecDestroy(x,ierr) call VecDestroy(b,ierr) call MatDestroy(A,ierr) ！

! 退出PETSc系统

call PetscFinalize(ierr) end

16

2．例二、（petsc-2.0.28/src/snes/examples/tutorial/ex5f.F）

！ 采用PETSc非线性解法器部件SNES并行求解二维规则区域上 ！ 非线性固体燃料点火Bratu方程

! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1 , ! u = 0 for x = 0, x = 1, y = 0, y = 1.

！ 程序运行： mpirun -np ex5f [-help] [all PETSc options] ! 命令行参数：

! -par ： SFI参数lambda ：0 <= lambda <= 6.81； ! -mx ： X方向网格点个数； ! -my ：Y方向网格点个数； ! -Nx ：X方向进程个数； ! -Ny ：Y方向进程个数； !

program main implicit none !

#include \"ex5f.h\"

SNES snes Vec x, r Mat J

ISLocalToGlobalMapping isltog

integer its, Nx, Ny, matrix_free, flg, N, ierr, m double precision lambda_max, lambda_min

external FormFunction, FormInitialGuess, FormJacobian

call PetscInitialize(PETSC_NULL_CHARACTER,ierr) call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr) call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)

lambda_max = 6.81 lambda_min = 0.0 lambda = 6.0 mx = 4 my = 4

call OptionsGetInt(PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr) call OptionsGetInt(PETSC_NULL_CHARACTER,'-my',my,flg,ierr)

call OptionsGetDouble(PETSC_NULL_CHARACTER,'-par',lambda,flg,ierr) if (lambda .ge. lambda_max .or. lambda .le. lambda_min) then if (rank .eq. 0) write(6,*) 'Lambda is out of range' SETERRA(1,0,' ')

17

endif N = mx*my

call SNESCreate(PETSC_COMM_WORLD,SNES_NONLINEAR_EQUATIONS, & & snes,ierr)

Nx = PETSC_DECIDE Ny = PETSC_DECIDE

call OptionsGetInt(PETSC_NULL_CHARACTER,'-Nx',Nx,flg,ierr) call OptionsGetInt(PETSC_NULL_CHARACTER,'-Ny',Ny,flg,ierr)

if (Nx*Ny .ne. size .and. & & (Nx .ne. PETSC_DECIDE .or. Ny .ne. PETSC_DECIDE)) then if (rank .eq. 0) then

write(6,*) 'Incompatible number of procs: Nx * Ny != size' endif

SETERRA(1,0,' ') endif

call DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX,mx,& & my,Nx,Ny,1,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)

call DACreateGlobalVector(da,x,ierr) call DACreateLocalVector(da,localX,ierr) call VecDuplicate(x,r,ierr)

call VecDuplicate(localX,localF,ierr)

call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, & & PETSC_NULL_INTEGER,ierr)

call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, & & PETSC_NULL_INTEGER,ierr)

xs = xs+1 ys = ys+1 gxs = gxs+1 gys = gys+1

ye = ys+ym-1 xe = xs+xm-1 gye = gys+gym-1 gxe = gxs+gxm-1

call SNESSetFunction(snes,r,FormFunction,PETSC_NULL_OBJECT,ierr)

18

call OptionsHasName(PETSC_NULL_CHARACTER,'-snes_mf',matrix_free, & & ierr) if (matrix_free .eq. 0) then if (size .eq. 1) then

call MatCreateSeqAIJ(PETSC_COMM_WORLD,N,N,5, & & PETSC_NULL_INTEGER,J,ierr) else

call VecGetLocalSize(x,m,ierr)

call MatCreateMPIAIJ(PETSC_COMM_WORLD,m,m,N,N,5, & & PETSC_NULL_INTEGER,3,PETSC_NULL_INTEGER,J,ierr) endif

call SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL_OBJECT, & & ierr)

call DAGetISLocalToGlobalMapping(da,isltog,ierr) call MatSetLocalToGlobalMapping(J,isltog,ierr)

endif

call SNESSetFromOptions(snes,ierr)

call FormInitialGuess(x,ierr) call SNESSolve(snes,x,its,ierr) if (rank .eq. 0) then write(6,100) its endif

100 format('Number of Newton iterations = ',i5)

if (matrix_free .eq. 0) call MatDestroy(J,ierr) call VecDestroy(x,ierr) call VecDestroy(r,ierr)

call VecDestroy(localX,ierr) call VecDestroy(localF,ierr) call SNESDestroy(snes,ierr) call DADestroy(da,ierr) call PetscFinalize(ierr) end

subroutine FormInitialGuess(X,ierr) implicit none

19

#include \"ex5f.h\" Vec X integer ierr

Scalar lx_v(0:1) PetscOffset lx_i ！

ierr = 0

call VecGetArray(localX,lx_v,lx_i,ierr)

call ApplicationInitialGuess(lx_v(lx_i),ierr) call VecRestoreArray(localX,lx_v,lx_i,ierr)

call DALocalToGlobal(da,localX,INSERT_VALUES,X,ierr) ！

return end

subroutine ApplicationInitialGuess(x,ierr) implicit none #include \"ex5f.h\"

Scalar x(gxs:gxe,gys:gye) integer ierr

integer i, j, hxdhy, hydhx

Scalar temp1, temp, hx, hy, sc, one

ierr = 0 one = 1.0

hx = one/(dble(mx-1)) hy = one/(dble(my-1)) sc = hx*hy*lambda hxdhy = hx/hy hydhx = hy/hx

temp1 = lambda/(lambda + one)

do 20 j=ys,ye

temp = dble(min(j-1,my-j))*hy do 10 i=xs,xe

if (i .eq. 1 .or. j .eq. 1 & & .or. i .eq. mx .or. j .eq. my) then x(i,j) = 0.0 else

x(i,j) = temp1 * & & sqrt(min(dble(min(i-1,mx-i)*hx),dble(temp)))

20

endif 10 continue 20 continue

return end

subroutine FormFunction(snes,X,F,dummy,ierr) implicit none #include \"ex5f.h\" SNES snes Vec X, F integer dummy integer ierr

Scalar lx_v(0:1), lf_v(0:1) PetscOffset lx_i, lf_i

call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr) call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr) call VecGetArray(localX,lx_v,lx_i,ierr) call VecGetArray(localF,lf_v,lf_i,ierr)

call ApplicationFunction(lx_v(lx_i),lf_v(lf_i),ierr) call VecRestoreArray(localX,lx_v,lx_i,ierr) call VecRestoreArray(localF,lf_v,lf_i,ierr)

call DALocalToGlobal(da,localF,INSERT_VALUES,F,ierr) call PLogFlops(11*ym*xm,ierr)

return end

subroutine ApplicationFunction(x,f,ierr) implicit none #include \"ex5f.h\"

Scalar x(gxs:gxe,gys:gye), f(gxs:gxe,gys:gye) integer ierr

Scalar two, one, hx, hy, hxdhy, hydhx, sc Scalar u, uxx, uyy integer i, j

one = 1.0 two = 2.0

21

hx = one/dble(mx-1) hy = one/dble(my-1) sc = hx*hy*lambda hxdhy = hx/hy hydhx = hy/hx

do 20 j=ys,ye do 10 i=xs,xe

if (i .eq. 1 .or. j .eq. 1 & & .or. i .eq. mx .or. j .eq. my) then f(i,j) = x(i,j) else

u = x(i,j)

uxx = hydhx * (two*u & & - x(i-1,j) - x(i+1,j))

uyy = hxdhy * (two*u - x(i,j-1) - x(i,j+1)) f(i,j) = uxx + uyy - sc*exp(u) endif 10 continue 20 continue

return end !

subroutine FormJacobian(snes,X,jac,jac_prec,flag,dummy,ierr) implicit none #include \"ex5f.h\"

SNES snes Vec X

Mat jac, jac_prec MatStructure flag integer dummy integer ierr Scalar lx_v(0:1) PetscOffset lx_i

call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr) call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr) call VecGetArray(localX,lx_v,lx_i,ierr)

call ApplicationJacobian(lx_v(lx_i),jac,jac_prec,ierr) ！

22

call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr) call VecRestoreArray(localX,lx_v,lx_i,ierr)

call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr) flag = SAME_NONZERO_PATTERN

call MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,ierr)

return end !

subroutine ApplicationJacobian(x,jac,jac_prec,ierr) implicit none #include \"ex5f.h\"

Scalar x(gxs:gxe,gys:gye) Mat jac, jac_prec integer ierr

integer row, col(5), i, j

Scalar two, one, hx, hy, hxdhy, hydhx, sc, v(5)

one = 1.0 two = 2.0

hx = one/dble(mx-1) hy = one/dble(my-1) sc = hx*hy hxdhy = hx/hy hydhx = hy/hx

do 20 j=ys,ye

row = (j - gys)*gxm + xs - gxs - 1 do 10 i=xs,xe row = row + 1 ! boundary points

if (i .eq. 1 .or. j .eq. 1 & & .or. i .eq. mx .or. j .eq. my) then

call MatSetValuesLocal(jac,1,row,1,row,one, & & INSERT_VALUES,ierr) ! interior grid points else

v(1) = -hxdhy v(2) = -hydhx

v(3) = two*(hydhx + hxdhy) & & - sc*lambda*exp(x(i,j))

23

v(4) = -hydhx v(5) = -hxdhy

col(1) = row - gxm col(2) = row - 1 col(3) = row col(4) = row + 1 col(5) = row + gxm

call MatSetValuesLocal(jac,1,row,5,col,v, & & INSERT_VALUES,ierr) endif 10 continue 20 continue

return end

24

六、PETSc的具体应用与比较

1． 求解二维规则区域上时间依赖非线性对流扩散方程：

 均匀网格有限差分离散+ 时间步进 +非线性系数静止 线性化 + Krylov子空间迭代 + 预条件技术 ；  算法与程序实现技术：

 QPJM：Qmrcgstab+点Jacobi预条件+手工MPI实现；  GBJP：Gmres + BILU预条件 + PETSc实现（SLES）；  BPJP：Bicgstab + 点Jacobi预条件 +PETSc实现（SLES）；  BBJP：Bicgstab + BILU预条件 + PETSc实现（SLES）；  微机机群（9 P-II 400 MHz，100Mbps 交换机）上运行时间（秒，400*400网格规模）： 方法 \\ P QPJM GBJP BPJP BBJP 2 1 （每线性迭代） 6704 （0.284） 3473 5300 （0.779） 2827 1084 （0.283） 612 560 （0.497） 337 4 1842 1478 332 190 6 1191 1124 231 141 8 979 702 171 95  Qmrcgstab与Gmres的数值性能相当，但仅为BiCGSTAB的1/5-1/6；

 BILU比点Jacobi能获取一倍的性能提高，但每次迭代时间为点Jacobi的1.7倍，且迭代次数随处理机个数增加而增加，幅度小于20%；

25

 与自身串行计算时间比较，四个方法的并行加速比大致相同；

 比较QMJM和BPJP的每次线性迭代时间（0.283秒）： PETSc的每次线性迭代子程序的并行性能可与手工MPI相当；

 代码开发周期：手工MPI = 20天，PETSc = 5 天；  代码长度：手工MPI = 2400行，PETSc = 400 行；  迭代方法选择明显优于手工：PETSc非常灵活，只需在命令行给出选项，即可适用任何网格规模、进程个数与划分、迭代方法等，代码不需重编译。

2．二维三温多物质非定常流体力学问题：

 算法：二维柱坐标 + 结构网格离散 + Lagrange方法 + 显示格式求解流体力学方程 + 隐式格式非线性块GS迭代求解三温能量方程 + 时间依赖；

 PVM并行实现：时间3个月，Master-Slave模式，含动态负载平衡方法，无重要参数和网格移动的实时可视化输出，代码7100行；

 PETSc实现：时间20天，仅使用“分布阵列（DA）、向量、图形等”最简单的部件，沿用原程序的数值迭代方法和计算子程序，SPMD模式，无动态负载平衡，含重

26

要参数和网格移动的实时可视化输出，代码6300行；  性能比较：同样的算法和固定网格划分的静态负载平衡方法，PVM程序与PETSc程序在微机机群和Origin 2000上的并行性能比较

表2 微机机群上一个时刻的模拟时间比较

网格 方法\\P PVM 40*53 PETSc PVM 160*53

表3 DSM并行机上一个时刻的模拟时间比较 网格 方法\\P PVM 40*53 PETSc PVM 160*53

PETSc 1 508 395 (22%) 2194 1725 (21%) 2 269 204 (24%) 1097 863 (21%) 4 166 8 105 16 - - 225 PETSc 1 1456 1525 (-5%) 6171 60 (-6%) 2 794 812 (-2%) 3226 3413 (-6%) 4 509 4 (4%) 1903 1968 (-3%) 6 403 375 (7%) 1455 1343 (6%) 8 362 319 (12%) 1125 1048 (7%) 124 84 (25%) (20%) 658 345 524 286 180 (20%) (17%) (20%) 27

 微机机群上，PETSc由于引入了面向对象的额外开销，而又没有利用PETSc内部的计算功能部件，故串行性能低于MPI版本（5%-6%），但随处理机个数的增加，并行性能逐渐高于PVM版本（7%），两个方面的原因：1）微机机群上，MPI通信性能高于PVM约5%-10%；2）PETSc的MPI通信开销可与手工相当，性能是较好的。

 DSM并行机上，PETSc版本的单机性能比PVM版本高22%，这是因为PETSc内部的编译优化选项专门针对微处理器调整到最优，而用户一般很难做到。注意：这里的PVM版本的采用编译优化选项已经与PETSc采用的一致，否则，性能相差更大！

 DSM并行机上，PETSc版本的并行性能基本保持为PVM版本的20%-25%，这主要是由于串行性能高20%引起的，也说明DSM并行机上，PVM的通信性能与MPI大致相同，而PETSc的通信通信性能也可与手工相当！

28

七、PETSc的优点、缺陷与建议

1．优点

 一般库软件具备的高性能、可移植优点，PETSc均具备；  面向对象技术使得PETSc内部功能部件的使用非常方便，接口简单而又适用面广，大量缩短开发周期和减少工作量；

 PETSc完全兼容MPI系统，即在PETSc程序中，可任意调用MPI函数；

 PETSc通信性能可与手工MPI程序相比较；

 PETSc提供强大而又丰富的运行时选项数据库、内部功能部件性能的实时监控系统、并行性能分析系统；  PETSc具有很好的软件可继承、可集成、可发展性，与其他库软件可以互操作；

2．缺陷：

 尽管PETSc不需要用户懂得面向对象技术，但要求用户必须基于它所提供的数据结构来组织并行程序的开发；  PETSc仍然要求用户对MPI并行程序设计模式（甚至MPI函数）有较好的理解；  学习（入门）有一定的难度； 3．建议：

29

 用户首先熟悉MPI的使用，最好对面向对象程序设计方法有一定的了解后，再学习使用PETSc；

 对结构、非结构网格问题，尤其涉及到隐格式线性代数方程组求解，采用PETSc可以大大减少重复工作，达到事半功倍的效果，且能保证程序的可移植性、高效性，且一旦PETSc升级，应用程序也自动升级；

 对大型应用软件，尤其是库软件的的设计，也可以借鉴PETSc的思想解决软件的可继承性、可维护性和可发展性问题；

30