PETSc并行程序设计方法.ppt

1 PETSc并行程序设计方法中科院计算机网络信息中心超级计算中心 2004年11月 北京 2 报告基本目标 全面介绍PETSc体系结构 基本特色和功能 促进PETSc应用程序开发方面的合作与交流 寻求PETSc软件开发方面的广泛建议和帮助 3 报告内容提要 一 PETSc的基本概况四 PETSc的范例简介二 PETSc的基本特色五 PETSc的其它功能三 PETSc的实例分析六 PETSc的简单总结 4 一 PETSc的基本概况 PETSc简单概述 PETSc体系结构 PETSc基本对象 PETSc核心组件 PETSc下载安装 5 1 概述 PETSc Portable ExtensibleToolkitforScientificComputation 是美国能源部ODE2000支持开发20多个ACTS工具箱之一 由Argonne国家实验室SatishBalay等人开发的可移植可扩展科学计算工具箱 主要用于高性能求解偏微分方程组及相关问题 PETSc所有消息传递均采用MPI标准实现 PETSc用C语言开发 遵循面向对象设计的基本特征 用户可以基于PETSc对象灵活开发应用程序 目前 PETSc支持Fortran77 90 C和C 编写的串行和并行代码 6 PETSc是系列软件和库的集合 三个基本组件SLES SNES和TS本身基于BLAS LAPACK MPI等库实现 同时又为TAO ADIC ADIFOR Matlab ESI等工具提供数据接口或互操作功能 并具有极好的可扩展性能 PETSc为用户提供了丰富的Krylov子空间迭代方法和预条件子 并提供错误检测 性能统计和图形打印等功能 PETSc最新版本为petsc 2 2 1 2003年8月发布 PETSc网站http www mcs anl gov petsc 7 2 PETSc的体系结构 不同于其它微分 代数方程解法器 PETSc为用户提供了一个通用的高层应用程序开发平台 基于PETSc提供的大量对象和解法库 用户可以灵活地开发自己的应用程序 还可随意添加和完善某些功能 如为线性方程求解提供预条件子 为非线性问题的牛顿迭代求解提供雅可比矩阵 为许多数值应用软件和数学库提供接口等 图1表示了PETSc在实现层次上的抽象 8 应用程序 BLAS LAPACK MPI PDE解法器 TS 时间步进 SNES 非线性方程解法器 SLES 线性方程解法器 KSP Krylov子空间方法 PC 预条件子 Draw 矩阵 向量 索引集 核心组件 数值方法 基本对象 运行环境 性能分析 图1PETSc实现的抽象层次 9 应用程序 用户在PETSc环境下基于PETSc对象和算法库编写的串行或并行应用程序 尽管PETSc程序完全在MPI上实现 但PETSc程序具有固定的框架结构 即有初始化 空间释放和运行结束等环境运行语句 PDE解法器 用户基于PETSc三个基本算法库 TS SNES和SLES 构建的偏微方程求解器 但它却不是PETSc的基本组件之一 TS 时间步进积分器 用于求解依赖时间或时间演化ODE方程 或依赖时间的离散化PDE方程 对于非时间演化或稳态方程 PETSc提供了伪时间步进积分器 TS积分器最 10 后依赖线性解法器SLES和非线性解法器SNES实现 PETSc为PVODE提供了接口 另外 TS用法简单方便 SNES 非线性解法器 为大规模非线性问题提供高效的非精确或拟牛顿迭代解法 SNES依赖于线性解法器SLES 并采用线性搜索和信赖域方法实现 SNES依赖于雅可比矩阵求解 PETSc既支持用户提供有限差分程序 同时又为用户提供了依赖ADIC等自动微分软件生成的微分程序接口 SLES 线性解法器 它是PETSc的核心部分 PETSc几乎提供了各种高效求解线性方程的解法器 既有串行求解也有并行求解 既有直接法求解也有迭代法求解 对于大规模 11 线性方程组 PETSc提供了大量基于Krylov子空间方法和预条件子的各种成熟而有效的迭代方法 以及其它通用程序和用户程序的接口 KSP Krylov子空间方法 广泛涉及Richardson方法 共扼梯度法 CG和BiCG 广义最小残差法 GMRES 最小二乘QR分解 LSQR 等 PC 预条件子 包括雅可比矩阵 分块雅可比矩阵 SOR SSOR方法 不完全Cholesky分解 不完全LU分解 可加性Schwartz方法 多重网格预条件子等 DRAW 应用程序的性能分析和结果显示 12 3 PETSc的基本对象 向量 Vectors 矩阵 Matrices 索引集 IndexSets 13 3 1向量 Vectors PETSc的基本数据对象 PETSc向量对象主要用于存储线性方程组的解和右端向量 对于规则的正交网格 PETSc自动对向量进行划分 并通过分布式存储向量 即DA对象 来管理 通过DA对象 用户可以简单地实现向量的分发 聚集 局部和全局之间的相互映像 边界点的通信等基本操作 DA对象隐藏了进程之间的通信 用户只需提供全局的向量结构和数值 但对于无结构网格 用户则需通过索引集 即IS对象 来实现向量的分发 聚集 映像 边界点的通信等基本操作 14 3 1 1向量对象的创建 赋值 聚集和释放VecCreateSeq 创建一个串行的PETSc向量VecCreateMPI 创建一个并行的PETSc向量VecCreate 创建一个PETSc向量对象VecSetSize 设置向量维数VecSetFromOptions 通过运行参数设置向量类型VecSet 将一个数值赋给向量的每个元素VecSetValues 分别给向量的每个元素插入或累加数值 15 VecAssemblyBegin 启动一个向量的创建VecAssemblyEnd 完成一个向量的创建VecDuplicate VecDuplicateVecs 复制一个 组向量VecDestroy VecDestroyVecs 释放一个 组向量 16 3 1 2向量对象的基本运算操作VecGetOwnershipRang 返回向量在划分区域 进程 中的下界和上界VecGetArray 获得向量在局部区域内的访存指针VecGetArray Vecv PetscScalar array VecRestoreArray 关闭在局部区域内的对向量的访存VecGetLocalSize 返回向量在局部区域内的维数VecGetSize 返回全局向量的维数 17 此外 还涉及向量的算术运算 如点积 VecDot 范数 VecNorm 最小值 VecMin 绝对值 VecAbs 倒数 VecReciprocal 各种乘积 VecAXPY VecAYPX VecWAXPY VecAXPBY VecScale 等 18 例一 如何通过DA存取向量元素DA首先获得局部向量的控制权 然后将局部向量的访存权交给局部数组 进行相应的读写等操作之后 再逐步释放对局部数组和局部向量的访存权和控制权 PetscScalar f u Veclocal global DAGetLocalVecor DAda Veclocal DAVecGetArray DAda Veclocal void u 19 DAVecGetArray DAda Vecglobal void f f i j u i j DAVecRestoreArray DAda Veclocal void u DAVecRestoreArray DAda Vecglobal void f DARestoreLocalVecor DAda Veclocal 20 3 1 3向量对象的排序PETSc排序与应用排序之间的差异主要来源于两个原因 局部序号和全局序号的不同 不同区域分解方法之间的差异 AOCreatBasic 定义一个新的排序映射AOCreatBasic MPI Commcomm intn int apordering int petscordering AO ao AOPetscToApplication 从PETSc排序得到应用排序的映射AOApplicationToPetsc 从应用排序得到PETSc排序的映射AOPetscToApplication AO ao intn int indices 21 AODestroy 释放一个排序映射AOCreatBasicIS 定义一个新的索引集排序映射AOPetscToApplicationIS 从PETSc索引集排序得到应用索引集排序的映射AOApplicationToPetscIS 从应用索引集排序得到PETSc索引集排序的映射 22 3 1 4规则网格的DA实现分布式数组 DA 是PETSc的一个重要特色 DA自动管理数据在局部进程之间的划分 消息传递和读写 用户只需提供数组的全局数值 并通过全局和局部之间的逻辑排序来引用数组元素 在PETSc学习过程中 用户要特别注意向量 Vector 和数组 Array 两个概念 前者着重强调独立于数组数值和存储之外的逻辑结构或排序上的抽象 而后者着重强调数组元素的访问和数值改写 23 DACreat2d 创建一个2维的DA数组DACreat2d MPI Commcomm DAPeriodicTypewrap DAStencilTypest intM intN intm intn intdof ints int lx int ly DA da DACreatGlobalVector 创建一个全局的DA向量DACreatLocalVector 创建一个局部的DA向量DAGlobalToLocalBegin 启动一个DA向量的分发DAGlobalToLocalBegin DAda Vecglobal InsertModeiora Veclocal DAGlobalToLocalEnd 完成一个DA向量的分发 24 DALocalToGlobal 完成一个DA向量的聚集DALocalToGlobal DAda Vecl InsertModemodel Vecg DALocalToLocalBegin 启动一个DA向量的局部分发DALocalToLocalBegin DAda Veclocal1 InsertModeiora Veclocal2 其中InsertMode有ADD VALUES和INSERT VALUES两种类型DALocalToLocalEnd 完成一个DA向量的局部分发DAGetScatter 获得一个DA在全局与局部 局部与局部之间的分发权 极少使用 25 DAGetLocalVector 获得一个局部DA向量的访问权DARestoreLocalVector 释放一个局部DA向量的访问权DAVecGetArray 获得一个局部DA数组的访存权DAVecRestoreArray 释放一个局部DA数组的访存权DAGetCorners 获得数组左边界的初始位置DAGetGhostCorners 获得数组伪边界的初始位置DAGetGlobalIndices 获得全局格点数 含伪边界 26 DAGetISLocalToGlobalMapping 获得局部到整体索引集之间的映射DASetLocalToGlobalMapping 设置一个向量的局部到整体之间的映射DAGetAO 从DA环境上获得一个应用排序DAGetAO DAda AO ao 27 3 2矩阵 Matrices PETSc的基本数据对象 与向量对象不同 一个PETSc的矩阵对象首先是在局部 各个进程 数据填充完成之后再对其进行全局聚集 然后再由PETSc对象统一管理和实现矩阵的各种运算和操作 当问题的计算规模较大时 稀疏矩阵的填充是影响程序性能的关键因素 PETSc的矩阵运算和操作包括矩阵的创建 插值 聚集 各种算术运算和释放 PETSc的各种矩阵运算和操作使用起来非常方便 用户无需关心矩阵的具体存储实现 PETSc还为用户提供了的单纯依赖向量来实现矩阵各种运算 或无矩阵运算 的程序接口 28 3 2 1矩阵对象的创建 聚集和释放MatCreate 创建一个PETSc矩阵对象MatSetValues 给一个矩阵赋值MatSetValues MatA intm im intn intin PetscScalar values INSERT VALUES ADD VALUES MatSetOption 设置一个矩阵的行或列主存储次序MatAssemblyBegin 启动一个矩阵的聚集MatAssemblyEnd 完成一个矩阵的聚集 29 MatGetOwnershipRang 获得矩阵的局部行划分的上 下界MatCreateSeqAIJ 创建一个串行压缩稀疏行格式矩阵对象MatCreateSeqAIJ PETSC COMM SELF intm intn intnz int nnz Mat A MatCreateMPIAIJ 创建一个并行压缩稀疏行格式矩阵对象MatCreateMPIAIJ MPI Comm intm intn intM intN intd nz int d nnz into nz int o nnz Mat A 30 MatCreateSeqDense 创建一个串行稠密格式的矩阵对象MatCreateMPIDense 创建一个并行稠密格式的矩阵对象MatDestroy 释放一个矩阵对象 31 例二 并行稀疏矩阵对象的存储格式 P1 d nz 2 d nnz 2 2 2 o nz 2 o nnz 2 2 2 P2 d nz 3 d nnz 3 3 2 o nz 2 o nnz 2 1 1 P3 d nz 1 d nnz 1 1 o nz 4 o nnz 4 4 32 3 2 2矩阵对象的基本运算操作MatZeroRows 初始化一个矩阵行MatZeroEntries 初始化一个矩阵MatConvert 矩阵变换MatGetRow 获得矩阵的一行元素的值MatRestorerow 释放对矩阵行的访问权MatCopy 复制一个矩阵此外 还涉及矩阵算术运算 包括矩阵向量乘积 MatMult 矩阵范数 MatNorm 和矩阵转置 MatTranspose 等 33 3 2 3 无矩阵 运算 无矩阵 运算 Matrix FreeMatrices 指不通过显式存储整个矩阵而是通过向量运算来实现矩阵的各种操作和运算的方法 PETSc给用户提供了一个虚拟的矩阵对象 环境 用户在这个环境下可以自由开发矩阵运算程序 MatCreateShell 创建一个虚拟的矩阵对象MatCreateShell MPI Comm intm intn intM intN void ctx Mat mat 34 UserMult 用户编写的矩阵向量乘积程序MatShellSetOperation 将一个用户程序封装到一个虚拟的矩阵对象中MatShellSetOperation Matmat MatOperationMATOP MULT void void int UserMult Mat Vec Vec 35 3 2 4矩阵对象的划分对于许多无结构网格的PDE求解 格点在各进程中的分布对计算性能具有非常重要的影响 目前PETSc不支持动态数据划分 动态负载平衡等技术 而是采取反复创建和释放的办法来优化网格的划分 另外 PETSc提供了对并行图形划分软件ParMETIS的接口 MatCreateMPIAdj 创建一个含邻接信息的并行矩阵对象MatPartitioningCreate 创建一个并行矩阵划分 36 MatPartitioningSetAdjacency 设置矩阵划分的邻接信息MatPartitioningSetFromoptions 运行参数设置矩阵划分MatPartitioningApply 启动一个并行矩阵划分MatPartitioningDestory 释放一个并行矩阵划分ISPartitioningToNumbering 获得矩阵局部索引集和排序 37 3 3索引集 IndexSets PETSc的基本对象之一 一系列数据操作对象的集合 专门用来管理无结构网格向量的分发 聚集 局部和全局之间的映像 边界点的通信等基本操作 它对于PETSc在应用程序中最终实现通信抽象具有重要的意义 38 ISCreatGeneral 创建一个索引集排序ISCreatGeneral MPI Comm intn int indices IS is ISCreatStride 创建一个具有跨度的索引集排序ISCreatStride MPI Comm intn intfirst intstride IS is ISDestroy 释放一个索引集排序ISGetSize 获得一个索引集排序的数据规模的大小ISStrideGetInfo 获得一个含跨度的索引集排序信息ISStrideGetInfo IS is int first int stride 39 ISGetIndices 获得一个索引集排序的所有索引列表ISRestoreIndices 释放一个索引集列表的存储空间ISCreatBlock 创建一个块索引集排序ISCreatBlock MPI Comm intbs intn int indices IS is ISBlockGetIndices 获得一个块索引集排序的所有索引列表ISBlockGetSize 获得一个块索引集排序的数据规模的大小ISBlockGetBlockSize 获得一个块索引集排序的块规模大小ISBlock 没有描述 40 VecScatterCreat 创建一个向量与向量之间的分发VecScatterCreat Vecx ISix Vecy ISiy VecScatter ctx VecScatterBegin 启动一个向量与向量之间的分发VecScatterEnd 完成一个向量与向量之间的分发VecScatterDestroy 释放一个向量与向量之间的分发VecCreatGhost 创建一个含伪边界点的PETSc向量对象VecCreatGhost MPI Comm intn intN intnghost int ghosts Vec vv 41 VecCreatGhostWithArray 创建一个含伪边界点的PETSc向量对象和数组VecCreatGhostWithArray MPI Comm intn intN intnghost int ghosts PetscScalar array Vec vv VecGhostGetLocalForm 获得一个含伪边界点向量的局部访问权VecGhostRestoreLocalForm 释放一个含伪边界点向量的局部访问权VecGhostUpdateBegin 启动一个含伪边界点向量的更新VecGhostUpdateEnd 完成一个含伪边界点向量的更新 42 4 PETSc的核心组件 线性方程解法器 SLES 非线性解法器 SNES 时间步进解法器 TS 43 4 1线性方程解法器 SLES SLES构成了PETSc最核心的部分 它不仅是几乎所有PDE方程求解器的基本内核 而且也是实现PETSc的其它两个核心组件SNES和TS的必不可少的部分 SLES求解线性方程组Ax b 1 4 1 其中解算子A是n n维非奇异矩阵 x是n维右端向量 b为n维解向量 下面从线性方程求解环境的创建 Krylov子空间方法和预条件子 PC 的选择 收敛性判据 LU直接求解等方面详细介绍SLES 44 4 1 1基本用法SLESCreate 创建一个线性方程求解环境SLESSetOperators 设置求解算子 矩阵 SLESSetOperators SLESsles MatAmat MatPCmat Mat structureflag 其中 flag 取 SAME NONZERO PATERN DIFFERENT NONZERO PATERN 或 SAME PRECONDITIONER 之一 SLESSetFromoptions 通过运行参数设置SLES运行选项 45 SLESSolve 启动一个线性方程求解器SLESSolve SLESsles Vecb Vecx int its SLESSetup 启动一个线性方程求解环境SLESDestroy 释放一个线性方程求解环境SLESGetPC 获得PC对象 环境的访问权SLESGetKSP 获得KSP对象 环境的访问权 46 4 1 2Krylov子空间方法Krylov子空间方法是求解大规模稀疏线性方程组和大型矩阵特征值问题最有效的方法之一 其收敛性质依赖于矩阵特征值或奇异值的分布 预处理技术是改善矩阵特征值分布和这类方法收敛性质的最有效手段 一个好的预条件子使得具体选取何种Krylov子空间方法在算法效率上没有太大差别 但是 如何构造一个好的预条件子本身是一个重要难题 因此 PETSc为用户提供了构造预条件子的用户程序接口 Krylov子空间方法涉及Richardson方法 共扼梯度法 CG和BiCG 广义最小残差法 GMRES 最小二乘QR分解 LSQR 等 47 KSPSetType 设置Krylov子空间方法的求解类型KSPSetType KSPksp KSPTypemethod 其中第二个参数为Krylov子空间方法的类型 包括 Richardson方法KSPRICHARDSON切比雪夫方法KSPCHEBYCHEV共扼梯度法KSPCG广义最小残差法KSPGMRES共扼转置拟最小残差法KSPTCQMR双共扼梯度平方法KSPBCGS 48 经典Gram Schmidt方法KSPCGS无转置拟最小残差法KSPTFQMR共扼残差法KSPCR最小二乘QR分解方法KSPLSQR双共扼梯度法KSPBICG无迭代直接求解方法KSPPREONLY 49 KSPRichardsonSetScale 设置Richardson方法的阻尼系数KSPRichardsonSetScale KSPksp doubledamping factor KSPChebychevSetEigenvalues 设置切比雪夫方法的参数KSPChebychevSetEigenvalues KSPksp doubleemax doubleemin 默认值emax 0 01 emin 100 0 没有描述KSPGMRESSetRestart 设置广义最小残差法的迭代步数KSPGMRESSetRestart KSPksp intmax steps KSPGMRESSetOrthogonalization 设置为经典的Gram Schmidt方法 其计算性能好 但数值稳定性差 50 KSPCGSetType 设置共扼梯度方法中的算子对称类型KSPCGSetType KSPksp KSPCGTypeKSP CG SYMMETRIC KSPSetInitialGuessNonzero 设置一个非零的初始猜值KSPSetInitialGuessNonzero KSPksp PetscTruthflg KSPSetPreconditionerSide 设置为右预条件子KSPSetPreconditionerSide KSPksp PCSidePC RIGHT 51 KSPSetNormType 设置为真实残差的范数类型KSPSetNormType KSPksp KSP UNPRECONDITIONED NORM 其中大多数方法默认为预处理后的残差范数KSPSetTolerances 设置最大误差界和最大迭代步数KSPSetTolerances KSPksp doublertol doubleatol doubledtol intmaxits 其中从第二个参数起依次为相对误差 默认值为rtol 10 5 绝对误差 默认值为atol 10 50 相对误差比率 默认值为dtol 105 和最大迭代步数 默认值为maxits 105 52 KSPSetConvergenceTest 启动用户收敛性测试程序的封装KSPSetConvergenceTest KSPksp int test KSPksp intit doublernorm KSPConvergedReason reason void ctx void ctx 其中ctx为用户程序运行环境KSPDefaultMonitor 输出每迭代步的残差KSPSingularValueMonitor 输出每迭代步预条件子的最大奇异值KSPTrueMonitor 输出每个迭代步的残差的详细信息 53 KSPLGMonitorCreate 创建一个图形输出对象 环境KSPLGMonitorCreate char display char title intx inty Intw inth PetscDrawLG lg KSPSetMonitor 设置输出模式KSPSetMonitor KSPksp KSPLGMonitor lg 0 KSPLGMonitorDestroy 释放一个图形输出对象 环境 54 KSPSetComputeEigenvalues 设置预条件子特征值的计算环境KSPComputeEigenvalues 计算预条件子的特征值KSPComputeEigenvalues KSPksp intn double realpart double complexpart double neig KSPComputeEigenvaluesExplicitly 用直接方法计算预条件子的特征值KSPComputeEigenvaluesExplicitly KSPksp intn double realpart double complexpart 55 KSPGetSolution 获得解向量的值KSPGetRhs 获得右边向量的值KSPBuildSolution 获得每个迭代步解向量的逼近值KSPBuildResidual 获得每个迭代步残差的逼近值 56 4 1 3预条件子PCSetType 设置预条件子的类型PCSetType PCpc PCTypemethod 其中第二个参数为预条件子的可选类型 包括 JacobiPCJAOBI块JacobiPCBJAOBISOR SSORPCSORSOR Eisenstat方法 串行 PCEISENSTAT 57 不完全Cholesky分解PCICC不完全LU分解 串行 PCILU可加性SchwarzPCASM线性解法器PCSLES组合预条件子PCCOMPOSITELU分解 串行 PCLUCholesky分解PCCholesky无预条件子PCNONE用户提供预条件子PCSHELL 58 PCILUSetLevels 设置ILU预条件子的优化级别PCILCCSetLevels 设置ILCC预条件子的优化级别PCILUSetReuseOrdering 复用上次ILU分解中的排序对象PCILUSetUseDropTolerance 没有描述PCILUSetUseDropTolerance PCpc doubledt doubledtcol intdtcount PCILUDTSetReuseFill 复用上次ILU分解中的填充模式PCILUSetUSeInPlace 设置ILU分解的本地存取方式PCILUSetAllowDiagonalFill 没有描述 59 MatCreateMPIRowbs 创建BlockSolve95矩阵稀疏行结构MatCreateMPIRowbs MPI Comm intm intM intnz int nnz Mat A PCSORSetOmega 设置SOR方法的松弛因子PCSORSetIterations 设置SOR方法的迭代步数PCSORSetSymmetric 设置SOR方法的对称类型PCLUSetUseInplace 设置LU分解的本地存取方式 60 PCBJacobiGetSubSLES 为子块获取SLES环境信息PCBJacobiGetSubSLES PCpc int n local int first local SLES subsles PCJacobiSetTotalBlocks 设置全局划分块数目PCBJacobiSetLocalBlocks 设置局部划分块数目PCASMGetSubSLES 为子块获取SLES环境信息 同上 PCASMSetTotalSubdomains 设置全局划分 索引集 PCASMSetLocalSubdomains 设置局部划分 索引集 61 PCASMSetType 设置可加性Schwarz预条件子的类型PCASMSetOverlap 设置重叠区域PCASMSetOverlap PCpc intoverlap 其默认值为1 PCShellSetApply 设置一个由用户提供预条件子的环境PCShellSetApply PCpc int apply void ctx Vec Vec void ctx PCShellSetSetup 设置由用户提供预条件子的启动程序PCShellSetSetup PCpc int setup void ctx 62 PCCompositeAddPC 添加一个新的预条件子及其类型PCCompositeAddPC PCpc PCTypetype 默认为可加性的复合类型PCCompositeSetType 设置预条件子的复合类型PCCompositeSetType PCpc PCCompositeTypetype 其中第二个参数为PC COMPOSITE ADDITIVE或PC COMPOSITE MULTIPLICATIVEPCCompositeSetUseTrue 设置为真实残差的范数类型 63 PCCompositeGetPC 设置单个预条件子的类型参数PCCompositeGetPC PCpc intcount PC subtype PCSLESGetSLES 将解法器设置为预条件子PCSLESSetUserTrue 将解法器算子设置为预条件子 64 MGSetLevels 设置多重网格预条件子的优化级别MGSetType 设置多重网格预条件子的类型MGSetCycles 设置W cycle参数MGSetNumberSmoothUp 设置向前平滑参数MGSetNumberSmoothDown 设置向后平滑参数MGGetCoarseSolve 设置粗网格求解器MGGetSmoother 设置平滑算子MGGetSmootherUp 设置向前平滑算子MGGetSmootherDown 设置向后平滑算子MGSetInterpolate 设置插值算子 65 MGSetRestriction 设置约束条件MGSetResidual 设置多重网格预条件子的残差MGSetRhs 设置多重网格预条件子的右边向量MGSetX 设置多重网格预条件子的解向量MGSetR 设置多重网格预条件子的残差 66 4 1 4奇异方程求解MatNullSpaceCreate 创建一个零子空间对象MatNullSpaceCreate MPI Comm PetscTruthhasconstants intdim Vec basis MatNullSpace nsp PCNullSpaceAttach 设置一个零子空间对象PCNullSpaceAttach PCpc MatNullSpacensp 67 4 2非线性方程解法器 SNES SNES非线性解法器基于牛顿迭代法 线性搜索和信赖域方法 并依赖线性解法器SLES实现 雅可比矩阵的求解是SNES解法器的重要组成部分 SNES求解以下形式的非线性方程组组F x 0 1 4 2 其中解算子F Rn Rn 68 4 2 1基本用法SNESCreate 创建一个非线性方程求解环境SNESSetType 设置非线性求解器的类型SNESSetType SNESsnes SNESTypemethod 其中第二个参数为非线性求解器的可选类型 包括 线性搜索SNESLS信赖域方法SNESTRJacobi检测SNESTEST 69 SNESSetFromOptions 通过运行参数设置SNES运行选项SNESSolve 启动一个非线性方程求解器SNESDestroy 释放一个非线性方程求解器SNESSetFunction 设置非线性函数SNESSetFunction SNESsnes Vecf int FormFunction SNESsnes Vecx Vecf void ctx void ctx SNESSetJacobian 设置雅可比矩阵SNESSetJacobian SNESsnes MatA MatB int FormJacobian SNESsnes Vecx MatA MatB MatStructure flag void ctx void ctx 70 4 2 2非线性解法器SNESSetLineSearch 设置线性搜索方法的求解环境SNESSetLineSearch SNESsnes int ls SNESsnes Vec Vec Vec Vec double double double void lsctx SNESSetTolerances 设置最大迭代步数和误差界SNESSetTolerances SNESsnes doubleatol doublertol doublestol intits intfcts stol xk 1 xk SNESSetTrustRegionTolerances 同上 信赖域方法 SNESSetTrustRegionTolerances SNESsnes doubletrtol 71 SNESSetConvergenceTest 设置用户编写的收敛性测试对象SNESSetConvergenceTest SNESsnes int test SNESsnes doublexnorm doublegnorm doublef SNESConvergedReasonreason void cctx void cctx 其中参数xnorm gnorm和f依次为当前迭代误差 当前步误差或当前步梯度范数 最优化问题 和函数值 SNESSetMonitor 设置输出模式SNESSetMonitor SNESsnes int mon SNESsnes intits doublenorm void mctx void mctx int monitordestroy void 其中参数its和norm依次为迭代次数和函数值 72 SNESGetSolution 从SLES环境获得解向量对象SNESGetSolution SNESsnes Vec x SNESGetFunction 从SLES环境获得函数对象SNESGetFunction SNESsnes Vec r void ctx int func SNESsnes Vec Vec void 73 4 2 3 无矩阵 方法在其迭代求解线性系统的过程中 SNES完全支持用户提供的 无矩阵 预条件子 同前面类似 PETSc首先给用户提供了一个虚拟的矩阵对象和环境 然后用户基于这个环境可以自由开发自己的雅可比矩阵乘积计算程序 MatCreateSNESMF 在SNES中创建一个虚拟的矩阵对象MatCreateMF 创建一个虚拟的矩阵对象MatSNESMFSetFunctionError 设置有限差分的逼近误差 74 MatSNESMFDefaultSetUmin 设置默认的有限差分增量MatSNESMFWPSetComputeNormA 消除有限差分右乘向量计算中的冗余计算MatSNESMFWPSetComputeNormU 消除有限差分方法解向量计算中的冗余计算MatSNESMFSetHHistory 设置有限差分增量MatSNESMFReSetHHistory 更新有限差分增量MatSNESMFGetH 获得有限差分增量 75 4 2 4有限差分雅可比逼近雅可比矩阵的计算包括两个部分 稀疏结构的着色或排序 分别按重新得到的行或列计算雅可比矩阵的所有非零元素 有限差分逼近使得在一次差分逼近中 可以同时计算出雅可比矩阵同一列中的所有元素 借助于右乘初始输入矩阵 可以在一次差分逼近中同时计算出处于不同列并相互正交的更多非零元素 有限差分逼近简单表示为 76 例一 稀疏雅可比矩阵的着色与右乘初始矩阵 77 MatGetColoring 给稀疏雅可比矩阵着色 仅适用于串行MatFDColoringGreate 创建一个与着色相关的数据结构ISColoringDestroy 释放一个与着色相关的数据结构 对象MatFDColoringFromOptions 通过运行参数创建与着色相关的有限差分计算对象SNESSetJacobian 在SNES环境中启动疏雅可比的着色和有限差分计算 78 4 3时间步进解法器 TS TS时间步进积分器 用于求解依赖时间或时间演化的ODE方程 或依赖时间的离散化后的PDE方程 TS主要求解如下时间依赖问题ut F u t 1 4 3 其中u为有限维解向量 上式通常为运用有限差分或有限元方法离散后的常微分方线性程组 对于非时间演化或稳态方程 PETSc提供了伪时间步进积分器 TS积分器最终依赖线性解法器SLES和非线性解法器SNES实现 PETSc为PVODE求解器提供了接口 79 4 3 1基本用法TSCreate 创建一个TS求解环境TSCreate MPI Commcomm TSProblemTypeproblemtype TS ts 其中第二个参数为算子类型 包括A A t 或F u t TSSetType 设置TS求解器的类型TSSetType TSts TSTypetype 其中第二个参数为求解器算子类型 包括TS EULER TS BULER和TS PSEUDO TSSetInitialTimeStep 设置初始时间和步长TSSetInitialTimeStep TSts doubletime doubledt 80 TSSetTimeStep 设置 更改时间步长TSGetTimeStep 获得时间步长TSSetDuration 设置最大时间步数TSSetDuration TSts doublemaxsteps doublemaxtime TSSetUp 启动TS求解环境TSDestroy 释放TS求解环境TSView 启动TS屏幕输出TSView TSts PetscViewerviewer 81 4 3 2求解时间依赖问题TSSetSolution 为TS求解器设置一个初始猜值TSSetRHSMatrix 设置右端算子 线性问题 TSSetRHSMatrix TSts MatA MatB int f TS double Mat Mat MatStructure void void fP 其中参数f为时间依赖问题 A t 的矩阵对象A及其预条件子B的计算函数 否则为PETSC NULL 参数fP为用户程序运行环境 82 TSSetRHSFunction 设置右端算子 非线性问题 TSSetRHSFunction TSts int f TS double Vec Vec void void fP TSSetRHSJacobian 设置右端算子的雅可比矩阵 非线性问题 TSSetRHSJacobian TSts MatA MatB int fjac TS double Vec Mat Mat MatStructure void void fP 83 4 3 3求解时间稳态问题TSPseudoSetTimeStep 设置TS求解器的时间步长TSPseudoSetTimeStep TSts int dt TS double void void dtctx 参数dt为用户用于求解下次时间步长的函数 TSPseudoincrementDtFromInitialDt 设置时间步长 从初始时间步长出发计算 而默认方式基于上次时间步长 TSPseudoSetTimeStepIncrement 设置时间步长的尺度因子 默认值为1 1 84 4 3 4PVODE求解器与TS的接口方式TSPVodeSetType 启动PVODE求解器TSPVodeGetPC PVODE从TS环境中获得预条件子TSPVodeSetTolerance 设置误差界TSPVodeSetTolerance TSts doubleabs doublerel TSPVodeSetGramSchmidtType 没有描述TSPVodeSetGMRESRestart 没有描述 85 5 PETSc的下载安装 在LINUX UNIX环境下安装PETSc之前 系统首先需要做如下配置 1 MPI的一种实现 2 BLAS和LAPACK的一个拷贝 3 可选软件包 然后 可按照以下步骤安装PETSc 1 设置环境变量PETSC DIR和PETSC ARCH 分别为PETSc主目录的路径和机器体系结构 例如setenvPETSC DIR home username petsc 2 1 3setenvPETSC ARCHrs6000 86 2 编辑文件 PETSC DIR bmake PETSC ARCH pack ages 以指定MPI LAPACK BLAS X windows等的路径及可选软件包 3 编辑文件 PETSC DIR bmake PETSC ARCH vari ables以从其缺省改变C C 或Fortran编译器的名字 4 在PETSc主目录中使用makeBOPT gall make log建立PETSc的调试版本或使用makeBOPT Oall make log来建立PETSc库的优化版本 标志BOPT确定建立哪种库 其它可选项是对C 版本的BOPT g c O c 和对复数版本的BOPT g complex O complex 87 二 PETSc的基本特色 计算能力 可兼容性 可扩展性 抽象数据类型 输出能力 88 1 计算能力 PETSc为用户提供了丰富的算法和函数库资源 三个求解器 SLES SNES和TS 构成了PETSc的核心组件 PETSc不仅为中小规模线性方程组的求解提供了高效的直接方法 还为大规模 稀疏 线性方程组的迭代求解提供了多种Krylov子空间方法和多种预条件子 89 2 可兼容性 一方面 PETSc具有很强的兼容能力 可在不同体系结构和不同操作系统环境高效运行 另一方面 PETSc本身基于高性能的线性代数库 BLAS和LAPACK 和MPI消息传递环境实现 同时又充分吸收和融入了其它优秀软件的优点 如无结构网格区域分解和雅可比矩阵求解等方面的功能 90 3 可扩展性 PETSc的可扩展性功能主要包括三个方面 计算性能的并行可扩展性 功能的可扩展性和计算能力的可扩展性 无论是在计算时间还是在浮点性能方面 PETSc提供的范例程序都有良好的线性加速比性能 面向对象的程序设计风格使得PETSc具有良好的功能扩展能力 作为一个高级应用程序开发平台 PETSc特别适合于用来开发大型应用程序 91