并行软件库介绍 赵永华 中国科学院计算机网络信息中心 超级计算中心 yhzhao@sccas.cn 各位评委,老师下午好: 下面我做正高级岗位聘任答辩报告

主要内容 自主并行软件包HPSEPS介绍 MUMPS并行软件包介绍 hypre并行软件包介绍 Parmetis并行软件包介绍 PETSc并行软件包介绍 我竞聘报告的内容包括……4个方面

自主并行软件包HPSEPS 开发者：赵永华　迟学斌等 中国科学院计算机网络信息中心 超级计算中心 任现职以来的研究成果

概述 是目前国际上唯一提供了稠密和稀疏特征问题等多个并行求解器 的并行软件包 ： 用户应用程序可以直接调用所需的并行求解器 HPSEPS（High Performance Symmetric Eigenproblem Software，高性能对称特征问题软件）是由中科院计算机网络信息中心/超级计算中心自主开发。主要用于高性能并行求解实对称/厄密矩阵特征问题、SVD奇异值分解、LSQR分解及其相关问题。 是目前国际上唯一提供了稠密和稀疏特征问题等多个并行求解器 的并行软件包 ： 稠密问题并行求解器优于国际通用软件包ScaLAPACK，并被鉴定专家评为国际领先水平（十五信息化建设专项鉴定专家）； 稀疏问题并行求解器在可扩展性优于国际通用软件包PARPACK； 形成了在千核以上对大规模问题的可扩展并行计算； 已用于多个具体应用问题的并行求解 　用户应用程序可以直接调用所需的并行求解器 　可将软件包的并行模块作为黑匣子嵌入应用程序

应用领域 　　　　大规模特征值问题是许多科学和工程应用：结构动力分析、量子化学、电子结构和材料科学等计算核心，成为国外许多超级计算机极力求解的计算问题之一。 随着计算模型变得越来越复杂，导致的稀疏矩阵规模达到千万阶甚至上亿。 量子计算中电子结构遇到的特征问题达到了千万阶以上规模，有时需要得到几百个特征对； 离子加速器产生的大规模稀疏特征问题可达到上亿阶； 日本地球模拟器上费米-Hubbard模型中遇到的大规模稀疏特征问题规模达到了几百亿阶。 　应用范围广 　是国外许多超级计算机极力求解的计算问题之一。需要扩展到几千处理器规模；具有求解上亿个自由度的计算能力；

HPSEPS功能 稠密对称特征问题并行求解器（基于不同的块算法） 稀疏对称特征问题并行求解器 SVD奇异值分解和LQSR并行求解器 显示重开始＋ deflate技术＋Lanczos 方法 LOPBCG方法 SVD奇异值分解和LQSR并行求解器 并已用于： 中科院化学所的二百万规模“有机共轭体系分子”的求解。 量子点（Quantum Dots）中较大规模稀疏本征问题的并行求解，形成了2048核的可扩展性并行计算； HPSEPS已被安装到中山大学光电材料与技术国家重点实验室的并行计算机系统上。求解了光子晶体中大规模问题，得到了非常好的效果。

软件包架构 HPSEPS的设计采用了自底向上的层次设计方法，每个子程序建立在一些基本例程之上。整个软件包由多层结构构成，其层次结构如下图。 稠密和稀疏矩阵存储分布方式：2D块循环 AIJ BAIJ BDIAG matrix-free 其它 稠密矩阵处理 LANCZOS迭代处理 LOBPCG迭代处理 SVD和PLSQR 分解并行器 求解器 厄密/对称稀疏矩阵 特征问题并行求解器 厄密/对称稠密矩阵特征问题并行求解器 并行QR、LU、LLT分解 HouseHolder并行转化 特征求解器有关的线性代数子程序，预条件：ILU AMG 块-Jacobi等 块-Jacobi等 不同的谱转换：（A-δI）-1 、 M-1A、（A-δM）M-1等 不同的并行正交化、重正交 以及 B-正交化等 管理通信、错误核查和一些辅助子程序 BLAS LAPACK Sparse Matrix-vector MPI OpenMP Multi-core User selects a parallel solver

软件包涉及到的主要算法 提出了各类大规模数值计算和对称特征问题有效并行算法和实现技术。 求解稠密特征问题的不同并行块算法； 求解稀疏特征问题的deflate Lanczos 并行算法； 基于最优化的各类预处理并行算法； 多级混合并行实现技术。 提出了许多有效的并行算法和实现技术：包括 其中，特征值问题求解在科学计算和工程中有广范应用 随着模型规模不断增加，需要将并行计算扩展到数千核，达到求解上亿自由度的计算能力

软件包功能模块接口 下面给出了HPSEPS提供的子程序的简要说明，其中在子程序名中出现的符号*代表z(复双精度)、c（复单精度）、d（双精度）或s（单精度）。 稠密对称矩阵特征问题的子程序主要模块和接口： 第一层是计算对称特征系统问题的一些驱动程序。包括： (1) p *gseps：广义对称/厄密特征问题并行求解（选定的特 征值和特征向量） (2) p*sseps: 标准对称/厄密特征问题并行求解（选定的特征 值和特征向量） 第二层包含特征问题并行求解器所需要的矩阵转换子程序、分解子程序和线性代数子程序等，主要包括： (1) p*syg2st：广义实对称特征问题转化为标准特征问题 (2) p*heg2st：广义Hermitian特征问题转化为标准特征问题 (3) p*trsm：并行计算含有多个右端项的实三角矩阵方程组 (4) p*htrsm：并行计算含有多个右端项的复三角矩阵方程组 (5) p*sytrd：Householder并行转换对称矩阵为三对角形式

(6) p*hetrd：Householder并行转换Hermitian矩阵为三对角形式。 (7) p*stebz：分而治之并行求解实对称三对角矩阵的特征值。 (8) p*steiz：逆迭代并行求解实对称三对角矩阵的特征向量 (9) p*t2s：回代转化并行求解标准特征问题的特征向量 (10) p*st2g：回代转化并行求解广义特征问题的特征向量 最后一层包含HPSEPS内部子程序、通信有关的子程序和一些数据处理、错误检测等管理工具。下面列出了其中部分主要子程序。 与通信有关的子程序： (1) mpi_init：MPI初始化子程序 (2) mpi_creat_cart：创建二维处理器网格通信器 (3) mpi_sub_col：创建一维行通信器 (4) mpi_sub_row：创建一维列通信器 与矩阵分布有关的子程序： (1) mat_2d: 矩阵的二维-块循环分布子程序，得到矩阵的数据 结构和在二维处理器网格上的分布信息

(2) indxg2l：得到存储全局矩阵元素（i, j）的处理器在二 维处理器网格中的逻辑坐标（row_i, col_j） (3) indxg2p：得到全局矩阵元素（i, j）在处理器器上的局 部子矩阵中对应的元素坐标(loc_i, loc_j) (4) mat _ div：三对角矩阵秩-2划分。 (5) *get_sub_mat：得到三对角矩阵秩-2后的子矩阵。 (6) pdist_A: 将矩阵按2D-块方式分布到二维处理器网格中 其它子程序 (1) p*gnrm: 广义特征问题特征向量余范数求解 (2) p*nrm: 标准特征问题特征向量余范数求解 (3) *lag_app_eigen：Laguerre迭代求解函数近似值 (4) *sort：数据排序子程序 (5) Mem_free：释放内存空间 (6) Comm_free：释放通信器

稀疏对称特征问题 在HPSEPS中，提供了基于显式重启-再正交和deflate技术的Lanczos算法的稀疏对称矩阵特征问题并行求解模块。 对于标准稀疏对称特征问题： Ax=λ x　 　 HPSEPS提供了两种求解方式： 　　　标准求解方式 OP=A 　　　 位移逆求解方式 OP=(A-δI )-1

对于广义稀疏对称特征问题： Ax=λMx 　首先将此问题转换为标准特征问题，HPSEPS提　 　　供了四种求解方式： (1) 标准逆方式 OP=M-1A (2) 位移逆方式 OP=(A-δM) -1M 　 为此，HPSEPS提供了不同的用户求解接口，为了给 　　用户使用该软件包提供更好的灵活性，软件包允许用户 　　提供不同的OP操作： 为了保持操作的有效性，矩阵－向量应保持输入向量和输出向量在处理器上分布的一致性。输入向量的第j个元素在处理器P上，输出向量的第j个元素也必须在处理器P上

主要模块和接口： 　　　(1) p*lancs：Lanczos框架接口。其通过调用不同的模块，完成矩阵的三对角分解、正交化处理，得到收敛的Ritz对等。 　 (2) p*getv：产生分布在不同处理器上的初始向量。 　　　(3) p*sletr：m-步Lanczos并行化处理和分解。 　 (4) p*orth：向量并行正交化过程 　　　(5) p*norm2：并行计算向量的2-范数。 　　 另外，针对一般性稀疏矩阵结构，HPSEPS提供了稀疏矩阵-向量积的并行求解模块。

使用HPSEPS编程的方法 HPSEPS为求解不同模式的矩阵特征问题提供了相应的模板。用户通过适当的修改这些模板，可以得到求解具体特征问题的程序。下面是使用HPSEPS软件包应遵循的一些步骤： 选择一个合适的驱动程序。 确定处理器的二维网格结构，分布矩阵到各处理器（稠密问题）。 修改问题依赖的变量。 核查计算结果的精度。

稠密特征问题：在深腾7000超级计算机，使用128，256，512，1024核并行求解 　　　　　　　　　　　　 　30000×30000和60000×60000　　　　　 　规模问题的全部本征对。 稀疏特征问题：求解问题规模大约为190万，得到5个最小本征值。2048个核上性能达到了较高的可扩展性能

劳伦斯-利弗莫尔国家实验室（LLNL）/ Hypre软件包 美国加州大学（UC） 劳伦斯-利弗莫尔国家实验室（LLNL）/ 应用科学计算中心（CASC）

软件包概述 Hypre （High Performance Preconditioners, 高性能预条件子）源于美国能源部和LLNL等在研究国防、环境、能源和生物科学中的物理现象时，开发的一些模拟代码 。主要用于大规模并行计算机上求解大型稀疏线性方程组，目的是为用户提供高级并行预条件子 ，Hypre具有强壮性、易用性、适应性和互动性，其主要特性为： 可扩展的预条件子：包括诸如结构化多重网格（SMG）和代数多重网格（AMG）等几类可扩展求解超大规模稀疏线性方程组的预条件子算法。 常用的迭代法实现： Hypre提供一些最常用的基于Krylov子空间迭代法。比如求解非对称矩阵的GMRES和求解对称矩阵的CG（包括PCG, CGNR, BiCGStab）。 直观的以网格为中心的界面： Hypre通过各种网格界面表示和处理稀疏矩阵，每个界面提供对一些求解器的访问，因此不需要用户去学习和创建复杂的数据结构

Hypre：数据结构、求解器和网格接口关系 第一层表示各种线性系统的网格界面， 第二层表示各种线性求解器（迭代法和预条件子） 第三层表示各种数据划分和矩阵向量存储策略。

网格接口 HYPRE为不同的应用提供了不同的接口, 该接口目前仅支持标量偏微分方程。 结构化（Struct）接口 ： 面向结构网格离散的应用.每个网格点的离散格式具有相同的模式。如有限差分（FDM） 半结构化（Sstruct）接口 面向半结构(semi-struct)网格离散的应用，如局部加密AMG、块结构网格上的应用, 如有限差分方法（FDM）, 有限体积方法（FVM） 基于有限元的无结构界面（FEI） 应用于有限元（FEM）得到的线性方程组 基于线性代数的非结构矩阵界面（IJ） 该接口以矩阵方式显式地表示线性代数方程组,是适用范围最广泛的接口。应用于稀疏线性方程组, 为网格界面的补充

迭代法与预条件子 迭代方法 Krylov解法器（CG，GMRES(缺省情形)，TFQMR，BiCGSTAB）； BoomerAMG（一个并行代数多重网格解法器）； 具有迭代加细(refinement)的SuperLU直接解法器（串行）。 预条件子 Diagonal：对角，块Jacobi预条件子(缺省情形)； PILUT：具有阈值（threshold）的并行不完全LU分解（PILU）； Euclid：并行ILU预条件子的扩展； SMG：半粗化（semi-coarsening）多重网格预条件子；二维和三维情形的光滑子(smoother)分别采用线松弛和面松弛 PFMG：半粗化多重网格预条件子，使用简单点松弛作为光滑子； BoomerAMG：并行代数多重网格（AMG）预条件子；用户可选择不同的并行粗化策略及松驰格式光滑子. ParaSails：并行稀疏近似逆预条件子

ParaSails用于计算优化问题： , 因此M为Frobenius范数下A的近似逆。如果A对称, 且有Cholesky分解：A=LLT, 求解 得到三角近似逆G, ； PILUT并行求解A的一个近似分解。矩阵A支持HYPRE 的ParCSR格式、PETSc 的矩阵形式和ISIS++ Row的矩阵形式。由于M是非对称的（即使A 是对称的），因此不适合作为对称矩阵的迭代法（如CG）的预条件子； Euclid是一种扩展性能较好的并行不完全LU分解（ILU）预条件子，它支持各种ILU（k）和 ILUT, 包括块Jacobi ILU（k），并行ILU（k），它比块Jacobi预条件子更有效

Hypre的网格界面与求解器的关系 X表示支持 网格接口与求解器的关系 HYPRE为不同的接口定义了不同的数据结构,并配以适合该接口的解法器 Hypre的网格界面与求解器的关系 X表示支持

多重网格MG 多重网格解法器是HYPRE的重要特色. 多重网格方法包含三个要素：光滑算子、限制算子和延拓算子 分片线性插值作为延拓 相邻点的加权平均作为限制 松弛迭代(如Gauss-Seidel、SSOR)等简单迭代作为光滑 HYPRE提供多个多重网格解法器 如AMS,SMG,PFMG,MLI, BoomerAMG. 这些可满足各种应用的需求。其中SMG和BoomerAMG是目前实际应用中使用最广泛的两个解法器.

SMG求解矩形网格下对流扩散方程的FDM, FEM或FVM离散得到的方程组。二维时SMG只在x方向半粗化, 在y方向用的是线光滑, 三维时则采用面光滑。而PFMG仅使用简单的点光滑，因此PFMG的健壮性不如SMG，但是它在作V-循环迭代时效率更高. BoomerAMG既可作为迭代法, 也可作为预条件子。用户可以根据情况选择不同的并行粗化技巧（比如CLJP粗化、经典RS粗化）和松弛策略（比如Gauss-Seidel松弛、Jacobi或加权Jacobi松弛）.

Hypre 使用方法 下面的例子是采用并行半粗化多重网格迭代法求解结构网格界面下的线性系统 /* Set up the grid and stencil */ HYPRE_StructGridCreate(MPI_COMM_WORLD, dim, &grid); HYPRE_StructGridSetExtents(grid, ilower, iupper); HYPRE_StructGridAssemble(grid); % 构造结构网格和模板 HYPRE_StructStencilCreate(dim, stencil_size, &stencil); HYPRE_StructStencilSetElement(stencil, 0, offset0); /* Set up the matrix, right-hand side, and initial guess*/ HYPRE_StructMatrixCreate(MPI_COMM_WORLD, grid, stencil, &A); HYPRE_StructMatrixInitialize(A); % 构造结构化矩阵 HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, nelts, elts, Avalues); HYPRE_StructMatrixAssemble(A); HYPRE_StructVectorCreate(MPI_COMM_WORLD, grid, &b); HYPRE_StructVectorInitialize(b); % 右端向量的初始化 HYPRE_StructVectorSetBoxValues(b, ilower, iupper, bvalues); ... HYPRE_StructVectorAssemble(b); HYPRE_StructVectorCreate(MPI_COMM_WORLD, grid, &x); HYPRE_StructVectorInitialize(x); % 解向量的初始化

HYPRE_StructVectorSetBoxValues(x, ilower, iupper, xvalues); HYPRE_StructVectorSetBoxValues(x, ilower, iupper, xvalues); ... HYPRE_StructVectorAssemble(x); /* Set up the solver */ HYPRE_StructPFMGCreate(MPI_COMM_WORLD, &solver); HYPRE_StructPFMGSetMaxIter(solver, 50); /* optional */ HYPRE_StructPFMGSetTol(solver, 1.0e-06); /* optional */ HYPRE_StructPFMGSetup(solver, A, b, x); %创建求解器(PFMG) /* Solve the linear system */ HYPRE_StructPFMGSolve(solver, A, b, x); %求解线性方程组 /* Get solution info and free up memory */ %返回结果并释放内存 HYPRE_StructVectorGetBoxValues(x, ilower, iupper, xvalues); HYPRE_StructPFMGDestroy(solver); HYPRE_StructGridDestroy(grid); HYPRE_StructStencilDestroy(stencil); HYPRE_StructMatrixDestroy(A); HYPRE_StructVectorDestroy(b); HYPRE_StructVectorDestroy(x);

算例 对流-反应-扩散方程 对流-反应-扩散(Convection-Reaction-Diffusion)方程: div (-K grad u + B u) + C u = F in Ω, 采用五点差分离散, 得到方程组: Au = b, 其中A = [Aii Aib ; Abi Abb], u = [ui ; ub], b = [bi ; u0]。 考虑到边界条件 u= u0 on ∂Ω, 即 ub= u0 . 于是 [Aii 0 ; 0 I] [ui ; ub] = [bi - Aibu0 ; u0] 。 在hypre-1.10.0b/src/examples/ex4.c中, Ω为单位正方形, 处理机网格为N ×N，每个处理机上的网格为n ×n, h=1/(Nn+1), 采用结构网格界面和5点差分离散, 并考虑边界条件。 相关系数: 扩散系数K=x2+exp(y); 对流系数B=1.0; 反应系数 C=10.0; 边界条件: u0 =(sin(5πx)+sin(5πy))/1000; 右端项: F= 2π2 sin(πx)sin(πy) 。 ex4.c提供4种迭代法：SMG、PFMG、PCG、GMRES, 后两种迭代法可以增加SMG、PFMG, 对角或块Jacobi等4种预条件子。

Struct和IJ两种界面下, 各求解器的迭代次数和运行时间(256×256,Np=4 ) 结构化网格界面 矩阵界面 求解器 SMG PCG + SMG AMG PCG + AMG PCG + Parasails CG 迭代次数 9 6 7 209 437 计算时间(s) 2.24 1.93 2.16 2.34 11.52 4.66 注: T(SMG) = T(SMG_setup) + T(SMG_solve), (cpu clock time), Np=4, tol=1.e-6 AMG的并行效率(网格规模为1024×1024)

由 CEC ESPRIT IV长期研究计划项目资助 MUMPS 由 CEC ESPRIT IV长期研究计划项目资助

MUMPS概述 MUMPS：多波前大规模并行稀疏直接解法器（A MUltifrontal Massively Parallel sparse direct Solver） MUMPS是一个通过直接方法求解线性方程组: Ax=b 的并行软件包，其中A是一个对称或非对称的稀疏方阵。 MUMPS基于多波前方法的直接求解方法。通过将矩阵A直接分解为A=LU或A=LDLT（对称矩阵）形式完成大规模线性方程的求解。

主要功能 求解不同类型的稀疏矩阵方程问题： 对称或非对称矩阵（部分主元法），复和实算术矩阵 提供了多种矩阵输入格式： 组装格式（assembled format） 分布式组装格式（distributed assembled format） 单元格式（elemental format） 迭代加密和向前误差分析; 部分分解和Schur补矩阵 提供了多个排序接口：AMD，AMF，PORD，METIS和SCOTCH

输入矩阵 矩阵类型 矩阵类型在初始化阶段（JOB=-1）由所有进程通过参数mumps par%SYM设定： 0: A是非对称型 1：A是对称正定型 3：A是一般对称矩阵 矩阵的输入格式 MUMPS提供了多种矩阵输入格式，这些由参数ICNTL(5)和ICNTL(18)控制。 单元格式： 矩阵由主进程(host)集中输入， 置 ICNTL(5)=1 ICNTL(18)=0 组装格式（assembled format） 矩阵由主进程(host)集中输入 结构由主进程提供 (analysis) 元素被分布到各处理器上 (numeric factorization)

主要计算步 MUMPS计算Ax=b通过三步完成: (1)分析(JOB=1) 主进程执行排序操作 主进程执行符号分解 (2) A=LU或A=LDLT分解(JOB=2) A被分布到各处理器 由主进程和一个或多个从进程对每个波前矩阵进行数值分解 (3) 求解(JOB=3) b由主进程分布到各处理器 x由分布到各处理器的因子计算得到 x被聚集到主进程或分布到各处理器

主要特性 每个处理阶段可独立调用 异步通信 使得计算和通信实现了重叠 动态调度 算法是自适应的，在执行时重分布任务和数据到适当的处理器

MUMPS应用 MUMPS用户包括学术界和工业界，目前用户数已超过1000个。 并且平均每天被下载一次。 典型应用包括： 结构力学和CAD 流体力学、磁流体和物理化学 地震波传播与成像，海洋模式 声学和电磁学传播 有限元分析，数值优化，仿真 航空、地球物理等

OSSAU code of French CEA/CESTA：2D / 3D structural mechanics code ODYSSEE code of French CEA/CESTA 已成功应用到具有3千万未知量的3D网格问题 Electro-magnetism code (Finite Element Meth. + Integral Equation) Complex double precision, Schur Compl. Fluid mechanics LU factorization with static pivoting (SuperLU approach like) Car body 148770 unknowns and 5396386 nonzeros MSC.Software

MUMPS使用方法及实例 组装格式程序 下面是使用MUMPS计算双精度问题的组装格式例子程序。两个文件mpi.h和mumps_struc.h必须被包含在程序中，MPI的初始化和终止通过调用MPI_INIT和MPI_FINALIZE完成。在程序中，首先设定 JOB=1对初始化MUMPS，由主进程读入求解的问题（N, NZ, IRN, JCN, A, 和HS）。通过设定JOB=6，然后调用MUMPS由所有进程完成问题的求解。最后设定JOB=-2，调用MUMPS完成数据结构的释放 PROGRAM MUMPS_EXAMPLE NCLUDE ’mpif.h’ INCLUDE ’dmumps_struc.h’ TYPE (DMUMPS_STRUC) id INTEGER IERR, I CALL MPI_INIT(IERR) C Define a communicator for the package 　 id%COMM = MPI_COMM_WORLD C Ask for unsymmetric code 　 id%SYM = 0 C Host working 　 id%PAR = 1 C Initialize an instance of the package 　 id%JOB = -1 　 CALL DMUMPS(id)

这样，对一个构造好的矩阵和右端向量： ， 我们给出如下输入： 5 ：N 12 : NZ 1 2 3.0 2 3 -3.0 4 3 2.0 C Define problem on the host (processor 0) 　 IF ( id%MYID .eq. 0 ) THEN 　　 READ(5,*) id%N 　　 READ(5,*) id%NZ 　　 ALLOCATE( id%IRN ( id%NZ ) ) 　　 ALLOCATE( id%JCN ( id%NZ ) ) 　　 ALLOCATE( id%A( id%NZ ) ) 　　 ALLOCATE( id%RHS ( id%N ) ) 　　 READ(5,*) ( id%IRN(I) ,I=1, id%NZ ) 　　 READ(5,*) ( id%JCN(I) ,I=1, id%NZ ) 　　 READ(5,*) ( id%A(I),I=1, id%NZ ) 　　 READ(5,*) ( id%RHS(I) ,I=1, id%N ) 　END IF C Call package for solution 　id%JOB = 6 　CALL DMUMPS(id) C Solution has been assembled on the host 　IF ( id%MYID .eq. 0 ) THEN 　　WRITE( 6, * ) ’ Solution is ’,(id%RHS(I),I=1,id%N) C Deallocate user data 　IF ( id%MYID .eq. 0 )THEN 　　DEALLOCATE( id%IRN ) 　　DEALLOCATE( id%JCN ) 　　DEALLOCATE( id%A ) 　　DEALLOCATE( id%RHS ) C Destroy the instance (deallocate internal data structures) 　id%JOB = -2 　CALL MPI_FINALIZE(IERR) 　STOP END 这样，对一个构造好的矩阵和右端向量： ，　　 　　我们给出如下输入： 5 　　　　　　　　　　　：N 12 　　　　　　　　　　: NZ 1 2 3.0 2 3 -3.0 4 3 2.0 5 5 1.0 2 1 3.0 1 1 2.0 5 2 4.0 3 4 2.0 2 5 6.0 3 2 -1.0 1 3 4.0 3 3 1.0 　　　　　　　　: A 20.0 24.0 9.0 6.0 13.0 :RHS　 我们将得到解RHS(i) = i, i = 1, . . . , 5.

Lawrence Livermore National Laboratory ParMETIS并行软件 Lawrence Livermore National Laboratory

ParMETIS ParMETIS：并行图划分和填充-约化矩阵排序(Parallel Graph Partitioning and Fill-reducing Matrix Ordering) ， 特别适合于大规模无结构网格的并行数值模拟。 ParMETIS基于MPI并行库，实现了用于无结构图划分、网格划分、计算稀疏矩阵的填充-约化次序等多种算法。 ParMETIS扩展了METIS所提供的功能并包含了特别适合于并行计算和大规模数值模拟的子程序。 ParMETIS中实现的算法基于并行多层k-路图划分算法、自适应再划分算法及并行多约束算法。

功能描述 图划分 快速计算非常大规模图的高质量划分； 利用几何信息 (当可用时) 降低划分时间而不损失质量； 可为多相及多物理计算划分图。 网格划分 直接计算非常大规模网格高质量划分, 无需应用程序创建基本图； 提供网格对偶图的高效并行程序。 图重划分 快速计算自适应加密网格的高质量再划分； 优化移去的顶点个数以及所得划分的边切割。 划分加细 改进由其它划分算法产生的划分的质量。 矩阵重排序 计算稀疏矩阵的填充-约化(fill-reducing)次序； 使用基于节点的嵌套剖分算法，此算法显示比其它流行重排序算法更优越。

子程序调用 ParMetis可以执行下列操作 ParMETIS_V3_PartKway ParMETIS_V3_PartGeom 无结构图划分 是否存在顶点坐标 ParMETIS_V3_PartKway 你有什么时间/质量权衡 ParMETISV3_PartGeomKway ParMETIS_V3_PartGeom ParMETIS_V3_Mesh2Dual ParMETIS_V3_PartMeshKway ParMETIS_V3_AdaptiveRepart ParMETIS_V3_RefineKway ParMETIS_V3_NodeND 网格划分 由网格构造图 自适应加密重划分图 精化划分质量 计算填充-约化次序

无结构图划分 图划分的并行子程序ParMETIS_V3_PartKway基于串行多层k-路分区算法。该算法已被证明能够迅速生成高质量的划分。它包括三个阶段： 图的粗化、初步划分、加密。 当顶点坐标可用时，PARMETIS也为非结构图划分提供了ParMETIS_PartGeom函数。 ParMETIS _PartGeom仅基于空间填充曲线方法计算一个分区。因此，它非常快（通常比ParMETIS PartGeomKway快5到10倍），但它的计算质量差。此子程序可用于空间填充曲线适合的分区技术中的某些计算（例如，N体计算）。

下面为一个无结构图的划分过程：

直接网格划分 PARMETIS3.0提供了一个新的例程ParMETIS_V3_PartMeshKway支持由网格（而不是图）作为输入的划分和重划分计算。 在内部，ParMETIS_V3的PartMeshKway使用mesh_to-graph子程 序，然后调用由ParMETIS_V3 PartKway和ParMETIS_V3 PartGeomKway使用的同样核心划分子程序 。 PARMETIS没有提供直接计算网格自适应重划分这样的例程。然而，它提供了例程ParMETIS_V3_Mesh2Dual，用来快速、并行地为一个给定的网格构造一个偶图（dual graph）。它可以用来为ParMETIS_V3_AdaptiveRepart子程序构建一个输入图形。 从本质上讲，ParMETIS_V3_PartMeshKway和ParMETIS_V3_Mesh2Dual承担着用户高效地编写一个网格到图的子程序责任。实验表明，这个例行通常占用了PARMETIS计算划分时约一半的运行时间

自适应加密网格 PARMETIS提供了重划分自适应加密网格的子程序ParMETIS_V3_ AdaptiveRepart。 该子程序假设网格已很好地分布在各处理器，但存在着不好的负载平衡。

划分加密 PARMETIS提供了用来改善已存在划分质量的子程序ParMETIS_V3_ RefineKway。一旦一个图被划分并被重新分布，可以调用ParMETIS_V3_ RefineKway,进一步改善划分的质量。 像ParMETIS_V3_AdaptiveRepart，该子程序假设图已很好的被分布在各处理器。 同ParMETIS_V3_AdaptiveRepart一样, ParMETIS_V3_RefineKway执行局部粗化。这两个子程序也使用用样加密算法。不同之处在于ParMETIS_V3_RefineKway没用像ParMETIS_V3_AdaptiveRepart一样对最粗的图进行初始平衡化分解。二是假设图已很好地分布，并且初始划分有好的平衡。这样对于已很好分布的图而言，ParMETIS_V3_RefineKway是很快的。

填充-约化(fill-reducing)次序 ParMETIS_V3_NodeND是PARMETIS提供的计算填充-约化次序的子程序。适合于基于Cholesky的直接分解算法。ParMETIS_V3_NodeND对图初始如何分布在各处理器没用要求。它能有效地随机分布图的填充-约化次序。 ParMETIS_V3_NodeND首先使用ParMETIS V3 PartKway计算高质量的划分并相应的进行图的重新分布。接下来计算logp层消去树。当图已经被分成P部分（P是处理器数），图被重新分布在各处理器。这样每个处理器接收唯一的子图，而一个多最小度算法用来对这些更小的子图进行排序。

Parmetis输入输出格式 在Parmetis中，所有与图有关的子程序的输入格式包括：图的邻接结构、顶点和边的权重、描述图如何被分布在各处理器上的数组。 图得结构采用压缩存储格式（CSR）,我们首先为串行图描述CSR存储结构，然后描述图如何被分布到各处理器上。下面为一个结构图的存储例图, 图(a)为一个简单的图，图(b)为串行CSR存储格式，而图(C)是分布式CSR格式。

测试结果 测试数据： 在深腾7000上对上述数据在1024核上所测试的结果： Graph scheme 32核 128核 256核 512核 Auto WF LMSR URA 2.35 2.31 2.42 1.25 1.16 1.20 0.62 0.73 0.37 0.58 0.55 0.35 0.45 0.49 mdual2 3.29 3.25 3.21 1.61 1.56 2.54 0.80 0.87 0.48 0.50 0.53 0.41 0.42 0.56 mrng3 11.11 11.13 11.32 5.55 5.54 5.71 2.88 2.86 3.06 1.47 1.64 0.93 0.96 1.11

PETSC：并行可扩展科学计算工具箱 Argonne 国家实验室

PETSC 越来越多的应用程序在PETSc环境上开发，并逐渐显示出PETSc在高效求解大规模数值模拟问题方面的优势和威力 　　　PETSC：并行可扩展科学计算工具箱（Parallel Extensible Toolkits for Scientific Computing） 国能源部ODE2000 支持开发的20 多个ACTS 工具箱之一，由Argonne 国家实验室开发的可移植可扩展科学计算工具箱，主要用于在分布式存储环境高效求解偏微分方程组及相关问题。PETSc所有消息传递通信均采用MPI标准实现。 提供大量面向对象的并行代数数据结构、解法器和相关辅助部件，适合并行可扩展求解PDE方程（有限差分、有限元、有限体离散的隐式和显示格式）。 基于MPI、BLAS库、LAPACK库 使用Fortran、C/C++开发 　　　越来越多的应用程序在PETSc环境上开发，并逐渐显示出PETSc在高效求解大规模数值模拟问题方面的优势和威力

PETSc：并行可扩展科学计算工具箱（Parallel Extensible Toolkits for Scientific Computing） 核心人员：美国数学与计算机部、Argonne国家重点实验室等等 基于MPI、BLAS库、LAPACK库 使用Fortran、C/C++开发 PETSc软件包含一个功能强大的工具集以在高性能计算机上数值求解偏微分方程及其相关问题 可移植性：CRAY T3D，T3E，Origin 2000, IBM SP, HP UX, ASCI Red, Blue Mountain, NOWs，LINUX，ALPHA等 公开源代码，免费下载 http://www.mcs.anl.gov/petsc

PETSc 的一些模块处理： 索引集，包括用于向量索引的置换，重新计数等 向量 矩阵（一般是稀疏的） 分布阵列（对正规的基于网格问题的并行化有用） Krylov 子空间方法 预条件子，包括多重网格和稀疏直接解法器 非线性解法器 解时间相关（非线性）PDEs 的时间步进解法器

体系结构 PETSc为用户提供了一个通用的层次化应用程序开发平台。基于PETSc提供的大量对象和解法库，用户可以灵活地开发自己的应用程序。 PDE解法器 PDE解法器 PDE解法器 PDE解法器 SNES(无约束优化、非线性解法器) SNES(无约束优化、非线性解法器) SNES(无约束优化、非线性解法器) SNES(无约束优化、非线性解法器) SNES(无约束优化、非线性解法器) SNES(无约束优化、非线性解法器) SNES(无约束优化、非线性解法器) SLES(线性方程解法器) SLES(线性方程解法器) SLES(线性方程解法器) SLES(线性方程解法器) SLES(线性方程解法器) SLES(线性方程解法器) SLES(线性方程解法器) SLES(线性方程解法器) KSP(Krylov子空间方法) KSP(Krylov子空间方法) KSP(Krylov子空间方法) KSP(Krylov子空间方法) KSP(Krylov子空间方法) KSP(Krylov子空间方法) KSP(Krylov子空间方法) KSP(Krylov子空间方法) KSP(Krylov子空间方法) KSP(Krylov子空间方法) KSP(Krylov子空间方法) DRAW DRAW DRAW DRAW DRAW DRAW DRAW DRAW DRAW PC(预条件) PC(预条件) PC(预条件) PC(预条件) PC(预条件) PC(预条件) PC(预条件) PC(预条件) PC(预条件) PC(预条件) 矩阵 矩阵 矩阵 矩阵 矩阵 矩阵 矩阵 矩阵 矩阵 矩阵 矩阵 矩阵 向量 向量 向量 向量 向量 向量 向量 向量 向量 向量 向量 向量 向量 向量 索引集 索引集 索引集 索引集 索引集 索引集 索引集 索引集 索引集 索引集 索引集 索引集 索引集 BLAS BLAS BLAS BLAS BLAS BLAS BLAS BLAS BLAS BLAS BLAS BLAS BLAS BLAS BLAS BLAS LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK LAPACK MPI MPI MPI MPI MPI MPI MPI MPI MPI MPI MPI MPI MPI MPI MPI MPI PETSc实现的层次结构

PETSc的数值组件 向量 非线性解法器 时间步法 Euler 方法 向后Euler方法 拟时间步 其它 Krylov 子空间方法 预条件子 牛顿迭代法 其它 线搜索 信赖域 时间步法 Euler 方法 向后Euler方法 拟时间步 其它 Krylov 子空间方法 GMR CG CGS Bi-CG-Sta TFQMR Richardson Chebyshev 其它 预条件子 加法Schwarz 块 Jacobi Jacobi ILU ICC LU 其它 向量 压缩稀疏行(AIJ) 块压缩稀疏(BAIJ) 块对角(BDiag) 稠密 其它 向量 索引集 索引 块索引 跨度 其它

PETSc的基本对象 向量 　　　　向量是最简单的PETSc对象。PETSc向量对象主要用于存储线性方程组的解和右端向量。PETSc 提供AO 对象来管理向量在全局和局部之间的索引、排序和映射。PETSc 还提供了两个对象DA和IS，来分别管理向量在规则正交网格和无结构网格上各进程之间的分发、聚集和边界点的数据通信等操作。 矩阵 　　　矩阵是PETSc的基本对象。PETSc同时提供了稠密矩阵和稀疏行矩阵的基本运算功能，以及一些特殊格式（如“无矩阵”实现，无结构网格划分等内容）和用户提供的某些功能扩展和实现。PETSc 的矩阵运算和操作主要包括矩阵的创建、插值、聚集、各种算术运算和释放。PETSc 的各种矩阵运算和操作使用起来非常方便，用户无需关心矩阵的具体存储实现。

PETSc的核心组件 PETSc的三个核心组件包括：线性方程求解器（SLES）、非线性方程求解器（SNES）和时间步进积分器（TS）。 　　　　构成了PETSc最核心的部分。它不仅是几乎所有PDE方程求解器的基本内核，而且也是实现PETSc 的其它两个核心组件SNES 和TS 的必不可少的部分。SLES求解线性方程组 Ax = b 　　　　　 　其中解算子A是n＊ｎ维非奇异矩阵，b是n维右端向量，x为n维解向量。 　　　　 SLES求解包括：线性方程求解环境的创建、Krylov子空间方法和预条件子（PC）的选择、收敛性判据、LU直接求解等。 其本用法： SLESCreate：创建一个线性方程求解环境 SLESSetOperators：设置求解算子（矩阵） SLESSetFromoptions：通过运行参数设置SLES运行选项 SLESSolve：启动一个线性方程求解器 SLESDestroy：释放一个线性方程求解环境 SLESSetup：启动一个线性方程求解器 SLESGetPC：获得PC对象/环境的访问权 SLESGetKSP：获得KSP对象/环境的访问权

非线性方程求解 　　　　SNES非线性解法器基于牛顿迭代法（线性搜索和信赖域方法），依赖线性解法器SLES实现。雅可比矩阵的求解是SNES解法器的重要组成部分。SNES求解以下形式的非线性方程组 　　　　　　　　　F(x) = 0 　　其中解算子F : Rn → Rn。 基本用法 SNESCreate：创建一个非线性方程求解环境 SNESSetType：设置非线性求解器的类型 SNESSetFromOptions：通过运行参数设置SNES运行选项 SNESSolve：启动一个非线性方程求解器 SNESDestroy：释放一个非线性方程求解器 SNESSetFunction：设置非线性函数 SNESSetJacobian：设置雅可比矩阵

时间步进积分 　　　TS 时间步进积分器，用于求解依赖时间或时间演化的ODE 方程，或依赖时间的离散化后的PDE方程。TS主要求解如下时间依赖问题 u t = F(u,t) 　其中u为有限维解向量，上式通常为运用有限差分或有限元方法离散后的常微分方线性程组。对于非时间演化或稳态方程，PETSc提供了伪时间步进积分器。TS积分器最终依赖线性解法器SLES和非线性解法器SNES来实现。PETSc为PVODE求解器提供了接口。 基本用法 TSCreate：创建一个TS求解环境 TSSetType：设置TS求解器的类型 TSSetInitialTimeStep：设置初始时间和步长 TSSetTimeStep：设置时间步长 TSGetTimeStep：获得时间步长 TSSetDuration：设置最大时间步数 TSSetUp：启动TS求解环境 TSDestroy：释放TS求解环境 TSView：启动TS屏幕输出

PETSc与其它软件 PETSc 可扩展性的另一个方面表现在其为非常广泛的一类数值软件和数学库提供了很方便的软件接口。主要包括以下几种类型： 线性代数求解器，如AMG、BlockSolve95、DSCPACK、hypre、ILUTP、LUSOL、SPAI、SPOOLES、SuperLU、SuperLU_Dist； 最优化软件，如TAO、Veltisto； 离散化和网格生成和优化工具包，如Overture、SAMRAI、SUMAA3d； 常微分方程求解器，如PVODE； 其它，如Matlab、ParMETIS。

PETSc编程 　　　　PETSc与MPI消息传递并行环境完全兼容。在这个意义上，用户可以在PETSc上开发任何基于消息传递的应用程序。但是PETSc 更希望用户将其视为一个在分布式计算环境中的PDE 数值模拟和科学计算的平台，用户基于PETSc 提供的大量线性方程和非线性方程求解器、丰富的数值迭代方法和各种预条件子，大量的对象和库资源，以及软件接口来开发和调试应用程序。总之，PETSc 给用户提供了丰富的算法资源和可编程环境。对于SNES 和TS求解器，用户通常还需提供计算函数及其雅可比矩阵的子程序。 　　　下面是一个典型的PETSc 程序范例，它使用有限差分来求解二维Laplace 问题。用SLES（Krylov 子空间方法和预条件子）求解线性方程组。

/*$Id: ex2.c,v 1.94 2001/08/07 21:30:54 bsmith Exp $*/ /* 运行方式: mpirun -np <procs> ex2 [-help] [all PETSc options] */ static char help[] = "Solves a linear system in parallel with SLES.\n\ Input parameters include: \n\ -random_exact_sol : use a random exact solution vector\n\ -view_exact_sol : write exact solution vector to stdout\n\ -m <mesh_x> : number of mesh points in x-direction\n\ -n <mesh_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) { 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",&m,PETSC_NULL); CHKERRQ(ierr); ierr =PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL); CHKERRQ(ierr);

26 /* - - - - - - - - - - - - - - - -- - - - - - - - - - -- - - - - - - - - - -- - - - - - - - - - -- - - - - - - - - - -- - 计算矩阵算子和右端向量 - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - - - - -- - - - - - - - - - - - - - - - - -- */ /*创建矩阵对象*/ 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; I<Iend; I++) { v = -1.0; i = I/n; j = I - i*n; if (i>0) {J = I - n; ierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} if (i<m-1) {J = I + n; if (j>0) {J = I - 1; if (j<n-1) {J = I + 1; v = 4.0; ierr = MatSetValues(A,1,&I,1,&I,&v,INSERT_VALUES);CHKERRQ(ierr);} /*矩阵集聚*/ 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",&flg); CHKERRQ(ierr); 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); } 27 ierr = MatMult(A,u,b); CHKERRQ(ierr); /*屏幕显示精确解*/ ierr = PetscOptionsHasName(PETSC_NULL,"-view_exact_sol",&flg); CHKERRQ(ierr); 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 <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>*/ ierr = SLESSetFromOptions(sles); CHKERRQ(ierr); /*求解线性方程组*/ ierr = SLESSolve(sles,b,x,&its);CHKERRQ(ierr);

/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 检测误差并释放内存空间 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ /*误差检测*/ 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",norm,its); CHKERRQ(ierr); /*释放存储空间*/ 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); 28 /*退出PETSc运行环境*/ ierr = PetscFinalize(); CHKERRQ(ierr); return 0; }