(北京应用物理与计算数学研究所).doc

PETSc简介莫 则 尧（北京应用物理与计算数学研究所）一、 PETSc的起源与现状二、 PETSc的成功应用典范三、 PETSc的体系结构四、 PETSc的核心组件五、 PETSc程序示例六、 PETSc的具体应用与比较七、 PETSc的优点与缺陷一、PETSc的起源与现状1 全称：l 并行可扩展科学计算工具箱l Parallel Extensible Toolkits for Scientific Computingl /petsc2起源：l 美国能源部ODE2000支持的20多个ACTS（Advanved Computing Test & Simulation program, 美国能源部超级计算中心/）工具箱之一，其中包括：n 能提供算法的工具： Aztec ：分布式存储并行机求解大规模稀疏线性代数系统库； Hypre ：线性系统求解预条件库； Opt+ ：串行非线性优化问题数值库； PETSc ：并行可扩展科学计算工具箱，提供大量面向对象的并行代数数据结构、解法器和相关辅助部件，适合并行可扩展求解PDE方程（有限差分、有限元、有限体离散的隐式和显示格式）； PVODE ：并行常微分方程求解库； ScaLAPACK ： SuperLU ：代数系统直接求解库；n 算法开发辅助工具： Global Arrays ：以共享存储并行程序设计风格简化分布存储并行机上程序设计； Overture ：网格生成辅助工具； POET ：并行面向对象环境与工具箱； POOMA ：并行面向对象方法与应用，提供大量适合有限差分和粒子类模拟方法的数据并行程序设计（HPF）的C+类；n 运行调试与支持工具：CUMULVS，Globus，PAWS，SILOON，TAU，Tulip；n 软件开发工具： ATLAS & PHiPAC ：针对当代计算机体系结构的高性能特征，自动产生优化的数值软件，可与手工调试的BLAS库相比较（？）； Nexus , PADRE, PETE；3现状l 时间 ：1995年现在；l 目前版本：PETSc-2.0.28 + patch，源代码公开（不包含用户自己加入的核心计算子程序）；l 核心人员：数学与计算机部，Argonne国家重点实验室，Satish Balay, William Gropp, Lois C.McInnes, Barry Smith；l 参研人员：相关访问学者（几十人次，不同组件实现）；l 可移植性：CRAY T3D，T3E，Origin 2000, IBM SP, HP UX, ASCI Red, Blue Mountain, NOWs，LINUX，ALPHA等；l 目前，已下载上百套；二、PETSc的成功应用典范1PETSc-FUN3D：l 参考：W.K.Anderson etc., Achieving high sustained performance in an unstructured mesh CFD applications, SC99. ()l FUN3D：四面体三维无结构网格离散、求解可压或不可压Euler和Navire-Stokes方程、串行程序、百万量级的非结构网格点， NASA Langley 研究中心W.K.Anderson开发，适合飞行器、汽车和水下工具的设计优化；l 核心算法：非线性方程拟时间步逼近定常解、隐格式离散、Newton线性化、Krylov子空间迭代算法、加性Schwarz预条件（每个子区域近似精确求解），具有很好的数值可扩展性（即非线性迭代次数不随处理机个数的增加而显著增加）；l 移植周期：五个月（初步1996.101997.3），包括熟悉FUN3D与网格预处理器、学习ParMetis无结构网格剖分工具并集成到PETSc中、加入和测试PETSc的新功能、优化FUN3D面向向量机的代码段到面向cache的代码段、PETSc移植（非常少的时间，小于20天），并行I/O与后处理阶段还没完成；l 并行性能：n 代码行从14400减少为3300行（77%），包含I/O；n 优化后，串行程序发挥各微处理器峰值性能的78%-26%；（附页1）n ONERA M6 Wing, 2.8百万个网格单元（11百万个未知量），5123072个 ASCI Red 节点（双Pentium Pro 333MHz，每节点一个进程），保持95%以上的并行效率，发挥峰值性能的22.48%；（附页2）n 其他并行机：CRAY T3E、Origin 2000、IBM SP； l 奖励：SC99 Gordon Bell最佳应用奖；2石油：21世纪新一代油藏数值模拟框架； （USA Texas 大学油藏数值模拟中心）3空气动力学数值模拟中多模型多区域耦合流场问题：（USA自然科学交叉学科重点项目）；4天体物理中恒星热核爆炸问题数值模拟；（USA Chicago大学）三、PETSc的体系结构PETSc层次TS时间步进PDE解法器SNES(无约束优化、非线性解法器)SLES(线性方程解法器)PC(预条件)KSP(Krylov子空间方法)DRAWIndex SetsVectorsMatricesLAPACKMPIBLAS四、PETSc的核心组件1程序设计技术l 面向对象程序设计风格 + 标准C语言实现；l 标准C语言（C+）和FORTRAN语言接口；l 强调以面向对象的数据结构为中心设计数值库软件，并组织科学数值计算程序的开发；l PETSc应用：a) 根据应用需求，通过调用PETSc子程序建立数据结构（如向量、规则网格阵列、矩阵等）；b) 调用PETSc各功能部件（如索引、排序、基于规则网格的拟边界数据分布与收集、线性解法器、非线性解法器、常微分方程解法器、简单的图形输出等）的子程序来对数据对象进行处理，从而获取PETSc提供的科学计算功能；2核心组件（component ，class）：l 向量（Vectors）： 创建、复制和释放指定长度和自由度的串行或MPI并行向量（数据段被自动分配到不同的进程）； 向量操作：元素的赋值与索引、类似于BLAS的向量运算、向量的可视化显示等；l 索引与排序（Index Sets，Ordering）： 向量和矩阵元素的局部与全局、自然与多色序号的对应关系； 建立和释放任意两个集合的元素之间的对应和映射关系； 适合无结构网格在进程间的的任意网格剖分，以及通过数据映射操作，完成相应无规则网格拟边界数据交换的消息传递；l 分布阵列（DA：Distributed Array）： 建立在规则网格之上，一维、二维和三维（i=1,.,L, j=1,.,N, k=1,.,N），自动或指定阵列在进程间的区域划分，并沿拟边界设置宽度任意（离散格式需求）的影象（ghost）数组，存储邻近进程在相应位置包含的网格点的数值； 元素索引点 阵列元素可包含多个自由度，且可以任意索引和访问，属于向量的一种特殊情形，多有向量操作均适合于它； 数值计算中应用最为广泛的数据结构； 局部序和全局序可以索引和映射；l 矩阵（Matrices）： 指定维数和自由度大小的串行和MPI并行矩阵的生成、复制和释放； 矩阵元素的索引与访问； 稀疏矩阵压缩存储格式：AIJ稀疏行、BAIJ块稀疏行、Bdiag块对角； 稠密矩阵； 类似BLAS的基本矩阵操作，以及矩阵元素的标准或可视化输出； 矩阵的隐式形成与使用； 基于无结构网格划分工具（ParMetis）的并行矩阵的形成与使用；l 线性代数方程解法器（SLES）： 基于稀疏矩阵与向量数据结构； SLES的建立、访问、设置和释放； 目前实现的解法器：Krylov子空间方法（GMRES、CG、CGS、BiCGSTAB、TFQMR、Richardson、Chebychev），预条件（Additive Schwarz、Block Jacobi（ILU）、Jacobi、serial ILU、serial ICC、serial LU）； 收敛性测试与监控（实时图形显示迭代误差下降趋势）；l 非线性代数方程与无约束优化方程解法器（SNES）： 基于稀疏矩阵、向量和SLES数据结构； SNES的建立、访问、设置和释放； Newton线性化：line serach、trust region； 收敛性测试与监控（实时图形显示迭代误差下降趋势）；l PDE或ODE时间依赖方程解法器（TS）： 基于稀疏矩阵、向量、SLES和SNES数据结构； TS的建立、访问、设置和释放； 方法：Euler、Backward Euler、拟时间步逼近定常解等；l 对象的打印、图形和可视化输出：l 选项数据库支持： 对所有用户指定的算法和功能部件的性能监控，可在MPI程序运行时由命令行参数输入，非常方便；mpirun np 4 example -ksp_type bcgs ksp_xmonitor 并行性能自动统计、输出（ log_summary）； 用户自定义选项（如网格规模、进程个数、图形可视化输出等）；3与其他库软件的功能互用与接口：l BlockSolve95（并行ICC（0）、ILU（0）预条件）；l ESSL（IBM 快速稀疏矩阵LU分解）；l Matlab（数据的图形和数值后处理）；l ParMeTis（并行无结构网格图剖分）；l PVODE（并行常微分积分）；l SPAI（并行稀疏近似逆预条件）；l SAMRAI，Overture（并行网格管理软件包）；五、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方向网格点个数 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) 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 -pc_type -ksp_monitor -ksp_rtol call SLESSetFromOptions(sles,ierr)！ ! 求解线性系统 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) end2例二、（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, ) 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) 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#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) 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 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)！ 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