2023年12月20日发(作者：六年级的第一学期数学试卷)
PETSc简介
莫 则 尧
（北京应用物理与计算数学研究所）
一、 PETSc的起源与现状
二、 PETSc的成功应用典范
三、 PETSc的体系结构
四、 PETSc的核心组件
五、 PETSc程序示例
六、 PETSc的具体应用与比较
七、 PETSc的优点与缺陷
1
一、PETSc的起源与现状
1． 全称：
并行可扩展科学计算工具箱
Parallel Extensible Toolkits for Scientific Computing
/petsc
2．起源：
美国能源部ODE2000支持的20多个ACTS（Advanved
Computing Test & Simulation program, 美国能源部超级计算中心/）工具箱之一，其中包括：
能提供算法的工具：
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 s, Barry
Smith；
参研人员：相关访问学者（几十人次，不同组件实现）；
可移植性：CRAY T3D，T3E，Origin 2000, IBM SP, HP
UX, ASCI Red, Blue Mountain, NOWs，LINUX，ALPHA等；
目前，已下载上百套；
4
二、PETSc的成功应用典范
1．PETSc-FUN3D：
参考：on etc., Achieving high sustained
performance in an unstructured mesh CFD applications,
SC’99. ()
FUN3D：四面体三维无结构网格离散、求解可压或不可压Euler和Navire-Stokes方程、串行程序、百万量级的非结构网格点， NASA Langley 研究中心on开发，适合飞行器、汽车和水下工具的设计优化；
核心算法：非线性方程拟时间步逼近定常解、隐格式离散、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
!
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 ( .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 ( .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 ( .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 ( .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\",flg,ierr)
! 命令行输入是否采用随机数产生器生成精确解
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\",flg,ierr)
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)
!
! 命令行输入确定迭代算法、预条件技术、收敛条件（残差下降比例）、
！ 是否可视化观察残差迭代下降历史
！ 省缺为GMRES（30）、块JACOBI（ILU（0））、10-6、否。
! -ksp_type
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
! 命令行参数：
! -par ： SFI参数lambda ：0 <= lambda <= 6.81；
! -mx
! -my
! -Nx
! -Ny
!
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网格规模）：
2
方法 P
1
（每线性迭代）
QPJM
6704
（0.284）
3473
GBJP
5300
（0.779）
2827
BPJP
1084
（0.283）
612
BBJP
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
1
1456
2
794
812
(-2%)
3226
3413
(-6%)
4
509
489
(4%)
1903
1968
(-3%)
6
403
375
(7%)
1455
1343
(6%)
8
362
319
(12%)
1125
1048
(7%)
1525
PETSc
(-5%)
PVM 6171
6540
(-6%)
160*53
PETSc
表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
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
更多推荐
网格,并行,迭代,性能,数值,矩阵
发布评论