您好,欢迎来到爱go旅游网。
搜索
您的当前位置:首页PETSc 用户指南

PETSc 用户指南

来源:爱go旅游网


超级计算环境基础并行软件平台建设与应用 并行软件开发小组系列测试报告之一

PETSc 用户指南

中科院计算机网络信息中心超级计算中心

http:://www.sccas.cn/sce Email: walls@sccas.cn

程强 迟学斌 冯仰德 王建 赵永华

NCIC-SC-001, SCCAS 2004年8月,北京

目录

前言

1 PETSc简介

1.1 概况 ……….……………………………………..……………………….1 1.2 体系结构 ………………………………………………………………… 1 1.3 基本特色 ………………………………………………………………….4 1.4 安装PETSc ……………………………………………………………….5

2 PETSc的基本对象

2.1 向量 ……………………………………………….……………….………7 2.1.1 创建和聚集 ………………………. ……………………….…….………7 2.1.2 基本运算操作 ………………………. ……………………….…….……8 2.1.3 索引和排序 …………………………. …………………….…….………8 2.1.4 规则网格与DA …………….……………………………….…….………9

2.1.5 无结构网格与IS …………………………………………………….……10

2.2 矩阵 ……………………………………………….………………. . . .… 11 2.1.1 创建和聚集 ………………………. ……………………….……………11 2.1.2 基本运算操作 ………………………. …………………….…….……...11 2.1.3 无矩阵运算 ………………………….…………………….…….………12 2.1.4 矩阵的划分…………………………………. ……………………….……12

3 PETSc的基本功能

3.1 线性方程求解……………………………………………………………13

3.1.1 基本用法 ………………………. ………………………. . . .….……… 13 3.1.2 Krylov子空间方法 ……………….………………………….…….……13 3.1.3 预条件子 …………………………….……………………….…….……14 3.1.4 奇异方程求解 ………………………….………………….…….………16 3.2 非线性方程求解 ………………………………………………………… 16 3.2.1 基本用法 ………………………. ………………………….…….………17 3.2.2 非线性解法器 ………………………. ……………………….…….……17 3.2.3 无矩阵方法 …………………………. …………………….…….………17

1

3.2.4 有限差分雅可比逼近 ……… ……….………………….… ….……… 18 3.3 时间步进积分 …………………………………………………………… 18 3.3.1 基本用法 ………………………. ………………………. . . .….……… 19 3.3.2 求解时间依赖问题 …………………. ……………………… . …….… 19 3.3.3 求解时间稳态问题 …………………. ……………………… . …….… 19 3.3.4 其它求解器 ………………… ………. …………………….……. …… 19 3.4 PETSc的其他功能 ……………………………………………………… 20 3.4.1 性能分析 ………………………. ……………………….……. . .………20 3.4.2 图形输出 ……………………………. ……………………….…….……22 3.4.3 调试和错误检测 ……………………. …………………….…….………23 3.5 PETSc与其它软件 ………………………………………………… …… 23 3.5.1 DMMG ………………………. …………………………….…….………24 3.5.2 ADIC/ADIFOR ……………………. ……………………….……. . ...… 24 3.5.3 Matlab …………………………. …………………….……. . . . . . .…… 24 3.5.4 ESI ………………………. ……………………….…….…. . . . . . . . . ..… 24

4 PETSc编程

4.1 PETSc程序范例 ………………………………………………………… 25 4.2 PETSc程序结构 ………………………………………………………… 28

5 PETSc范例测试

5.1 线性方程求解 ……………………………………………………………29 5.1.1 范例简介 …………………………. ……………………………………29 5.2 非线性方程求解 ……………..………………………………………… 30 5.2.1 范例简介 ………………………………………………………………. 30 5.3 时间步进积分 ……………………………………………………………31 5.3.1 范例简介 …………………………………………………………….… 31

6 PETSc测试总结 参考文献

2

前言

通过计算手段进行重大科学发现,已普遍为人们所共识。从传统科学与工程领域如航空航天、地震预报、气候预测、大型水利建设和石油地质勘探,到大型基因组测试、新药设计和新材料合成等新兴科学研究领域,无处不需要大规模数值模拟和科学计算。在当今社会,科学计算已经逐渐成为影响和关系到一个国家经济发展、科技进步和国家安全的关键性环节。

2002年10月,中国科学院确立了“十五”信息化建设规划项目“超级计算环境建设与应用”,而“基础并行软件平台建设”又为其核心开发内容之一。它以基础并行软件环境开发为主要目标,通过广泛搜集目前流行的数值和并行应用软件,进行深层次研究和开发,改进和提高计算性能,并最终移植到特定计算环境如正在建设的国家网格节点上,主要供国内从事科学计算的广大科技人员共享。同时,我们也正逐步开发一系列具有重要应用价值的数值软件和面向网格的科学计算平台,并率先在国内开展自动微分算法的应用基础研究,以及逐渐开发一系列自动微分软件包和数值微分计算环境。

在此背景下,我们将陆续推出PETSc、TAO、FFTW、DOUG等系列并行软件在深腾6800上的测试报告,这包括应用和开发两个方面的目标。一方面,我们希望通过这些相对成熟的并行软件在一些典型应用上的成功来改变或导向人们传统的“从无到有”编程思路,我们希望科学家与工程师们更多的依赖这些成熟的、高性能的软件资源以便更高层次开发其应用程序,从而获得可比较性的性能提高。另一方面,任何软件的“成熟性”和可靠性都只能是相对的,不同的应用和不同的体系结构会有不同的代价和性能,以及体系结构的不断发展和变化与应用之间的相互作用对任何软件的生命期都可能是致命的。在典型应用的驱使下,我们希望能够基于这些并行软件做一些高层次的研究和开发。 基础并行软件平台建设与应用开发计划包括以下三个方面:

1. 数值并行软件的移植及其性能优化。主要包括可移植、可扩展科学计算工具箱(PETSc)、大型稀疏线性方程组并行迭代求解(AZTEC)、高级最优化工具箱(TAO)、非线性与微分/代数方程解法器(SUNDIALS)、FFTW、LAPACK、ScalLAPACK、稀疏矩阵特征值问题并行求解(PARPACK)和无结构网格上的区域分解(DOUG)等内容。

2. 自主开发一些高性能数值与并行计算软件包。这些软件包主要包括并行特征值求解器(PSEPS)、自动微分转换系统(DFT)、自动伴随生成器系统(ADG)和并行多维FFT软件包(PMDFFT)等内容。

3

3. 举办用户培训和相关学术会议,推动并行软件的应用和研究开发。主要包括PETSc、TAO、DOUG等基础并行数值软件的应用,以及LINUX/UNIX/网格应用、并行计算方法及其相关内容。 最后,我们热忱欢迎国内外所有来自不同学科、不同专业领域的科学家和工程技术人员就更为广泛的科学和工程计算问题提出建议。我们的联系方式: 单位:中国科学院计算机网络信息中心超级计算中心 网站:http://www.sccas.cn/sce 电话:(86-10) 5881-2132 邮件:walls@sccas.cn

地址:北京中关村南4街4号

通信:北京349信箱,SC,100080

4

1 PETSc简介

1.1 概况

PETSc (Portable, Extensible Toolkit for Scientific Computation) 是美国能源部ODE2000支持开发的20多个ACTS工具箱之一,由Argonne国家实验室开发的可移植可扩展科学计算工具箱,主要用于在分布式存储环境高效求解偏微分方程组及相关问题。PETSc所有消息传递通信均采用MPI标准实现。 PETSc用C语言开发,遵循面向对象设计的基本特征,用户基于PETSc对象可以灵活开发应用程序。目前,PETSc支持Fortran 77/90、C和C++编写的串行和并行代码。

PETSc是系列软件和库的集合,三个基本组件SLES、SNES和TS本身基于BLAS、LAPACK、MPI 等库实现,同时为TAO、ADIC/ADIFOR、Matlab、ESI 等工具提供数据接口或互操作功能,并具有极好的可扩展性能。PETSc为用户提供了丰富的Krylov子空间迭代方法和预条件子,并提供错误检测、性能统计和图形打印等功能。

如今,越来越多的应用程序在PETSc环境上开发,并逐渐显示出PETSc在高效求解大规模数值模拟问题方面的优势和威力。

PETSc最新版本为petsc-2.2.1,PETSc网站:http://www.mcs.anl.gov/petsc。 目前,PETSc 3正在开发中。

1.2 体系结构

不同于其它微分/代数方程解法器,PETSc为用户提供了一个通用的高层应用程序开发平台。基于PETSc提供的大量对象和解法库,用户可以灵活地开发自己的应用程序,还可随意添加和完善某些功能,如为线性方程求解提供预条件子、为非线性问题的牛顿迭代求解提供雅可比矩阵、为许多数值应用软件和数学库提供接口等。图1表示了PETSc在实现层次上的抽象,图2具体列举了PETSc的基本数值部件。这里做简要说明。

应用程序:用户在PETSc环境下基于PETSc对象和算法库编写的串行或并行应用程序。尽管PETSc完全在MPI上实现,但PETSc程序具有固定的框架结构,即有初始化、空间释放和运行结束等环境运行语句。

PDE解法器:用户基于PETSc的三个基本算法库(TS、SNES和SLES)构建的偏微方程求解器。但它却不是PETSc的基本组件。

1

TS:时间步进积分器,用于求解依赖时间或时间演化的ODE方程,或依赖时间的离散化后的PDE方程。对于非时间演化或稳态方程,PETSc提供了伪时间步进积分器。TS积分器最终依赖线性解法器SLES和非线性解法器SNES来实现。PETSc为PVODE库提供了接口。另外,TS的用法非常简单方便。 SNES:非线性解法器,为大规模的非线性问题提供高效的非精确或拟牛顿迭代解法。SNES依赖于线性解法器SLES,并采用线性搜索和信赖域方法实现。SNES非常依赖于雅可比矩阵求解,PETSc既支持用户提供的有限差分程序,同时又为用户提供了依赖ADIC等自动微分软件生成的微分程序接口。

SLES:线性解法器,它是PETSc的核心部分。PETSc几乎提供了各种高效求解线性方程的解法器,既有串行求解也有并行求解,既有直接法求解也有迭代法求解。对于大规模的线性方程,PETSc提供了大量基于Krylov子空间方法和预条件子的各种成熟而有效的迭代方法,以及其它通用程序和用户程序的接口。

应用程序 PDE解法器 TS (时间步进) SNES (非线性方程解法器) SLES (线性方程解法器) KSP (Krylov子空间方法) PC (预条件子) Draw 索引集 矩阵 向量 BLAS LAPACK MPI

图1:PETSc实现的层次结构

KSP:Krylov子空间方法,广泛涉及Richardson方法,共扼梯度法(CG和BiCG),广义最小残差法(GMRES),最小二乘QR分解(LSQR)等。 PC:预条件子,包括雅可比矩阵,分块雅可比矩阵,SOR/SSOR方法,不完全Cholesky分解,不完全LU分解,可加性Schwartz方法,多重网格预条件子等。 DRAW:应用程序的性能分析和结果显示。

2

非线性解法器

牛顿迭代法 线搜索 信赖域 其它 Euler 方法 时间步法 向后Euler方法 拟时间步 其它 Krylov子空间方法 GMRCG CGS Bi-CG-StaTFQMR Richardson Chebyshev 其它 预条件子 加法Schwarz 块Jacobi Jacobi ILU ICC LU 其它 矩阵 压缩稀疏行(AIJ) 块压缩稀疏(BAIJ) 块对角(BDiag) 稠密 其它

向量 索引 索引集 块索引 跨度 其它 图2:PETSc的数值组件

矩阵:PETSc的基本数据对象。与向量对象不同,一个PETSc的矩阵对象首先是在局部(各个进程)数据填充完成之后再对其进行全局聚集,然后再由PETSc对象统一管理和实现矩阵的各种运算和操作。当问题的计算规模较大时,稀疏矩阵的填充将是影响程序性能的关键因素。PETSc还专门为用户提供了的单纯依赖向量来实现矩阵基本运算(或无矩阵运算)的接口。

向量:PETSc的基本数据对象。对于规则的正交网格,PETSc自动对向量进行划分,并通过分布式存储向量(即DA对象)来管理。通过DA对象,用户可以简单地实现向量的分发、聚集、局部和全局之间的相互映像、边界点的通信等基本操作。DA对象隐藏了进程之间的通信,用户只需提供全局的向量结构和数值。但对于无结构网格,用户则可以通过索引集(即IS对象)来实现向量的分发、聚集、映像、边界点的通信等基本操作。

索引集:它是一系列数据操作对象的集合,专门用来管理无结构网格向量的分发、聚集、局部和全局之间的映像、边界点的通信等基本操作。

3

1.3 基本特色

众所周知,用户通过PETSc来开发应用程序往往具有相当的难度。一方面,它需要用户本身具有较高的数值计算方法方面的专业知识和并行计算方法方面的编程技巧。另一方面,总的来说PETSc只是一个高级的应用程序开发环境,它为许多软件(库)和用户程序提供接口,用户只有充分熟悉和利用现有的软件资源和数学库的基础上才有可能开发出高效的应用程序。尽管如此,PETSc仍然因为其具有其它软件不可比拟的优点吸引着越来越多的用户用它来开发应用程序。下面我们一一介绍和分析PETSc的这些优点或特色。

计算能力:PETSc为用户提供了丰富的算法和库资源。三个求解器(SLES、SNES和TS)构成了PETSc的核心组件。PETSc不仅为中小规模线性方程组的求解提供了高效的直接方法,还为大规模(稀疏)线性方程组的迭代求解提供了多种Krylov子空间方法和多种预条件子。

可兼容性:一方面,PETSc具有很强的兼容能力,可在不同体系结构和不同操作系统环境高效运行。另一方面,PETSc本身基于高性能的线性代数库(BLAS和LAPACK)和MPI消息传递环境实现,同时又充分吸收和融入了其它优秀软件的优点,如无结构网格区域分解和雅可比矩阵求解等方面的功能。 可扩展性:PETSc的可扩展性功能主要包括三个方面:计算性能的并行可扩展性、功能的可扩展性和计算能力的可扩展性。无论是在计算时间还是在浮点性能方面,PETSc提供的范例程序都有良好的线性加速比性能。面向对象的良好设计风格使得PETSc具有良好的功能扩展能力。作为一个高级应用程序开发平台,PETSc特别适合于用来开发大型应用程序。

抽象数据类型:PETSc基于面向对象技术实现,具有所有面向对象软件的可移植性、可继承性和可扩展性等基本程序特征。PETSc 的向量、矩阵等基本数据对象完全采用抽象数据类型实现,尽量对用户屏蔽数据对象的区域分解和存储等细节。所有PETSc格点数据对象的划分、初始化和存取等基本操作都由DA对象来管理和相应PETSc库函数实现。用户基于PETSc对象可以灵活开发其应用程序。PETSc对象和组件为构造大规模应用程序奠定一个良好的基础。 输出能力:PETSc具有良好的性能分析和图形输出功能。同时,PETSc还具有高可用性,并具有很强的错误诊断能力。

总之,PETSc在无论是在计算能力、设计风格还是在可兼容性和可扩展性等方面都显示出极大的优越性。PETSc不但为科学与工程计算领域的科学家和工程师提供了强大的(大规模)偏微方程求解工具,而且也为模型科学应用和高效算法设计提供了一个丰富的试验平台和计算环境。它使得算法的扩展和应用程序的个性化实现都更为容易。另外,PETSc的这种设计风格增强了代码的再利用性和编程的灵活性,同时将并行性问题与算法的选择分离开来。

4

1.4 安装PETSc

这里,我们以PETSc 2.1.3为例,介绍如何在LINUX/UNIX环境下安装PETSc。在安装PETSc之前, 系统首先需要做如下配置:

1)MPI的一种实现。对由厂家提供并已经安装了MPI的并行机,推荐使用厂家提供的MPI实现, 因为通常比使用免费版本具有更好的性能。否则推荐MPICH, 可在http://www.mcs.anl.gov/mpi/mpich免费下载获得。由于PETSc本身包含了MPI的一个简捷版本,因此只对串行运行PETSc感兴趣的用户除外。

2)BLAS和LAPACK的一个拷贝。许多机器本身提供了BLAS或LAPACK等数学库。例如,DEC alpha提供了DXML, 而IBM rs6000提供了ESSL版本。检查你使用的机器配置,如果这些库尚未在你的机器上安装,可在如下网站上获得 ftp://info. mcs.anl.gov/ pub/petsc/fblaslapack.tar.gz。建议你尽可能使用厂家提供的BLAS库。

3)可选软件包。PETSc提供了许多软件和库的接口,用户可以根据不同的应用需要安装相应的软件和库。 然后,你就可以按以下步骤逐步安装PETSc:

1)在http://www-unix.mcs.anl.gov/petsc/petsc-2/download/index.html上尽可能获取最新版本的petsc.tar.gz,并用如下方式解包 gunzip -c petsc.tar.gz | tar xof –

缺省时, 将自动创建petsc-2.1.3或其最新版本目录。 PETSc版本的补丁程序可通过http://www-unix.mcs.anl.gov/petsc/petsc-patches.html获得。 2)设置环境变量PETSC_DIR和PETSC_ARCH,分别为PETSc主目录的路径和机器体系结构。例如 setenv PETSC_DIR/home/username/petsc-2.1.3 setenv PETSC_ARCH rs6000

3)编辑文件${PETSC_DIR}/bmake/${PETSC_ARCH}/packages,以指定MPI, LAPACK,BLAS,X-windows等的路径及可选软件包。

注意: 如果只安装单机运行的PETSc版本, 则不必安装MPI,用户只需要在${PETSC_DIR}/bmake/${PETSC_ARCH}/packages中设置

MPI_LIB = ${PETSC_DIR}/lib/lib${BOPT}/${PETSC_ARCH}/libmpiuni.a MPI_INCLUDE = -I${PETSC_DIR}/src/sys/src/mpiuni MPIRUN = ${PETSC_DIR}/src/sys/src/mpiuni/mpirun

5

4)或许还需编辑文件${PETSC_DIR}/bmake/${PETSC_ARCH}/variables以从其缺省改变C, C++, 或Fortran编译器的名字: Solaris使用GNU编译器用PETSC_ARCH=solaris_gnu IBM rs6000使用GNU编译器用PETSC_ARCH=rs6000_gnu CRAY t3d: 确信环境变量TARGET被设置给cray-t3d

5) 在PETSc主目录中使用 make BOPT=g all >& make_log 来建立PETSc的调试版本或使用 make BOPT=O all >& make_log 来建立PETSc库的优化版本。标志BOPT确定建立哪种库。其它可选项是对C++版本的BOPT= [g_c++, O_c++]和对复数版本的BOPT=[g_complex,O_complex] 在HP-UX机器上我们强烈推荐使用Gnu make (如果还没用,可安装之),而不用PETSc提供的makefile脚本。在PETSc构建文件(makefiles)中使用gnumake, 需在文件${PETSC_DIR}/bmake/hpux/base定义OMAKE为gnumake的路径。

6)检查make_log,以确定在安装过程中是否发生任何错误。查阅一般问题的帮助信息http://www-unix.mcs.anl.gov/petsc/petsc-patches.html。 另外,用户还可以通过configure安装PETSc。请参考PETSc网站,这里从略。

6

2 PETSc的基本对象

本章从对象的创建、划分、功能和索引等几个方面简单介绍PETSc的两个基本对象:向量和矩阵。所涉及的内容和顺序基本与PETSc用户手册一致。

2.1 向量

向量是最简单的PETSc对象。PETSc向量对象主要用于存储线性方程组的解和右端向量。PETSc提供AO对象来管理向量在全局和局部之间的索引、排序和映射。PETSc还提供了两个对象DA和IS,来分别管理向量在规则正交网格和无结构网格上各进程之间的分发、聚集和边界点的数据通信等操作。

2.1.1 创建和聚集

VecCreatSeq:创建一个串行的PETSc向量 VecCreatMPI:创建一个并行的PETSc向量 VecCreat:创建一个PETSc向量名 VecSetSizes:设置向量维数

VecSetFromOptions:通过运行参数设置向量数据类型 VecSet:将一个数值赋给向量的每个元素

VecSetValues:分别给向量的每个元素插入或累加数值 VecAssemblyBegin:启动一个向量的创建 VecAssemblyEnd:完成一个向量的创建 VecView:屏幕打印向量的值 VecDuplicate:复制一个向量 VecDuplicateVecs:复制一组向量 VecDestroy:释放一个向量 VecDestroyVecs:释放一组向量

VecCreatSeqWithArray:创建一个串行的PETSc向量(用户程序) VecCreatMPIWithArray:创建一个并行的PETSc向量(用户程序)

7

2.1.2 基本运算操作

VecGetOwnershipRang:返回向量在局部区域(进程)中的下界和上界 VecGetArray:返回向量在局部区域内的访存指针 VecRestoreArray:关闭在局部区域内的对向量的访存 VecGetLocalSize:返回向量在局部区域内的维数 VecGetSize:返回全局向量的维数

关于向量的算术运算,如点积(VecDot)、范数(VecNorm)、最小值(VecMin)等,请参考PETSc用户手册(Pages 38)。

2.1.3 排序

AOCreatBasic:定义一个新的排序映射

AOPetscToApplication:从PETSc排序得到应用排序的映射 AOApplicationToPetsc:从应用排序得到PETSc排序的映射 AOCreatBasicIS:定义一个新的索引集排序映射

AOPetscToApplicationIS:从PETSc索引集排序得到应用索引集排序的映射 AOApplicationToPetscIS:从应用索引集排序得到PETSc索引集排序的映射 AOView:查询一个排序映射 AODestroy:释放一个排序映射

ISLocalToGlobalMappingCreat:建立一个局部排序与全局排序之间的映射 ISLocalToGlobalMappingApply:获得局部排序到全局排序的映射 ISLocalToGlobalMappingApplyIS:获得局部索引集排序到全局排序的映射 ISLocalToGlobalMappingDestroy:释放一个局部排序与全局排序的映射 ISGlobalToLocalMappingApply:获得全局排序到一系列局部排序的映射 VecSetLocalToGlobalMapping:设置从一系列局部排序与全局排序的映射 VecSetValuesLocal:得到一系列局部排序

在2.1.5,我们还将详细介绍主要用于无结构网格划分情形下的索引集排序。

8

2.1.4 规则网格的DA实现

分布式数组(DA)是PETSc的一个重要特色。DA自动管理数据在局部进程之间的划分、消息传递和读写。用户只需提供数组的全局数值,并通过全局和局部之间的逻辑排序来引用数组元素。在PETSc学习过程中,用户要特别注意向量(Vector)和数组(Array)两个概念,前者着重强调独立于数组数值和存储之外的逻辑结构或排序上的抽象,而后者着重强调数组元素的访问和数值改写。

DACreat1d:创建一个1维的DA数组 DACreat2d:创建一个2维的DA数组 DACreat3d:创建一个3维的DA数组

DACreatGlobalVector:创建一个全局的DA向量 DACreatLocalVector:创建一个局部的DA向量 DAGlobalToLocalBegin:启动一个DA向量的分发 DAGlobalToLocalEnd:完成一个DA向量的分发 DALocalToGlobal:完成一个DA向量的聚集 DALocalToLocalBegin:启动一个DA向量的局部分发 DALocalToLocalEnd:完成一个DA向量的局部分发 DAGetScatter:启动一个DA向量的分发

DAGetLocalVector:获得一个局部DA向量的访问权 DARestoreLocalVector:释放一个局部DA向量的访问权 DAVecGetArray:获得一个局部DA数组的访存权 DAVecRestoreArray:释放一个局部DA数组的访存权 DAGetCorners:获得数组左边界的初始位置 DAGetGhostCorners:获得数组伪边界的初始位置 DAGetGlobalIndices:获得全局格点数(含伪边界)

DAGetISLocalToGlobalMapping:获得局部到整体索引集之间的映射 DASetLocalToGlobalMapping:设置一个向量的局部到整体的映射 MatSetLocalToGlobalMapping:设置一个矩阵的局部到整体的映射 DAGetAO:从DA环境上获得一个应用排序

9

2.1.5 无结构网格的IS实现

索引集(IS)是PETSc的基本对象之一。它管理和实现无结构网格上向量之间的各种数据分发与聚集、伪边界点数据交换等基本操作。它对于PETSc在应用程序中最终实现通信抽象具有重要的意义。

ISCreatGeneral:创建一个索引集排序

ISCreatStride:创建一个具有跨度的索引集排序 ISDestroy:释放一个索引集排序

ISGetSize:获得一个索引集排序的数据规模的大小 ISStrideGetInfo:获得一个索引集排序的跨度值 ISGetIndices:获得一个索引集排序的所有索引列表 ISRestoreIndices:释放一个索引集列表的存储空间 ISCreatBlock:创建一个块索引集排序

ISBlockGetIndices:获得一个块索引集排序的所有索引列表 ISBlockGetSize:获得一个块索引集排序的数据规模的大小 ISBlockGetBlockSize:没有说明(见在线手册) ISBlock:没有说明(见在线手册)

VecScatterCreat:创建一个向量与向量之间的分发 VecScatterBegin:启动一个向量与向量之间的分发 VecScatterEnd:完成一个向量与向量之间的分发 VecScatterDestroy:释放一个向量与向量之间的分发 VecCreatGhost:创建一个含伪边界点的PETSc向量

VecCreatGhostWithArray:创建一个含伪边界点的PETSc向量(用户程序) VecGhostGetLocalForm:获得一个含伪边界点向量的局部访问权 VecGhostUpdateBegin:启动一个向量伪边界点的更新 VecGhostUpdateEnd:完成一个向量伪边界点的更新

VecGhostRestoreLocalForm:释放一个含伪边界点向量的局部访问权

10

2.2 矩阵

矩阵是PETSc的基本对象。PETSc同时提供了稠密矩阵和稀疏行矩阵的基本运算功能,以及一些特殊格式(如“无矩阵”实现,无结构网格划分等内容)和用户提供的某些功能扩展和实现。PETSc的矩阵运算和操作主要包括矩阵的创建、插值、聚集、各种算术运算和释放。PETSc的各种矩阵运算和操作使用起来非常方便,用户无需关心矩阵的具体存储实现。

2.2.1 创建和聚集

MatCreate:创建一个PETSc矩阵对象 MatSetValues:给一个矩阵赋值

MatSetOption:设置一个矩阵的存储格式 MatAssemblyBegin:启动一个矩阵的聚集 MatAssemblyEnd:完成一个矩阵的聚集

MatGetOwnershipRang:获得一个矩阵的局部行划分的上界和下界 MatCreateSeqAIJ:创建一个串行压缩稀疏行格式的矩阵对象 MatCreateMPIAIJ:创建一个并行压缩稀疏行格式的矩阵对象 MatCreateSeqDense:创建一个串行稠密格式的矩阵对象 MatCreateMPIDense:创建一个并行稠密格式的矩阵对象

2.2.2 基本运算操作

MatZeroRows:初始化一个矩阵行 MatZeroEntries:初始化一个矩阵 MatConvert:矩阵变换

MatGetRow:获得矩阵的一行元素的值 MatRestorerow:释放对矩阵行的访问权 MatCopy:复制一个矩阵 MatView:屏幕输出一个矩阵

11

关于矩阵的算术运算,如矩阵向量乘积(MatMult)、矩阵范数(MatNorm)、矩阵转置(MatTranspose)、获得矩阵对角元素(MatGetDiagonal)等,请参考PETSc用户手册(Pages 57)。

2.2.3 “无矩阵”运算

所谓“无矩阵”运算(Matrix-Free Matrices),是指不通过显示存储整个矩阵而是通过向量运算来实现矩阵的各种操作和运算的方法。PETSc给用户提供了一个虚拟的矩阵对象/环境,用户在这个环境下可以自由开发矩阵运算程序。

MatCreateShell:创建一个虚拟的矩阵对象 UserMult:用户编写的矩阵向量乘积程序

MatShellSetOperation:将一个用户程序封装到一个虚拟的矩阵对象中

2.2.4 矩阵的划分

对于许多无结构网格的PDE求解,格点在各进程中的分布对计算性能具有非常重要的影响。目前PETSc不支持动态数据划分,动态负载平衡等技术,而是采取反复创建和释放的办法来优化网格的划分。另外,PETSc提供了对并行图形划分软件ParMETIS的接口。

MatCreateMPIAdj:创建一个含邻接信息的并行矩阵对象 MatPartitioningCreate:创建一个并行矩阵划分 MatPartitioningSetAdjacency:设置矩阵划分的邻接信息 MatPartitioningSetFromoptions:通过运行参数设置矩阵划分 MatPartitioningApply:启动一个并行矩阵划分 MatPartitioningDestory:释放一个并行矩阵划分 MatDestroy:释放一个含邻接信息的并行矩阵对象 ISPartitioningToNumbering:获得矩阵局部的索引集和排序

12

3 PETSc的基本功能

本章主要介绍PETSc的三个核心组件、性能分析和接口功能。PETSc的三个核心组件包括线性方程求解器(SLES)、非线性方程求解器(SNES)和时间步进积分器(TS)。所涉及的内容和顺序基本与PETSc用户手册(版本2.2.0)一致。

3.1 线性方程求解

SLES构成了PETSc最核心的部分。它不仅是几乎所有PDE方程求解器的基本内核,而且也是实现PETSc的其它两个核心组件SNES和TS的必不可少的部分。SLES求解线性方程组

Axb (3.1)

其中解算子A是nn维非奇异矩阵,b是n维右端向量,x为n维解向量。本节从线性方程求解环境的创建、Krylov子空间方法和预条件子(PC)的选择、收敛性判据、LU直接求解等方面详细介绍SLES。

3.1.1 基本用法

SLESCreate:创建一个线性方程求解环境 SLESSetOperators:设置求解算子(矩阵)

SLESSetFromoptions:通过运行参数设置SLES运行选项 SLESSolve:启动一个线性方程求解器 SLESDestroy:释放一个线性方程求解环境 SLESSetup:启动一个线性方程求解器 SLESGetPC:获得PC对象/环境的访问权 SLESGetKSP:获得KSP对象/环境的访问权

3.1.2 Krylov子空间方法

KSPSetType:设置Krylov子空间方法的类型 KSPRichardsonSetScale:没有描述

13

KSPChebychevSetEigenvalues:没有描述 KSPGMRESSetRestart:没有描述

KSPGMRESSetOrthogonalization:没有描述 KSPCGSetType:设置共扼梯度方法的对称类型 KSPSetInitialGuessNonzero:设置一个非零的初始猜值 KSPSetPreconditionerSide:设置预条件子的位置 KSPSetNormType:设置范数类型 KSPSetTolerances:设置最大迭代步数

KSPSetConvergenceTest:启动用户收敛性测试程序的封装 KSPDefaultMonitor:输出每迭代步的残差

KSPSingularValueMonitor:输出每迭代步预条件子的最大奇异值 KSPTrueMonitor:输出每个迭代步详细的残差信息 KSPLGMonitorCreate:创建一个图形输出对象/环境 KSPSetMonitor:设置输出模式

KSPLGMonitorDestroy:释放一个图形输出对象/环境 KSPSetComputeEigenvalues:设置预条件子特征值的计算环境 KSPComputeEigenvalues:计算预条件子的特征值

KSPComputeEigenvaluesExplicitly:用直接方法计算预条件子的特征值 KSPGetSolution:获得解向量的值 KSPGetRhs:获得右边向量的值

KSPBuildSolution:获得每个迭代步解向量的逼近值 KSPBuildResidual:获得每个迭代步残差的逼近值

3.1.3 预条件子

PCSetType:设置预条件子的类型

PCILUSetLevels:设置ILU预条件子的优化级别

14

PCILCCSetLevels:设置ILCC预条件子的优化级别 PCILUSetReuseOrdering:设置复用ILU分解中的排序对象 PCILUSetUseDropTolerance:没有描述 PCILUDTSetReuseFill:没有描述 PCILUSetUSeInPlace:没有描述 PCILUSetAllowDiagonalFill:没有描述 MatCreateMPIRowbs:没有描述

PCSORSetOmega:设置SOR方法的松弛因子 PCSORSetIterations:设置SOR方法的迭代步数 PCSORSetSymmetric:设置SOR方法的对称类型 PCLUSetUseInplace:没有描述 PCBJacobiGetSubSLES:没有描述 PCASMGetSubSLES:没有描述 PCJacobiSetTotalBlocks:没有描述 PCASMSetTotalSubdomains:没有描述 PCASMSetType:没有描述

PCBJacobiSetLocalBlocks:没有描述 PCASMSetLocalSubdomains:没有描述 PCASMSetOverlap:没有描述

PCShellSetApply:启动一个用户提供预条件子的环境 PCShellSetSetup:设置一个用户提供预条件子的环境 PCCompositeAddPC:累加一个新的预条件子 PCCompositeSetType:设置预条件子的复合类型 PCCompositeSetUserTrue:没有描述 PCCompositeGetPC:没有描述 PCILUSetFill:没有描述

PCSLESGetSLES:将解法器算子设置为预条件子

15

PCSLESSetUserTrue:没有描述

MGSetLevels:设置多重网格预条件子的优化级别 MGSetCycles:没有描述

MGSetNumberSmoothUp:没有描述 MGSetNumberSmoothDown:没有描述 MGGetCoarseSolve:设置粗网格算子 MGGetSmoother:设置平滑算子 MGGetSmootherUp:没有描述 MGGetSmootherDown:没有描述 MGSetInterpolate:没有描述 MGSetRestriction:没有描述

MGSetResidual:设置多重网格预条件子的残差 MGSetRhs:设置多重网格预条件子的右边向量 MGSetX:设置多重网格预条件子的解向量 MGSetR:设置多重网格预条件子的残差

3.1.4 奇异方程求解

MatNullSpaceCreate:没有描述 PCNullSpaceAttach:没有描述

3.2 非线性方程求解

SNES非线性解法器基于牛顿迭代法(线性搜索和信赖域方法),依赖线性解法器SLES实现。雅可比矩阵的求解是SNES解法器的重要组成部分。SNES求解以下形式的非线性方程组

F(x)0 (3.2)

其中解算子F:RnRn。

16

3.2.1 基本用法

SNESCreate:创建一个非线性方程求解环境 SNESSetType:设置非线性求解器的类型

SNESSetFromOptions:通过运行参数设置SNES运行选项 SNESSolve:启动一个非线性方程求解器 SNESDestroy:释放一个非线性方程求解器 SNESSetFunction:设置非线性函数 SNESSetJacobian:设置雅可比矩阵

3.2.2 非线性解法器

SNESSetLineSearch:设置线性搜索方法的求解环境 SNESSetTolerances:设置最大迭代步数和误差界 SNESSetTrustRegionTolerances:同上

SNESSetConvergenceTest:设置用户编写的收敛测试 SNESSetMonitor:设置输出模式 SNESGetSolution:获得解向量的值 SNESGetFunction:获得函数值

3.2.3 “无矩阵”方法

在其迭代求解线性系统的过程中,SNES完全支持用户提供的“无矩阵”预条件子。同前面类似,PETSc首先给用户提供了一个虚拟的矩阵对象和环境,然后用户基于这个环境可以自由开发自己的预条件子计算程序。

MatCreateSNESMF:在SNES环境中创建一个虚拟的矩阵对象 MatCreateMF:没有描述

MatSNESMFSetFunctionError:设置有限差分的逼近误差

17

MatSNESMFDefaultSetUmin:设置默认的有限差分增量 MatSNESMFWPSetComputeNormA:设置有限差分增量 MatSNESMFWPSetComputeNormU:设置有限差分增量 MatSNESMFSetHHistory:没有描述 MatSNESMFReSetHHistory:没有描述 MatSNESMFGetH:没有描述 MatSNESMFKSPMonitor:没有描述 MatSNESMF:没有描述

3.2.4 有限差分雅可比逼近

雅可比矩阵的计算包括两个部分:稀疏结构的着色或排序;分别按重新得到的行或列计算雅可比矩阵的所有非零元素。

MatGetColoring:给稀疏雅可比矩阵着色

MatFDColoringGreate:创建一个与着色相关的数据结构/对象 ISColoringDestroy:释放一个与着色相关的数据结构/对象 MatFDColoringFromOptions:没有描述

SNESSetJacobian:在SNES环境中启动疏雅可比的着色和有限差分计算 DAGetColoring:没有描述

3.3 时间步进积分

TS时间步进积分器,用于求解依赖时间或时间演化的ODE方程,或依赖时间的离散化后的PDE方程。TS主要求解如下时间依赖问题

utF(u,t) (3.1)

其中u为有限维解向量,上式通常为运用有限差分或有限元方法离散后的常微分方线性程组。对于非时间演化或稳态方程,PETSc提供了伪时间步进积分器。TS积分器最终依赖线性解法器SLES和非线性解法器SNES来实现。PETSc为PVODE求解器提供了接口。

18

3.3.1 基本用法

TSCreate:创建一个TS求解环境 TSSetType:设置TS求解器的类型 TSSetInitialTimeStep:设置初始时间和步长 TSSetTimeStep:设置时间步长 TSGetTimeStep:获得时间步长 TSSetDuration:设置最大时间步数 TSSetUp:启动TS求解环境 TSDestroy:释放TS求解环境 TSView:启动TS屏幕输出

3.3.2 求解时间依赖问题

TSSetSolution:设置初始条件 TSSetRHSMatrix:设置右端算子 TSSetRHSFunction:设置右端算子

TSSetRHSJacobian:设置右端算子的雅可比矩阵

3.3.3 求解时间稳态问题

TSPseudoSetTimeStep:设置初始时间和步长 TSPseudoincrementDtFromInitialDt:设置时间步长 TSPseudoSetTimeStepIncrement:设置时间步长

3.3.4 其它求解器

TSPVodeSetType:启动PVODE求解器

TSPVodeGetPC:PVODE从TS环境中获得预条件子 TSPVodeSetTolerance:设置最大迭代步数和误差界 TSPVodeSetGramSchmidtType:没有描述 TSPVodeSetGMRESRestart:没有描述

19

3.4 PETSc的其它功能

本节简单介绍PETSc的性能分析、图形打印、错误检测方面的功能。PETSc不仅能为应用程序提供完整的计算时间、存储代价和浮点性能方面的信息,还能为计算的每个阶段和每个程序对象提供详细的信息。另外,PETSc在软件接口方面的功能我们将在下一节详细介绍。

3.4.1 性能分析

PETSc的性能分析(Profiling)功能主要包括计算时间、浮点性能、存储开销和通信性能几个方面的评价。PETSc提供了以下三个运行选项:

-log_summary:屏幕输出完整的性能信息 -log_info:屏幕输出使用对象的详细信息 -log_trace:屏幕输出过程执行的跟踪信息

注意使用-log_info选项对程序执行性能的影响可能较大,而-log_ summary 选项则对程序执行性能的测试影响较小。另外,-log_info和-log_trace选项为用户调试程序提供了方便。图3.1是SLES中的范例ex2在曙光2000机群上使用-log_summary选项后屏幕输出的结果。为简单起见,这里删去了一些辅助信息。

Norm of error 0.00305772 iterations 137

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- PETSc Performance Summary: ex2 on a rs6000 named warrior33. D2000-2 with 8 processors, by walls Tue Apr 27 09:27:32 2004 Using Petsc Version 2.1.3, Patch 0, Released May 31, 2002

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- Max Max/Min Avg Total Time (sec): 4.580e+00 1.01035 4.545e+00

Objects: 7.000e+01 1.00000 7.000e+01

Flops: 2.313e+07 1.00157 2.312e+07 1.850e+08

Flops/sec: 5.102e+06 1.01194 5.087e+06 4.070e+07

MPI Messages: 2.860e+02 2.00000 2.502e+02 2.002e+03 MPI Message Lengths: 2.918e+05 2.00000 1.020e+03 2.043e+06 MPI Reductions: 3.650e+01 1.00000

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- Summary of Stages: --------- Time ---------- --------- Flops --------- ----- Messages ----- --Message Lengths -- ----- Reductions –---

Avg %Total Avg %Total counts %Total Avg %Total counts %Total

0: Main Stage: 4.5447e+00 100.0% 1.8496e+08 100.0% 2.002e+03 100.0% 1.020e+03 100.0% 2.920e+02 100.0%

20

---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------Event Count Time (sec) Flops/sec ------- Global -------- ------- Stage ---------- Total

Max Ratio Max Ratio Max Ratio Mess Avg len Reduct %T%F %M %L %R %T %F %M %L %RMflop/s

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- MatMult 142 1.0 4.3690e-01 1.7 1.01e+07 1.7 2.0e+03 1.0e+03 0.0e+00 7 11 99 100 0 7 11 99 100 0 48

MatSolve 142 1.0 3.2932e-01 1.6 1.25e+07 1.6 0.0e+00 0.0e+00 0.0e+00 5 11 0 0 0 5 11 0 0 0 62

MatLUFactorNum 1 1.0 9.8193e-03 1.5 2.74e+06 1.5 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 14

MatILUFactorSym 1 1.0 1.5052e-02 1.8 0.00e+00 0.0 0.0e+00 0.0e+00 1.0e+00 0 0 0 0 0 0 0 0 0 0 0

MatAssemblyBegin 1 1.0 4.3116e-02 4.0 0.00e+00 0.0 0.0e+00 0.0e+00 2.0e+00 1 0 0 0 1 1 0 0 0 1 0

MatAssemblyEnd 1 1.0 1.9479e+00 1.0 0.00e+00 0.0 1.4e+01 5.1e+02 7.0e+00 4 3 0 1 0 2 43 0 1 0 20

MatGetOrdering 1 1.0 1.5409e-02 6.2 0.00e+00 0.0 0.0e+00 0.0e+00 2.0e+00 0 0 0 0 1 0 0 0 0 1 0

VecMDot 137 1.0 8.8890e-01 1.6 1.52e+07 1.6 0.0e+00 0.0e+00 1.4e+02 17 36 0 0 47 17 36 0 0 47 74

VecNorm 143 1.0 4.4378e-01 2.3 3.04e+06 2.3 0.0e+00 0.0e+00 1.4e+02 8 3 0 0 49 8 3 0 0 49 11

VecScale 142 1.0 2.1669e-02 1.8 2.39e+07 1.8 0.0e+00 0.0e+00 0.0e+00 0 1 0 0 0 0 1 0 0 0 107

VecCopy 4 1.0 1.8249e-03 1.6 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 0 0 0 0 0 0 0 0 0 0 0

VecSet 149 1.0 4.4439e-02 1.6 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 1 0 0 0 0 1 0 0 0 0 0

VecAXPY 10 1.0 3.2451e-03 1.4 1.80e+07 1.4 0.0e+00 0.0e+00 0.0e+00

VecMAXPY 142 1.0 3.8137e-01 1.8 4.11e+07 1.8 0.0e+00 0.0e+00 0.0e+00 6 38 0 0 0 6 38 0 0 0 185

VecScatterBegin 142 1.0 6.7344e-02 2.3 0.00e+00 0.0 2.0e+03 1.0e+03 0.0e+00 1 0 99 100 0 1 0 99 100 0 0

VecScatterEnd 142 1.0 7.6762e-02 3.9 0.00e+00 0.0 0.0e+00 0.0e+00 0.0e+00 1 0 0 0 0 1 0 0 0 0 0

PCSetUp 2 1.0 1.5821e-01 2.3 2.57e+05 2.3 0.0e+00 0.0e+00 3.0e+00 3 0 0 0 1 3 0 0 0 1 1

PCSetUpOnBlocks 1 1.0 6.5269e-02 2.4 6.41e+05 2.4 0.0e+00 0.0e+00 3.0e+00 1 0 0 0 1 1 0 0 0 1 2 PCApply 142 1.0 3.9893e-01 1.6 1.00e+07 1.6 0.0e+00 0.0e+00 0.0e+00 7 11 0 0 0 7 11 0 0 0 51

KSPGMRESOrthog 137 1.0 1.1071e+00 1.2 1.82e+07 1.2 0.0e+00 0.0e+00 1.4e+02 22 71 0 0 47 22 71 0 0 47 119

SLESSetup 2 1.0 1.9429e-01 2.1 1.90e+05 2.1 0.0e+00 0.0e+00 3.0e+00 3 0 0 0 1 3 0 0 0 1 1

SLESSolve 1 1.0 2.1358e+00 1.1 1.14e+07 1.1 2.0e+03 1.0e+03 2.8e+02 46 100 99 99 96 46 100 99 99 96 86

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- Memory usage is given in bytes:

Object Type Creations Destructions Memory Descendants' Mem. Matrix 4 4 376772 227284 Index Set 5 5 25828 0 Map 13 13 1664 0 Vec 41 41 639316 0

Vec Scatter 1 1 3504 0

Krylov Solver 2 2 16872 571064 Preconditioner 2 2 180 340024 SLES 2 2 0 928140

===============================================================================================================

图3.1 一个PETSc程序使用-log_summary选项后的性能输出信息

21

图3.1给出了一个SLES应用程序的时间代价和性能分析,主要包括三个部分:整体时间代价和性能;每个计算阶段或对象的时间代价和性能;以及每个PETSc对象的存储开销。关于PETSc的其它性能分析选项和接口功能,请参考PETSc使用手册,由于这些不是我们讨论的重点,这里从略。

3.4.2 图形输出

PetscViewerDrawOpenX:创建一个绘图对象/环境 PetscViewerDrawGetDraw:获得绘图对象的控制权 PetscDrawOpenX:打开一个绘图窗口

PetscDrawSetCoordinates:改变绘图窗口的坐标位置 PetscDrawSetViewPort:改变绘图窗口的大小 PetscDrawLine:绘一条直线 PetscDrawFlush:没有描述

PetscDrawSetDoubleBuffer:没有描述 PetscDrawSynchronlizedFlush:清除缓冲区 PetscDrawPause:设置图形的显示时间 PetscDrawString:字符显示

PetscDrawStringVertical:纵向字符显示 PetscDrawStringSetSize:设置字符大小 PetscDrawStringGetSize:获得字符大小

PetscDrawLGCreate:创建一个曲线绘图对象/环境 PetscDrawLGAddPoint:描绘一个点 PetscDrawLGAddPoints:描绘一组点 PetscDrawLGDraw:启动一个曲线绘图对象 PetscDrawLGDestroy:释放一个曲线绘图对象 PetscDrawLGReset:重置一个曲线绘图对象/环境 PetscDrawLGSetLimits:设置图形显示的范围

22

PetscDrawLGGetAxis:获得坐标刻度值 PetscDrawAxisSetColors:设置坐标颜色 PetscDrawAxisSetLabels:设置坐标刻度 KSPLGMonitor:图形收敛性显示

KSPLGMonitorCreate:创建图形收敛性显示对象

3.4.3 调试和错误检测

PetscError:检查栈内所有错误信息输出对象

PetscPushErrorHandler:向栈添加一个错误信息输出对象 PetscPopErrorHandler:忽略上一个错误 PetscTraceBackErrorHandler:没有描述

PetscAbortErrorHandler:当遇到出错时结束程序运行 PetscAttachErrorHandler:当遇到出错时启动调试程序 PetscPushSingnalHandler:捕捉错误信息 PetscSetFPTrap:设置浮点错误检测 SETERRQ:预先设置错误信息 PetscMalloc:存储检测

CHKERRQ:错误检测和出错信息输出

3.5 PETSc与其它软件

PETSc可扩展性的另一个方面表现在其为非常广泛的一类数值软件和数学库提供了很方便的软件接口。主要包括以下几种类型:1.线性代数求解器,如AMG、BlockSolve95、DSCPACK、hypre、ILUTP、LUSOL、SPAI、SPOOLES、SuperLU、SuperLU_Dist;2.最优化软件,如TAO、Veltisto;3.离散化和网格生成和优化工具包,如Overture、SAMRAI、SUMAA3d;4.常微分方程求解器,如PVODE;5. 其它,如Matlab、ParMETIS。这里,我们将重点介绍PETSc手册(版本2.1.6)中描述的几个软件和库,即DMMG、ADIC/ADIFOR、Matlab和ESI。

23

3.5.1 DMMG

PETSc为DMMG提供了高端访问的接口,并最终通过DA对象来实现在各处理器间的数据划分。目前DMMG仅提供了分段线性插值方法和线性多重网格预条件子,但在实际应用程序开发中很容易得到扩展。

关于PETSc如何引用DMMG的用法请参见PETSc用户手册(版本2.1.6,103)。

3.5.2 ADIC/ADIFOR

ADIC/ADIFOR是自动微分工具中的典型代表。PETSc通过DMMG高端访问接口为引用ADIC/ADIFOR生成的求解雅可比矩阵(通常是稀疏的或“无矩阵”结构)的程序提供了方便。

3.5.3 Matlab

PETSc为用户使用Matlab提供了三种方式:PETSc输出一系列文件给Matlab;从正在运行的PETSc程序中传递数据给Matlab;在PETSc和Matlab之间交互传递数据。其中后两种方式需要Matlab提供的实时指令或脚本文件。

关于PETSc如何引用Matlab的详细用法参见PETSc用户手册(版本2.1.6,107)。

3.5.4 ESI

方程求解器接口环境(ESI)是一个通用的可扩展的线性解法器接口。PETSc为ESI提供了两种接口方式:在PETSc应用程序中使用ESI对象和在ESI应用程序中使用PETSc对象。

24

4 PETSc编程

PETSc与MPI消息传递并行环境完全兼容。在这个意义上,用户可以在PETSc上开发任何基于消息传递的应用程序。但是,PETSc更希望用户将其视为一个在分布式计算环境中的PDE数值模拟和科学计算的平台,用户基于PETSc提供的大量线性方程和非线性方程求解器、丰富的数值迭代方法和各种预条件子,大量的对象和库资源,以及软件接口来开发和调试应用程序。

总之,PETSc给用户提供了丰富的算法资源和可编程环境。对于SNES和TS求解器,用户通常还需提供计算函数及其雅可比矩阵的子程序。下面我们结合一个并行SLES程序范例来介绍PETS的编程思路和程序结构。

4.1 PETSc程序范例

图4.1是一个典型的PETSc程序范例,它使用有限差分来求解二维Laplace问题。用SLES(Krylov子空间方法和预条件子)求解一个线性方程组。其代码见${PETSC_DIR}/src/sles/examples/tutorial/ex2.c。

/*$Id: ex2.c,v 1.94 2001/08/07 21:30:54 bsmith Exp $*/

/* 运行方式: mpirun -np ex2 [-help] [all PETSc options] */ static char help[] = \"Solves a linear system in parallel with SLES.\\n\\ Input parameters include: \\n\\

-random_exact_so l : use a random exact solution vector\\n\\ -view_exact_sol : write exact solution vector to stdout\\n\\ -m : number of mesh points in x-direction\\n\\ -n : number of mesh points in y-direction\\n\\n\";

/* Include \"petscsles.h\" so that we can use SLES solvers. Note that this file automatically includes:

petsc.h - base PETSc routines petscvec.h - vectors petscsys.h - system routines petscmat.h - matrices

petscis.h - index sets petscksp.h - Krylov subspace methods

petscviewer.h - viewers petscpc.h - preconditioners */

#include \"petscsles.h\" #undef __FUNCT__ #define __FUNCT__ \"main\"

int main(int argc,char **args)

25

{

Vec x,b,u; /* 近似解,右端向量和分析解 */ Mat A; /* 线性算子/矩阵 */ SLES sles; /* 线性解法器 */ PetscRandom rctx; /* 随机数发生器环境 */ PetscReal norm; /* 解误差的范数 */

int i, j, I, J, Istart, Iend, ierr, m = 8, n = 7, its;

PetscTruth flg;

PetscScalar v, one = 1.0, neg_one = -1.0; KSP ksp;

PetscInitialize(&argc, &args, (char *)0, help);

ierr =PetscOptionsGetInt(PETSC_NULL,\"-m\ ierr =PetscOptionsGetInt(PETSC_NULL,\"-n\

/* - - - - - - - - - - - - - - - -- - - - - - - - - - -- - - - - - - - - - -- - - - - - - - - - -- - - - - - - - - - -- - - - - - - - - - - - - - - - - - 计算矩阵算子和右端向量

- - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - - -- */

/* 创建矩阵对象 */

ierr=MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,&A); CHKERRQ(ierr);

ierr=MatSetFromOptions(A); CHKERRQ(ierr);

/*获得局部划分的上下界*/

ierr = MatGetOwnershipRange(A,&Istart,&Iend); CHKERRQ(ierr);

/*给矩阵的每个元素赋值*/ for (I=Istart; Iif (i>0) {J = I - n;

ierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} if (iierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} if (j>0) {J = I - 1;

ierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} if (jierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} v = 4.0;

ierr = MatSetValues(A,1,&I,1,&I,&v,INSERT_VALUES);CHKERRQ(ierr);}

/* 矩阵集聚 */

26

ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

/* 创建向量对象 */

ierr =VecCreate(PETSC_COMM_WORLD,&u); CHKERRQ(ierr); ierr =VecSetSizes(u,PETSC_DECIDE,m*n); CHKERRQ(ierr); ierr =VecSetFromOptions(u); CHKERRQ(ierr); ierr =VecDuplicate(u,&b); CHKERRQ(ierr); ierr =VecDuplicate(b,&x); CHKERRQ(ierr);

/* 设置精确解和右端向量 */

ierr = PetscOptionsHasName(PETSC_NULL,\"-random_exact_sol\if(flg){

ierr = PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rctx); CHKERRQ(ierr); ierr = VecSetRandom(rctx,u); CHKERRQ(ierr);

ierr = PetscRandomDestroy(rctx); CHKERRQ(ierr); } else{

ierr = VecSet(&one,u); CHKERRQ(ierr); }

ierr = MatMult(A,u,b); CHKERRQ(ierr);

/* 屏幕显示精确解 */

ierr = PetscOptionsHasName(PETSC_NULL,\"-view_exact_sol\ if(flg) {ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- - - - - - - - 创建线性算子并设置运行选项

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -- */

/*创建线性解法器环境和设置解算子*/

ierr = SLESCreate(PETSC_COMM_WORLD,&sles);CHKERRQ(ierr);

ierr =SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN); CHKERRQ(ierr);

/*预置运行参数*/

ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr);

ierr=KSPSetTolerances(ksp,1.e-2/((m+1)*(n+1)),1.e-50,PETSC_DEFAULT,PETSC_DEFAULT); CHKERRQ(ierr);

/* 设置运行选项,例如

-ksp_type -pc_type -ksp_monitor -ksp_rtol */

ierr = SLESSetFromOptions(sles); CHKERRQ(ierr);

/*求解线性方程组*/

ierr = SLESSolve(sles,b,x,&its);CHKERRQ(ierr);

27

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 检测误差并释放内存空间

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

/* 误差检测 */

ierr = VecAXPY(&neg_one,u,x);CHKERRQ(ierr); ierr = VecNorm(x,NORM_2,&norm); CHKERRQ(ierr);

/* 打印收敛信息 */

ierr = PetscPrintf(PETSC_COMM_WORLD,\"Norm of error %A iterations %d\\n\

/* 释放存储空间 */

ierr = SLESDestroy(sles); CHKERRQ(ierr); ierr = VecDestroy(u); CHKERRQ(ierr); ierr = VecDestroy(x); CHKERRQ(ierr); ierr = VecDestroy(b); CHKERRQ(ierr); ierr = MatDestroy(A); CHKERRQ(ierr);

/* 退出PETSc运行环境 */

ierr = PetscFinalize(); CHKERRQ(ierr); return 0; }

图4.1 并行PETSc程序范例

4.2 PETSc程序结构

PETSc的所有对象和库都是在MPI消息传递并行环境下开发。所以,一个标准的PETSc程序本身就是一个标准的MPI并行程序。唯一不同的是初始化和结束一个程序的运行需要PETSc来管理。具体地,一个标准的PETSc程序通过引用PetscInitialize函数来初始化或启动整个程序的运行,并通过引用PetscFinalize函数来结束整个程序的运行。

PETSc主要提供了三个核心的解法器,分别针对线性方程组、非线性方程组和时间依赖方程组的求解。而任何一个解法器的使用都包括创建、初始化、执行和释放空间几个阶段。以SLES解法器为例,PETSc程序一般分四个必要阶段:1.创建解算子、解向量和右端向量,并填充解算子(矩阵)、设置初始解和右端向量;2.创建SLES环境,设置Krylov子空间方法和预条件子;3.启动SLES解法器;4.释放所有存储空间

28

5 PETSc范例和应用程序测试

PETSc的范例测试是在曙光2000和深腾6800两个机群系统(关于这两个系统的详细介绍见附录)上完成的。测试分为三个部分:范例可行性测试,范例性能测试和PETSc的应用开发测试。其中性能测试主要包括三个方面:时间代价、浮点性能和计算可扩展性。

在PETSc测试中,我们使用PETSc的-log_summary选项功能来统计性能,同时每个试验重复测试六次,选取一个计算时间代价较小同时附近数据出现频度相对比较高的统计值。这样做的好处在于在机器一般运行状态下可以获得尽可能高的计算可重复性。因此,以下数据仅供参考。

按照PETSc三个主要功能部件(SLES、SNES和TS)的顺序,逐一介绍我们最近关于PETSc版本2.1.6的测试结果,供PETSc用户参考。

5.1 线性方程求解

5.1.1 范例简介

在PETSc版本2.1.6中,总共提供了二十多个SLES范例(C程序),基本覆盖了SLES求解功能的诸多方面,如稀疏线性方程组求解、不同Krylov子空间方法和预条件子的使用、五点格式有限差分PDE方程求解等。这里列举并简单说明如下,以供读者参考。

ex01.c:串行求解一个三对角稀疏线性方程组

ex02.c:并行求解一个线性方程组,二维拉普拉斯方程,5点格式有限差分 ex03.c:并行求解一个线性方程组,矩阵的块方式插入 ex04.c:并行求解一个线性方程组,块对角矩阵,可选预条件子 ex05.c:并行求解两个线性方程组,相同的预条件子和相同的矩阵结构 ex07.c:并行求解一个线性方程组,块雅可比预条件子 ex08.c:并行求解一个线性方程组,ASM预条件子 ex09.c:用不同方法并行求解两个线性方程组

ex10.c:并行求解一个线性方程组,从文件中读入矩阵和向量元素 ex11.c:并行求解一个线性方程组,Helmholtz方程,五点格式有限差分

29

ex12.c:并行求解一个线性方程组, 添加新的预条件子类型,二维拉普拉斯

方程,五点格式有限差分 ex13.c:用SLES串行求解一维Poisson问题

ex15.c:并行求解一个线性方程组,用户预条件子环境的设置 ex16.c:并行求解一系列线性方程组,反复输入一组右端向量 ex22.c:串行求解三维拉普拉斯方程,多重网格预条件子 ex23.c:并行求解一个三对角稀疏线性方程组,同ex1.c ex25.c:串行求解一维拉普拉斯方程,多重网格预条件子 ex26.c:并行求解一维拉普拉斯方程,使用ESI接口

ex27.c:并行求解一个线性方程组,从文件中读入矩阵和向量元素

5.2 非线性方程求解

在PETSc提供的SNES范例和标准测试程序中,绝大多数计算即使在网格规模非常大时其总的计算代价也非常小,同时由于内存的局限也不可能过分通过加细网格来增大计算规模,因此SNES范例的并行可扩展性能非常难于测试。除少数算例(如范例中的EX21、EX24)依赖特定计算环境外,大多范例和全部的标准测试均可正确运行。

5.2.1 范例简介

在PETSc版本2.1.6中,总共提供了二十多个SNES范例,基本覆盖了SNES求解功能的各个方面,这里列举并说明如下,供读者参考。

ex01.c:串行求解一个双变量的非线性方程,牛顿迭代法 ex02.c:串行求解uxxu2f,牛顿迭代法 ex03.c:并行求解uxxu2f,牛顿迭代法

ex05.c:并行求解二维Bratu方程,五点格式有限差分离散,DA实现 ex5s.c:在共享存储环境上(如SGI系统)并行求解二维Bratu方程 ex06.c:串行求解uxxu2f,“无矩阵”牛顿Krylov子空间方法,用户

提供预条件子

30

ex14.c:并行求解三维Bratu方程,七点格式有限差分离散,DA实现 ex18.c:并行求解二维辐射输运方程,多重网格方法,DA实现 ex19.c:并行求解二维驱动空腔流方程,多重网格方法,DA实现 ex20.c:并行求解三维辐射输运方程,多重网格方法,DA实现 ex21.c:串行求解PDE最优化问题

ex22.c:串行求解PDE最优化问题,多重网格方法,DA实现 ex23.c:串行求解PDE问题,中央差分

ex24.c:串行求解PDE最优化问题,多重网格方法,DA实现,伴随方法 ex25.c:并行求解曲面极小化问题,多重网格方法,DA实现 ex26.c:一维CHI平衡问题的Grad-Schafranov求解器

5.3 时间步进积分

由于与SNES相同原因,TS提供的范例和标准测试程序的并行可扩展性测试结果这里也没有给出。除范例中的EX6不能正常编译和标准测试中的EX4不能正常运行外,其余均可正确编译和运行。但是,就并行扩展性能而言,许多TS范例是不成功的。

5.3.1 范例简介

在PETSc版本2.1.6中,总共提供了七个TS范例,基本覆盖了TS求解功能的各个方面,这里列举并说明如下,供读者参考。

ex1.c:用伪时间步进法串行求解Bratu方程

ex2.c:并行求解时间依赖非线性PDE方程,隐式方法 ex3.c:串行求解热传导方程,时间依赖线性PDE方程 ex4.c:并行求解热传导方程,时间依赖线性PDE方程 ex5.c:串行求解热传导方程,时间依赖线性PDE方程 ex6.c:串行求解热传导方程,时间依赖线性PDE方程 ex7.c:并行求解二维时间依赖PDE方程,DA实现

31

附录

1. 曙光2000系统

曙光2000-II系统由82个节点和一个控制台组成。其中72个用作计算节点,配置Motorala MTX双PowerPC604r CPU,333MHz主频;8个节点用作数据库服务,配置IBM43P260双Power3 CPU,200 MHz主频;2个服务器节点一个IBM 43P/260、另一个曙光天演GT355T;1个普通Power604e 300MHz曙光天演GT140工作站,用作系统控制台。系统总计含82个节点,164个CPU,每个节点32KB一级缓存,512KB二级缓存,512M内存,4.5GB 专用硬盘。总共111.7Gflop/s峰值速度,50GB内存,441GB专用硬盘,223GB共享硬盘。

2. 深腾6800系统

深腾6800整体为峰值5万亿次的面向网格的超级计算机系统,含265个四CPU节点,1060个主频为1.3GHz的安腾2处理芯片(1024个用于计算),系统网络为Qsnet,点对点通信带宽大于每秒300MB,延迟时间小于7微秒。每个节点有256KB一级缓存,3MB的二级缓存,8/16GB内存,73GB SCSI硬盘。系统总内存2.6TB,磁盘存储总容量为80TB。LINPACK测试每秒4.183万亿次,效率达到78.5%,在全球TOP500排名第14位。

32

参考文献

[1] Satish Balay, et. al.. PETSc 2.0 users manual. Technical Report ANL-95/11- Revision 2.1.3, Argonne National Laboratory, May 2003. [2] Steve Benson, Lois Curfman McInnes, Jorge J. Moré, and Nason Sarish. Tao Users Manual. Mathematics and Computer Science Division Technical Report ANL/MCS-TM-242-Revision 1.5. 2003. [3] 程强, 迟学斌, 冯仰德等. PETSc:可移植可扩展科学计算工具箱,中科院计算机网络信息中心超级计算中心,2004年7月 [4] 并行软件开发小组,中科院计算机网络信息中心超级计算中心, 超级计算环境基础并行软件平台建设与应用网站:http://www.sccas.cn/~walls/sce

33

因篇幅问题不能全部显示,请点此查看更多更全内容

Copyright © 2019- igat.cn 版权所有

违法及侵权请联系:TEL:199 1889 7713 E-MAIL:2724546146@qq.com

本站由北京市万商天勤律师事务所王兴未律师提供法律服务