索引与排序
在许多应用中希望用自由度、单元、节点等的一种或多种“排序”(或编号)。在并行环境中做这些更为复杂,因为每一台处理机不能保持不同排序之间映射的完整列表。另外,PETSc线性代数程序中所使用的排序(常是连续范围)可能不对应于应用的“自然”排序。PETSc提供一些实用程序允许你清晰并有效地处理不同排序,即AO对象。
用户通过AOCreateBasic()创立AO对象后,可用如下命令之一 AOPetscToApplication(AO ao,int n,int *indices); AOApplicarionToPetsc(AO ao,int n,int *indices);
就可以通过在“PETSc”和“应用”两种排序之间映射来得到所需要的排序。
而对于局部编号方案到PETSc整体编号方案映射索引,PETSc提供ISLocalToGlobalMapping语句实现。这通过下面的程序完成: ISLocalToGlobalMappingCreate(MPI_Comm,N,*globalnum, * ctx); ISLocalToGlobalMappingApply(ctx, n, *in, *out); ISLocalToGlobalMappingApplyIS(ctx,IS isin,isout); ISLocalToGlobalMappingDestroy(ctx);
这里N表示局部索引的个数,globalnum包含每个局部编号的整体编号,ISLocalToGlobalMapping是所得的PETSc对象,它包含需要用ISLocalToGlobalMappingApply()或ISLocalToGlobalMappingApplyIS()应用映射所需的信息。
分布阵列DA
分布阵列(DA)组件的设计是PETSc的一大特色,它是专门用于逻辑正规矩形网格的网格管理工具。DA的主要设计目的是为了获得软件包对矩形网格区域分裂(必须按行列划分)问题的自适应能力,对相关方程的并行求解过程完成数据通信的自动处理。DA的设计实现得益于矩形网格的逻辑正规性,即一旦确定了区域的划分以及网格上的微分方程离散格式,则计算过程中所需要的全部通信信息是能够统一定位的,考虑到矩形网格使用的普遍性,这种高度集成的型板使用模式具有很高的应用价值。
有了DA对象,用户处理基于矩形网格区域分裂的并行计算问题变得相对简单。一开始,用户可以基于网格规模和处理机个数创建一个DA对象,这个对象将为各当前进程记录下整体和局部网格的丰富信息,其后即可以通过使用DA对象直观的创建与网格分布对应的全局向量和本地局部向量,而不必重新关注向量的存储分布与排序,也可以通过DA对象方便的获取向量信息和网格信息等。特别的,当局部计算需要用到非局部信息的时候,用户可以通过调用DA操作隐式实现数据通信的自动完成,而不必显式的定位通信的操作地址,有些情况下PETSc甚至可以主动管理通信,让用户彻底的摆脱对通信的安排。
无结构网格管理与IS
PETSc利用一个索引集(IS)的概念来简化定义在非结构化网格问题的映象点校正的向量分散与聚集。一个索引集是一个整数索引集合的推广,用于分散、聚集及在向量和矩阵上类似的运算。
对几乎所有的非结构网格来说,通过处理机操作调入进行网格的部分分配和存储器在性能方面有较大的冲突。在大多数PDE计算中,在数值计算之前处理机通过一种“预处理”的方法来完成网格划分和分配。但这并不意味着它需要用一种单独的、串行的程序来完成,在实际的程序中在建立并行网格数据结构之前就已经完成。PETSc为ParMETIS提供了接口,允许用并行的方式去划分。PETSc不能为动态再划分提供正确地支持,通过在两台处理机间移动数据来进行调入平衡。对需要网格精确加细的问题,PETSc需用重建数据结构来逼近。
6
2.2.2矩阵Mat
与向量相同,矩阵结构(Mat)也是微分方程数值计算的基本概念之一,但相比而言,矩阵在具体设计内容上比向量要复杂得多,这源于矩阵结构在实际应用中具有更丰富的应用策略。由于没有单一的矩阵格式能对所有问题保持最佳,PETSc提供了多种矩阵实现方式,支持多种稀疏存储方式(包括稠密存储)、多种分布存储处理方式、多种矩阵索引和装配方式等。
PETSc所采用面向对象技术所带来的功能函数的多态特质在矩阵这种数据结构上得到了非常明显的体现,即在统一的矩阵类型定义下,对各种不同格式矩阵的不同操作细节,实现一种外部接口一致、内部执行各异。矩阵对象只要在创建初期按照当前需求正确设置矩阵的内部状态,此后对矩阵的使用是非常一般化的,PETSc会根据相应的矩阵信息自动实现对底层执行路径的判断,而无须用户关心更多的调用细节。这种接口与执行分离的函数多态思想贯穿于PETSc的各个组件,对于简化软件使用模式具有非常重要的意义。一个典型的例子是矩阵乘法:
MatMatMult(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 这一操作接口无论对稀疏矩阵、稠密矩阵,还是对串行矩阵、并行矩阵,在形式上始终如一,只要相关矩阵规模与分布是协调的,相应的操作总能被正确执行,这与用户关心概念、不关心细节的应用期望完全相符。
在编程过程中,矩阵对象的使用与向量类似,也需要一个创建、赋值、使用、注销的过程,矩阵对象一经创建,它就具备了若干基础状态,如存储格式、并行划分方式等,这些基本状态与运行过程中随时被加入的其它特殊状态,是在多态函数中寻找唯一确定的执行路径的关键。由于矩阵状态类型对程序性能的影响依赖于具体问题的具体情况,因此要想获得最佳性能,用户必须主动选择最合理的矩阵状态类型,PETSc并不试图让默认的选择方案总是最优,而是致力于提供一个简便的选择接口。
PETSc矩阵提供的运算函数如下表所示: 函数名 MatAXPY(Mat Y,Scalar *a,Mat X,MatStructure); MatMult(Mat A, Vec x, Vec y); MatMultAdd(Mat A, Vec x, Vec y, Vec z); MatMultTrans(Mat A, Vec x, Vec y, Vec z); MatMultTransAdd(Mat A, Vec x, Vec y, Vec z); MatNorm(Mat A, NormType type, double *r); MatDiagonalScale(Mat A, Vec l, Vec r); MatScale(Mat A, Scalar a); MatConvert(Mat A, MatType type, Mat *B); MatCopy(Mat A, Mat B, MatStructure); MatGetDiagonal(Mat A, Vec x); MatTranspose(Mat A,MatReuse Mat *B); MatZeroEntries(Mat A); MatShift(Mat Y, Scalar a); 表2:PETSc 矩阵运算
PETSc也支持不通过显示存储整个矩阵而是通过向量运算来实现矩阵的各种操作和运
7
运算 Y = Y + a * X y = A * x z = y + A * x y = AT * x z = yATx r =||A||type A = diag(l) * A * diag(r) A = a * A B = A B = A x = diag(A) B = AT A = 0 Y = Y + a * I
算的方法,即“无矩阵”运算(Matrix-Free Matrices),用户在这个环境下可以自由开发矩阵运算程序:
MatCreateShell():创建一个虚拟的矩阵对象 UserMult():用户编写的矩阵向量乘积程序
MatShellSetOperation():将一个用户程序封装到一个虚拟的矩阵对象中。
8
PETSc方程求解器
PETSc工具箱的方程组求解器组件定义在各种数据结构组件之上,包括线性方程组求解器KSP,非线性方程组求解器SNES,时间依赖方程解法器TS,其中线性方程组求解器是实现关键的第一步,后两者都是以此为基础的。基于这三个求解器用户可以开发自己的应用程序。
线性求解器KSP
线性求解器KSP是PETSc的核心部分,它对软件包的所有线性方程组解法器,包括串行和并行,直接法和迭代法,提供了统一和高效的访问。对于线性方程组Ax = b ,只要给出定义该方程组的数据结构A与b,线性求解器KSP都可以通过一条简单的函数调用KSPSolve()来获得求解。
Krylov子空间方法与一个预条件子结合是大多数迭代求解线性方程组之现代数值代码的核心。PETSc的线性求解器模块就是由子空间法模块KSP与预条件子PC模块构成。在版本2.20之前PETSc提供SLES对象作为最高功能的集成体,线性求解器均调用SLES接口,后该对象被移除,直接由对象KSP完成该功能。特别地,KSP对象对直接解法也提供和迭代解法相同的调用序列。
一个典型的求解线性方程组的求解代码如下:
KSPCreate(comm, *ksp); 创建一个KSP求解器的对象。
KSPSetOperators(ksp, Amat, Pmat, flag) ; 设置线性方程组的矩阵。
KSPSetFromOptions(ksp) ; 设置运行时选项,主要是KSP和PC选项。 KSPSolve(ksp, b, x) ; 运行求解器求解。 KSPDestroy(ksp) ; 销毁对象。
对KSP软件包的一般使用来说,上面的过程就足够了。求解器缺省的解法器是重启动GMRES,单机情形预条件是ILU(0),多机情形预条件是块Jacobi方法(每个处理机一块,每个块用ILU(0)求解)。用户可通过KSPGetPC(ksp,*pc)来提取PC预条件子的环境上下文,再调用任意一个PC或KSP程序去修改相应的缺省选项。
PETSc的Krylov子空间方法KSP对象提供了丰富的迭代方法供用户选择,为设置将要使用的Krylov子空间方法,你需要调用命令:
KSPSetType(ksp,KSPTypemethod);
或在程序运行时使用-ksp_type ksptype进行修改。
Krylov子空间法对象缺省的收敛检验是基于剩余的l2-范数。命令: KSPSetTolerances(ksp, rtol, atol, dtol,maxits);
可以设置收敛条件。运行时选项-ksp_monitor或其它监视器则显示剩余的范数。
表3总结了目前KSP所提供的迭代方法,对应数据库选项名和缺省收敛监控器。
9
方法 KSPType 选项数 据库名 richardson chebychev cg bicg gmres bcgs cgs tfqmr tcqmr cr lsqr preonly 缺省收敛 监控器+ true true true true precond precond precond precond precond precond precond precond Richardson KSPRICHARDSON Chebychev KSPCHEBYCHEV Conjugate Gradient KSPCG BiConjugate Gradient KSPBICG Generalized Minimal Residual KSPGMRES BiCGSTAB KSPBCGS Conjugate Gradient Squared KSPCGS Transpose-Free Quasi-Minimal Residual (1) KSPTFQMR Transpose-Free Quasi-Minimal Residual (2) KSPTFCQMR Conjugate Residual KSPCR Least Squares Method KSPLSQR Shell for no KSP method KSPPREONLY + true – 表示真实剩余范数, precond – 表示预条件剩余范数 表3:KSP 方法. 缺省时所有方法使用左预条件
Krylov空间方法常与一个预条件子结合使用,表4中总结了PETSc支持的基本预条件方法,详细的使用可查阅相关资料。
方法 Jacobi Block Jacobi SOR (and SSOR) SOR with Eisenstat trick Incomlete Cholesky Incomplete LU Additive Schwarz Linear Solver Combination of Preconditioners LU Cholesky No preconditioning Shell for user-defined PC PCType PCJACOBI PCBJACOBI PCSOR PCEISENSTAT PCICC PCILU PCASM PCSLES PCCOMPOSITE PCLU PCCholesy PCNONE PCSHELL 选项数据库名 jacobi bjacobi sor eisenstat icc ilu asm sles composite lu cholesky none shell 表4:PETSc预条件子
对于奇异线性方程组,即系数矩阵具有零空间,例如离散具有Neumann边界条件的Laplace算子,PETSc也提供工具来帮助求解这些方程组。用户需用PETSc的Vecs的阵列形式的正交基来存储它,然后通过MatNullSpaceCreate()创建一个MatNullSpace对象,再通过语句PCNullSpaceAttach()然后将零空间传递给PC对象。
10