《最新并行软件库介绍PPT课件》由会员分享，可在线阅读，更多相关《最新并行软件库介绍PPT课件（68页珍藏版）》请在金锄头文库上搜索。
1、并行软件库介绍主要内容主要内容 l自主并行软件包HPSEPS介绍lMUMPS并行软件包介绍lhypre并行软件包介绍lParmetis并行软件包介绍lPETSc并行软件包介绍2软件包功能模块接口软件包功能模块接口 下面给出了HPSEPS提供的子程序的简要说明，其中在子程序名中出现的符号*代表z(复双精度)、c（复单精度）、d（双精度）或s（单精度）。稠密对称矩阵特征问题的子程序主要模块和接口：稠密对称矩阵特征问题的子程序主要模块和接口：u第一层是计算对称特征系统问题的一些驱动程序。包括： (1) p *gseps：广义对称/厄密特征问题并行求解（选定的特 征值和特征向量） (2) p*ssep
2、s: 标准对称/厄密特征问题并行求解（选定的特征 值和特征向量）u第二层包含特征问题并行求解器所需要的矩阵转换子程序、分解子程序和线性代数子程序等，主要包括： (1) p*syg2st：广义实对称特征问题转化为标准特征问题 (2) p*heg2st：广义Hermitian特征问题转化为标准特征问题 (3) p*trsm：并行计算含有多个右端项的实三角矩阵方程组 (4) p*htrsm：并行计算含有多个右端项的复三角矩阵方程组 (5) p*sytrd：Householder并行转换对称矩阵为三对角形式9 (6) p*hetrd：Householder并行转换Hermitian矩阵为三对角形式。
3、(7) p*stebz：分而治之并行求解实对称三对角矩阵的特征值。 (8) p*steiz：逆迭代并行求解实对称三对角矩阵的特征向量 (9) p*t2s：回代转化并行求解标准特征问题的特征向量 (10) p*st2g：回代转化并行求解广义特征问题的特征向量u最后一层包含HPSEPS内部子程序、通信有关的子程序和一些数据处理、错误检测等管理工具。下面列出了其中部分主要子程序。与通信有关的子程序： (1) mpi_init：MPI初始化子程序 (2) mpi_creat_cart：创建二维处理器网格通信器 (3) mpi_sub_col：创建一维行通信器 (4) mpi_sub_row：创建一维列
4、通信器与矩阵分布有关的子程序： (1) mat_2d: 矩阵的二维-块循环分布子程序，得到矩阵的数据 结构和在二维处理器网格上的分布信息10 (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-块方式分布到二维处理器网格中u其它子程序 (1)
5、p*gnrm: 广义特征问题特征向量余范数求解 (2) p*nrm: 标准特征问题特征向量余范数求解 (3) *lag_app_eigen：Laguerre迭代求解函数近似值 (4) *sort：数据排序子程序 (5) Mem_free：释放内存空间 (6) Comm_free：释放通信器11稀疏对称特征问题稀疏对称特征问题在HPSEPS中，提供了基于显式重启-再正交和deflate技术的Lanczos算法的稀疏对称矩阵特征问题并行求解模块。HPSEPS为标准稀疏对称特征问题和广义稀疏对称标准特征问题求解提供了不同的方法和接口:p对于标准稀疏对称特征问题： Ax= x HPSEPS提供了两种求
6、解方式： 标准求解方式 OP=A 位移逆求解方式 OP=(A-I )-1 12p对于广义稀疏对称特征问题： Ax=Mx 首先将此问题转换为标准特征问题，HPSEPS提供了四种求解方式： (1) 标准逆方式 OP=M-1A (2) 位移逆方式 OP=(A-M) -1M 为此，HPSEPS提供了不同的用户求解接口，为了给用户使用该软件包提供更好的灵活性，软件包允许用户提供不同的OP操作： 为了保持操作的有效性，矩阵向量应保持输入向量和输出向量在处理器上分布的一致性。输入向量的第j个元素在处理器P上，输出向量的第j个元素也必须在处理器P上13主要模块和接口：主要模块和接口：(1) p*lancs：L
7、anczos框架接口。其通过调用不同的模块，完成矩阵的三对角分解、正交化处理，得到收敛的Ritz对等。 (2) p*getv：产生分布在不同处理器上的初始向量。(3) p*sletr：m-步Lanczos并行化处理和分解。 (4) p*orth：向量并行正交化过程(5) p*norm2：并行计算向量的2-范数。 另外，针对一般性稀疏矩阵结构，HPSEPS提供了稀疏矩阵-向量积的并行求解模块。14使用使用HPSEPS编程的方法编程的方法 HPSEPS为求解不同模式的矩阵特征问题提供了相应的模板。用户通过适当的修改这些模板，可以得到求解具体特征问题的程序。下面是使用HPSEPS软件包应遵循的一些步
8、骤：选择一个合适的驱动程序。确定处理器的二维网格结构，分布矩阵到各处理器（稠密问题）。修改问题依赖的变量。核查计算结果的精度。15 稠密特征问题：稠密特征问题：在深腾7000超级计算机，使用128，256，512，1024核并行求解 3000030000和6000060000 规模问题的全部本征对。稀疏特征问题：稀疏特征问题：求解问题规模大约为190万，得到5个最小本征值。2048个核上性能达到了较高的可扩展性能16Hypre软件包软件包 美国加州大学（UC）劳伦斯-利弗莫尔国家实验室（LLNL）/应用科学计算中心（CASC）17软件包概述软件包概述 Hypre （High Performan
9、ce Preconditioners, 高性能预条件子）源于美国能源部和LLNL等在研究国防、环境、能源和生物科学中的物理现象时，开发的一些模拟代码 。主要用于大规模并行计算机上求解大型稀疏线性方程组，目的是为用户提供高级并行预条件子 ，Hypre具有强壮性、易用性、适应性和互动性，其主要特性为：n可扩展的预条件子：可扩展的预条件子：包括诸如结构化多重网格（SMG）和代数多重网格（AMG）等几类可扩展求解超大规模稀疏线性方程组的预条件子算法。n常用的迭代法实现：常用的迭代法实现：Hypre提供一些最常用的基于Krylov子空间迭代法。比如求解非对称矩阵的GMRES和求解对称矩阵的CG（包括PC
10、G, CGNR, BiCGStab）。n直观的以网格为中心的界面：直观的以网格为中心的界面：Hypre通过各种网格界面表示和处理稀疏矩阵，每个界面提供对一些求解器的访问，因此不需要用户去学习和创建复杂的数据结构18Hypre：数据结构、求解器和网格接口关系 第一层表示各种线性系统的网格界面，第二层表示各种线性求解器（迭代法和预条件子）第三层表示各种数据划分和矩阵向量存储策略。 19网格接口 HYPRE为不同的应用提供了不同的接口, 该接口目前仅支持标量偏微分方程。n 结构化结构化（Struct）接口接口 ： 面向结构网格离散的应用.每个网格点的离散格式具有相同的模式。如有限差分（FDM）n 半
11、结构化半结构化（Sstruct）接口 面向半结构(semi-struct)网格离散的应用，如局部加密AMG、块结构网格上的应用, 如有限差分方法（FDM）, 有限体积方法（FVM）n 基于有限元的无结构界面基于有限元的无结构界面（FEI） 应用于有限元（FEM）得到的线性方程组n 基于线性代数的非结构矩阵界面基于线性代数的非结构矩阵界面（IJ） 该接口以矩阵方式显式地表示线性代数方程组,是适用范围最广泛的接口。应用于稀疏线性方程组, 为网格界面的补充20迭代法与预条件子迭代法与预条件子n迭代方法迭代方法Krylov解法器（CG，GMRES(缺省情形)，TFQMR，BiCGSTAB）；Boome
12、rAMG（一个并行代数多重网格解法器）；具有迭代加细(refinement)的SuperLU直接解法器（串行）。n预条件子预条件子Diagonal：对角，块Jacobi预条件子(缺省情形)；PILUT：具有阈值（threshold）的并行不完全LU分解（PILU）；Euclid：并行ILU预条件子的扩展；SMG：半粗化（semi-coarsening）多重网格预条件子；二维和三维情形的光滑子(smoother)分别采用线松弛和面松弛PFMG：半粗化多重网格预条件子，使用简单点松弛作为光滑子；BoomerAMG：并行代数多重网格（AMG）预条件子；用户可选择不同的并行粗化策略及松驰格式光滑子.
13、ParaSails：并行稀疏近似逆预条件子21ParaSails用于计算优化问题： , 因此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），它比
14、块Jacobi预条件子更有效22网格接口与求解器的关系 Hypre的网格界面与求解器的关系的网格界面与求解器的关系X表示支持表示支持 HYPRE为不同的接口定义了不同的数据结构,并配以适合该接口的解法器23多重网格MGn多重网格解法器是HYPRE的重要特色.n多重网格方法包含三个要素多重网格方法包含三个要素：光滑算子、限制算子和延拓算子：光滑算子、限制算子和延拓算子分片线性插值作为延拓相邻点的加权平均作为限制松弛迭代(如Gauss-Seidel、SSOR)等简单迭代作为光滑nHYPRE提供多个多重网格解法器 如AMS,SMG,PFMG,MLI, BoomerAMG. 这些可满足各种应用的需求。
15、其中SMG和BoomerAMG是目前实际应用中使用最广泛的两个解法器. 24SMG求解矩形网格下对流扩散方程的FDM, FEM或FVM离散得到的方程组。二维时SMG只在x方向半粗化, 在y方向用的是线光滑, 三维时则采用面光滑。而PFMG仅使用简单的点光滑，因此PFMG的健壮性不如SMG，但是它在作V-循环迭代时效率更高.BoomerAMG既可作为迭代法, 也可作为预条件子。用户可以根据情况选择不同的并行粗化技巧（比如CLJP粗化、经典RS粗化）和松弛策略（比如Gauss-Seidel松弛、Jacobi或加权Jacobi松弛）.25Hypre使用方法使用方法下面的例子是采用并行半粗化多重网格迭
16、代法求解结构网格界面下的线性系统 /* 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, offset
17、0);/* 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_C
18、OMM_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); % 解向量的初始化 26HYPRE_StructVectorSetBoxValues(x, ilower, iu
19、pper, 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 l
20、inear 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_StructMatrixDest
21、roy(A); HYPRE_StructVectorDestroy(b); HYPRE_StructVectorDestroy(x);27算例p对流对流-反应反应-扩散方程扩散方程对流-反应-扩散(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
22、- 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(5x)+sin(5y)/1000; 右端项: F= 22 sin(x)sin(y) 。ex4.c提供4种迭代法：SMG、PFMG、PCG、GMRES, 后两种迭代法可以增加SMG、PFMG, 对角或块Jacobi等4种预条件子。28St
23、ruct和IJ两种界面下, 各求解器的迭代次数和运行时间(256256,Np=4 )界面结构化网格界面矩阵界面求解器SMGPCG + SMGAMGPCG + AMGPCG + ParasailsCG迭代次数9697209437计算时间(s)2.241.932.162.3411.524.66注注:T(SMG)=T(SMG_setup)+T(SMG_solve),(cpuclocktime),Np=4,tol=1.e-6AMG的并行效率的并行效率(网格规模为网格规模为10241024) 29 MUMPS由 CEC ESPRIT IV长期研究计划项目资助30MUMPS概述 MUMPS：多波前大规模并
24、行稀疏直接解法器（A MUltifrontal Massively Parallel sparse direct Solver）nMUMPS是一个通过直接方法求解线性方程组: Ax=b 的并行软件包，其中A是一个对称或非对称的稀疏方阵。nMUMPS基于多波前方法的直接求解方法。通过将矩阵A直接分解为A=LU或A=LDLT（对称矩阵）形式完成大规模线性方程的求解。31主要功能l求解不同类型的稀疏矩阵方程问题： 对称或非对称矩阵（部分主元法），复和实算术矩阵 l提供了多种矩阵输入格式：组装格式（assembled format）分布式组装格式（distributed assembled forma
25、t）单元格式（elemental format）l迭代加密和向前误差分析; l部分分解和Schur补矩阵l提供了多个排序接口：AMD，AMF，PORD，METIS和SCOTCH32输入矩阵输入矩阵 n矩阵类型矩阵类型 矩阵类型在初始化阶段（JOB=-1）由所有进程通过参数mumps par%SYM设定： 0: A是非对称型 1：A是对称正定型 3：A是一般对称矩阵n矩阵的输入格式矩阵的输入格式 MUMPS提供了多种矩阵输入格式，这些由参数ICNTL(5)和ICNTL(18)控制。l单元格式：单元格式： 矩阵由主进程(host)集中输入， 置 ICNTL(5)=1 ICNTL(18)=0l组装格
26、式（组装格式（assembledformat）矩阵由主进程(host)集中输入结构由主进程提供 (analysis)元素被分布到各处理器上 (numeric factorization)33主要计算步主要计算步n MUMPS计算Ax=b通过三步完成:(1)分析(JOB=1) 主进程执行排序操作 主进程执行符号分解(2) A=LU或A=LDLT分解(JOB=2) A被分布到各处理器 由主进程和一个或多个从进程对每个波前矩阵进行数值分解(3) 求解(JOB=3) b由主进程分布到各处理器 x由分布到各处理器的因子计算得到 x被聚集到主进程或分布到各处理器34主要特性n每个处理阶段可独立调用n异步通
27、信 使得计算和通信实现了重叠n动态调度 算法是自适应的，在执行时重分布任务和数据到适当的处理器35MUMPS应用应用MUMPS用户包括学术界和工业界，目前用户数已超过1000个。 并且平均每天被下载一次。典型应用包括：结构力学和CAD流体力学、磁流体和物理化学地震波传播与成像，海洋模式声学和电磁学传播有限元分析，数值优化，仿真航空、地球物理等36 nOSSAU code of French CEA/CESTA：2D / 3D structural mechanics codenODYSSEE code of French CEA/CESTA 已成功应用到具有3千万未知量的3D网格问题lElec
28、tro-magnetism code (Finite Element Meth. + Integral Equation)lComplex double precision, Schur Compl.nFluid mechanicslLU factorization with static pivoting (SuperLU approach like)nCarbody148770unknownsand5396386nonzerosMSC.Software37MUMPS使用方法及实例使用方法及实例 p组装格式程序组装格式程序 下面是使用MUMPS计算双精度问题的组装格式例子程序。两个文件mpi
29、.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(IE
30、RR)C Define a communicator for the package id%COMM = MPI_COMM_WORLDC Ask for unsymmetric code id%SYM = 0C Host working id%PAR = 1C Initialize an instance of the package id%JOB = -1 CALL DMUMPS(id)38C Define problem on the host (processor 0) IF ( id%MYID .eq. 0 ) THEN READ(5,*) id%N READ(5,*) id%NZ A
31、LLOCATE( 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 IFC Call package for solutionid%JOB = 6CALL DM
32、UMPS(id)C Solution has been assembled on the hostIF ( id%MYID .eq. 0 ) THENWRITE( 6, * ) Solution is ,(id%RHS(I),I=1,id%N)END IFC Deallocate user dataIF ( id%MYID .eq. 0 )THENDEALLOCATE( id%IRN )DEALLOCATE( id%JCN )DEALLOCATE( id%A )DEALLOCATE( id%RHS )END IFC Destroy the instance (deallocate intern
33、al data structures)id%JOB = -2CALL DMUMPS(id)CALL MPI_FINALIZE(IERR)STOPEND这样，对一个构造好的矩阵和右端向量：，我们给出如下输入：5 ：N12 : NZ1 2 3.02 3 -3.04 3 2.05 5 1.02 1 3.01 1 2.05 2 4.03 4 2.02 5 6.03 2 -1.01 3 4.03 3 1.0 : A20.024.09.06.013.0 :RHS我们将得到解RHS(i) = i, i = 1, . . . , 5.39ParMETIS并行软件并行软件LawrenceLivermoreNat
34、ionalLaboratory 40ParMETIS ParMETIS：并行图划分和填充-约化矩阵排序(Parallel Graph Partitioning and Fill-reducing Matrix Ordering) ， 特别适合于大规模无结构网格的并行数值模拟。nParMETIS基于MPI并行库，实现了用于无结构图划分、网格划分、计算稀疏矩阵的填充-约化次序等多种算法。nParMETIS扩展了METIS所提供的功能并包含了特别适合于并行计算和大规模数值模拟的子程序。nParMETIS中实现的算法基于并行多层k-路图划分算法、自适应再划分算法及并行多约束算法。 41功能描述图划分
35、快速计算非常大规模图的高质量划分；利用几何信息 (当可用时) 降低划分时间而不损失质量； 可为多相及多物理计算划分图。 网格划分 直接计算非常大规模网格高质量划分, 无需应用程序创建基本图； 提供网格对偶图的高效并行程序。图重划分快速计算自适应加密网格的高质量再划分；优化移去的顶点个数以及所得划分的边切割。划分加细 改进由其它划分算法产生的划分的质量。 矩阵重排序 计算稀疏矩阵的填充-约化(fill-reducing)次序； 使用基于节点的嵌套剖分算法，此算法显示比其它流行重排序算法更优越。 42子程序调用子程序调用 ParMetis可以执行下列操作可以执行下列操作无结构图划分是否存在顶点坐标
36、ParMETIS_V3_PartKway你有什么时间/质量权衡ParMETISV3_PartGeomKwayParMETIS_V3_PartGeomParMETIS_V3_Mesh2DualParMETIS_V3_PartMeshKwayParMETIS_V3_AdaptiveRepartParMETIS_V3_RefineKwayParMETIS_V3_NodeND网格划分由网格构造图自适应加密重划分图精化划分质量计算填充-约化次序43无结构图划分无结构图划分 图划分的并行子程序ParMETIS_V3_PartKway基于串行多层k-路分区算法。该算法已被证明能够迅速生成高质量的划分。它包括
37、三个阶段： 图的粗化、初步划分、加密。当顶点坐标可用时，PARMETIS也为非结构图划分提供了ParMETIS_PartGeom函数。 ParMETIS _PartGeom仅基于空间填充曲线方法计算一个分区。因此，它非常快（通常比ParMETIS PartGeomKway快5到10倍），但它的计算质量差。此子程序可用于空间填充曲线适合的分区技术中的某些计算（例如，N体计算）。 44下面为一个无结构图的划分过程：45直接网格划分直接网格划分 PARMETIS3.0提供了一个新的例程ParMETIS_V3_PartMeshKway支持由网格（而不是图）作为输入的划分和重划分计算。l 在内部，Par
38、METIS_V3的PartMeshKway使用mesh_to-graph子程 序，然后调用由ParMETIS_V3 PartKway和ParMETIS_V3 PartGeomKway使用的同样核心划分子程序 。l PARMETIS没有提供直接计算网格自适应重划分这样的例程。然而，它提供了例程ParMETIS_V3_Mesh2Dual，用来快速、并行地为一个给定的网格构造一个偶图（dual graph）。它可以用来为ParMETIS_V3_AdaptiveRepart子程序构建一个输入图形。 l 从本质上讲，ParMETIS_V3_PartMeshKway和ParMETIS_V3_Mesh2Du
39、al承担着用户高效地编写一个网格到图的子程序责任。实验表明，这个例行通常占用了PARMETIS计算划分时约一半的运行时间 46自适应加密网格自适应加密网格 PARMETIS提供了重划分自适应加密网格的子程序ParMETIS_V3_ AdaptiveRepart。 该子程序假设网格已很好地分布在各处理器，但存在着不好的负载平衡。47划分加密划分加密 PARMETIS提供了用来改善已存在划分质量的子程序ParMETIS_V3_ RefineKway。一旦一个图被划分并被重新分布，可以调用ParMETIS_V3_ RefineKway,进一步改善划分的质量。l 像ParMETIS_V3_Adapti
40、veRepart，该子程序假设图已很好的被分布在各处理器。l同ParMETIS_V3_AdaptiveRepart一样, ParMETIS_V3_RefineKway执行局部粗化。这两个子程序也使用用样加密算法。不同之处在于ParMETIS_V3_RefineKway没用像ParMETIS_V3_AdaptiveRepart一样对最粗的图进行初始平衡化分解。二是假设图已很好地分布，并且初始划分有好的平衡。这样对于已很好分布的图而言，ParMETIS_V3_RefineKway是很快的。48填充填充-约化约化(fill-reducing)次序次序 ParMETIS_V3_NodeND是PARME
41、TIS提供的计算填充-约化次序的子程序。适合于基于Cholesky的直接分解算法。ParMETIS_V3_NodeND对图初始如何分布在各处理器没用要求。它能有效地随机分布图的填充-约化次序。 ParMETIS_V3_NodeND首先使用ParMETIS V3 PartKway计算高质量的划分并相应的进行图的重新分布。接下来计算logp层消去树。当图已经被分成P部分（P是处理器数），图被重新分布在各处理器。这样每个处理器接收唯一的子图，而一个多最小度算法用来对这些更小的子图进行排序。49Parmetis输入输出格式输入输出格式在Parmetis中，所有与图有关的子程序的输入格式包括：图的邻接结
42、构、顶点和边的权重、描述图如何被分布在各处理器上的数组。图得结构采用压缩存储格式（CSR）,我们首先为串行图描述CSR存储结构，然后描述图如何被分布到各处理器上。下面为一个结构图的存储例图, 图(a)为一个简单的图，图(b)为串行CSR存储格式，而图(C)是分布式CSR格式。50测试结果测试结果 测试数据： 在深腾7000上对上述数据在1024核上所测试的结果： Graphscheme32核128核256核512核1024核AutoAutoAutoWFLMSRURA2.352.312.421.251.161.200.620.620.730.370.580.550.350.450.49mdual
43、2mdual2mdual2WFLMSRURA3.293.253.211.611.562.540.800.800.870.480.500.530.410.420.56mrng3mrng3mrng3WFLMSRURA11.1111.1311.325.555.545.712.882.863.061.471.471.640.930.961.1151PETSC：并行可扩展科学计算工具箱Argonne 国家实验室52PETSC PETSC：并行可扩展科学计算工具箱（Parallel Extensible Toolkits for Scientific Computing）n国能源部ODE2000 支持开发
44、的20 多个ACTS 工具箱之一，由Argonne 国家实验室开发的可移植可扩展科学计算工具箱，主要用于在分布式存储环境高效求解偏微分方程组及相关问题。PETSc所有消息传递通信均采用MPI标准实现。n提供大量面向对象的并行代数数据结构、解法器和相关辅助部件，适合并行可扩展求解PDE方程（有限差分、有限元、有限体离散的隐式和显示格式）。n基于MPI、BLAS库、LAPACK库n使用Fortran、C/C+开发越来越多的应用程序在PETSc环境上开发，并逐渐显示出PETSc在高效求解大规模数值模拟问题方面的优势和威力53PETSc：并行可扩展科学计算工具箱（Parallel Extensible
45、 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 54PETSc的一些模块处理：索引集，包括用于向量索引的置
46、换，重新计数等向量矩阵（一般是稀疏的）分布阵列（对正规的基于网格问题的并行化有用）Krylov 子空间方法预条件子，包括多重网格和稀疏直接解法器非线性解法器解时间相关（非线性）PDEs 的时间步进解法器55体系结构体系结构PETSc为用户提供了一个通用的层次化应用程序开发平台。基于PETSc提供的大量对象和解法库，用户可以灵活地开发自己的应用程序。PDE解法器SNES(无约束优化、非线性解法器)SLES(线性方程解法器)TS时间步进矩阵向量索引集KSP(Krylov子空间方法)PC(预条件)DRAWBLASLAPACKMPI应用程序PDE解法器应用程序PDE解法器应用程序PDE解法器LAPAC
47、KMPIBLASLAPACKMPIBLASLAPACKMPI向量BLASLAPACKMPI索引集向量BLASLAPACKMPI矩阵索引集向量BLASLAPACKMPIKSP(Krylov子空间方法)矩阵索引集向量BLASLAPACKMPIPC(预条件)KSP(Krylov子空间方法)矩阵索引集向量BLASLAPACKMPIDRAWPC(预条件)KSP(Krylov子空间方法)矩阵索引集向量BLASLAPACKMPISLES(线性方程解法器)DRAWPC(预条件)KSP(Krylov子空间方法)矩阵索引集向量BLASLAPACKMPISNES(无约束优化、非线性解法器)SLES(线性方程解法器)
48、DRAWPC(预条件)KSP(Krylov子空间方法)矩阵索引集向量BLASLAPACKMPISNES(无约束优化、非线性解法器)SLES(线性方程解法器)DRAWPC(预条件)KSP(Krylov子空间方法)矩阵索引集向量BLASLAPACKMPISNES(无约束优化、非线性解法器)SLES(线性方程解法器)DRAWPC(预条件)KSP(Krylov子空间方法)矩阵索引集向量BLASLAPACKMPITS时间步进SNES(无约束优化、非线性解法器)SLES(线性方程解法器)DRAWPC(预条件)KSP(Krylov子空间方法)矩阵索引集向量BLASLAPACKMPI应用程序TS时间步进SNE
49、S(无约束优化、非线性解法器)SLES(线性方程解法器)DRAWPC(预条件)KSP(Krylov子空间方法)矩阵索引集向量BLASLAPACKMPITS时间步进SNES(无约束优化、非线性解法器)SLES(线性方程解法器)DRAWPC(预条件)KSP(Krylov子空间方法)矩阵索引集向量BLASLAPACKPETSc实现的层次结构实现的层次结构56PETSc的数值组件非线性解法器非线性解法器牛顿迭代法其它线搜索信赖域时间步法时间步法Euler方法向后Euler方法拟时间步其它Krylov子空间方法子空间方法GMRCGCGSBi-CG-StaTFQMRRichardsonChebyshev其
50、它预条件子预条件子加法Schwarz块 JacobiJacobiILUICCLU其它向量向量压缩稀疏行(AIJ)块压缩稀疏(BAIJ)块对角(BDiag)稠密其它索引集索引集索引块索引跨度其它向量57PETSc的基本对象的基本对象n向量向量向量是最简单的PETSc对象。PETSc向量对象主要用于存储线性方程组的解和右端向量。PETSc 提供AO 对象来管理向量在全局和局部之间的索引、排序和映射。PETSc 还提供了两个对象DA和IS，来分别管理向量在规则正交网格和无结构网格上各进程之间的分发、聚集和边界点的数据通信等操作。n矩阵矩阵矩阵是PETSc的基本对象。PETSc同时提供了稠密矩阵和稀疏
51、行矩阵的基本运算功能，以及一些特殊格式（如“无矩阵”实现，无结构网格划分等内容）和用户提供的某些功能扩展和实现。PETSc 的矩阵运算和操作主要包括矩阵的创建、插值、聚集、各种算术运算和释放。PETSc 的各种矩阵运算和操作使用起来非常方便，用户无需关心矩阵的具体存储实现。58PETSc的核心组件PETSc的三个核心组件包括：线性方程求解器（SLES）、非线性方程求解器（SNES）和时间步进积分器（TS）。n线性方程求解线性方程求解构成了PETSc最核心的部分。它不仅是几乎所有PDE方程求解器的基本内核，而且也是实现PETSc 的其它两个核心组件SNES 和TS 的必不可少的部分。SLES求解
52、线性方程组Ax = b 其中解算子A是n维非奇异矩阵，b是n维右端向量，x为n维解向量。 SLES求解包括：线性方程求解环境的创建、Krylov子空间方法和预条件子（PC）的选择、收敛性判据、LU直接求解等。l其本用法：SLESCreate：创建一个线性方程求解环境SLESSetOperators：设置求解算子（矩阵）SLESSetFromoptions：通过运行参数设置SLES运行选项SLESSolve：启动一个线性方程求解器SLESDestroy：释放一个线性方程求解环境SLESSetup：启动一个线性方程求解器SLESGetPC：获得PC对象/环境的访问权SLESGetKSP：获得KSP
53、对象/环境的访问权59n非线性方程求解非线性方程求解SNES非线性解法器基于牛顿迭代法（线性搜索和信赖域方法），依赖线性解法器SLES实现。雅可比矩阵的求解是SNES解法器的重要组成部分。SNES求解以下形式的非线性方程组F(x) = 0其中解算子F : Rn Rn。l基本用法基本用法SNESCreate：创建一个非线性方程求解环境SNESSetType：设置非线性求解器的类型SNESSetFromOptions：通过运行参数设置SNES运行选项SNESSolve：启动一个非线性方程求解器SNESDestroy：释放一个非线性方程求解器SNESSetFunction：设置非线性函数SNESSe
54、tJacobian：设置雅可比矩阵60n时间步进积分时间步进积分TS 时间步进积分器，用于求解依赖时间或时间演化的ODE 方程，或依赖时间的离散化后的PDE方程。TS主要求解如下时间依赖问题u t = F(u,t)其中u为有限维解向量，上式通常为运用有限差分或有限元方法离散后的常微分方线性程组。对于非时间演化或稳态方程，PETSc提供了伪时间步进积分器。TS积分器最终依赖线性解法器SLES和非线性解法器SNES来实现。PETSc为PVODE求解器提供了接口。l基本用法基本用法TSCreate：创建一个TS求解环境TSSetType：设置TS求解器的类型TSSetInitialTimeStep：
55、设置初始时间和步长TSSetTimeStep：设置时间步长TSGetTimeStep：获得时间步长TSSetDuration：设置最大时间步数TSSetUp：启动TS求解环境TSDestroy：释放TS求解环境TSView：启动TS屏幕输出61PETSc与其它软件与其它软件PETSc 可扩展性的另一个方面表现在其为非常广泛的一类数值软件和数学库提供了很方便的软件接口。主要包括以下几种类型：线性代数求解器，如AMG、BlockSolve95、DSCPACK、hypre、ILUTP、LUSOL、SPAI、SPOOLES、SuperLU、SuperLU_Dist；最优化软件，如TAO、Veltist
56、o；离散化和网格生成和优化工具包，如Overture、SAMRAI、SUMAA3d；常微分方程求解器，如PVODE；其它，如Matlab、ParMETIS。62PETSc编程编程PETSc与MPI消息传递并行环境完全兼容。在这个意义上，用户可以在PETSc上开发任何基于消息传递的应用程序。但是PETSc 更希望用户将其视为一个在分布式计算环境中的PDE 数值模拟和科学计算的平台，用户基于PETSc 提供的大量线性方程和非线性方程求解器、丰富的数值迭代方法和各种预条件子，大量的对象和库资源，以及软件接口来开发和调试应用程序。总之，PETSc 给用户提供了丰富的算法资源和可编程环境。对于SNES
57、和TS求解器，用户通常还需提供计算函数及其雅可比矩阵的子程序。下面是一个典型的PETSc 程序范例，它使用有限差分来求解二维Laplace 问题。用SLES（Krylov 子空间方法和预条件子）求解线性方程组。63/*$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.nInput parameters include:
58、 n-random_exact_sol : use a random exact solution vectorn-view_exact_sol : write exact solution vector to stdoutn-m : number of mesh points in x-directionn-n : number of mesh points in y-directionnn;/* Include petscsles.h so that we can use SLES solvers. Note that this file automatically includes:pe
59、tsc.h - base PETSc routines petscvec.h - vectorspetscsys.h - system routines petscmat.h - matricespetscis.h - index sets petscksp.h - Krylov subspace methodspetscviewer.h - viewers petscpc.h - preconditioners */#include petscsles.h#undef _FUNCT_#define _FUNCT_ mainint main(int argc,char *args)Vec x,
60、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 =PetscOpti
61、onsGetInt(PETSC_NULL,-m,&m,PETSC_NULL); CHKERRQ(ierr);ierr =PetscOptionsGetInt(PETSC_NULL,-n,&n,PETSC_NULL); CHKERRQ(ierr);6426/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -计算矩阵算子和右端向量- - - - - - - - - - - - - - - - - - - - - - -
62、- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */*创建矩阵对象*/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); CH
63、KERRQ(ierr);/*给矩阵的每个元素赋值*/for (I=Istart; I0) J = I - n;ierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);if (i0) J = I - 1;ierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);if (jn-1) J = I + 1;ierr = MatSetValues(A,1,&I,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);v = 4.0;ierr =
64、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);ie
65、rr =VecSetFromOptions(u); CHKERRQ(ierr);ierr =VecDuplicate(u,&b); CHKERRQ(ierr);ierr =VecDuplicate(b,&x); CHKERRQ(ierr);65/*设置精确解和右端向量*/ierr = PetscOptionsHasName(PETSC_NULL,-random_exact_sol,&flg); CHKERRQ(ierr);if(flg)ierr = PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rctx); CHKERRQ(ierr);i
66、err = VecSetRandom(rctx,u); CHKERRQ(ierr);ierr = PetscRandomDestroy(rctx); CHKERRQ(ierr); elseierr = VecSet(&one,u); CHKERRQ(ierr); 27ierr = MatMult(A,u,b); CHKERRQ(ierr);/*屏幕显示精确解*/ierr = PetscOptionsHasName(PETSC_NULL,-view_exact_sol,&flg); CHKERRQ(ierr);if(flg) ierr = VecView(u,PETSC_VIEWER_STDOU
67、T_WORLD);CHKERRQ(ierr);/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -创建线性算子并设置运行选项创建线性算子并设置运行选项- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68、 - - - - - - - - - - - - - - - - - - - - - - - - - - - */*创建线性解法器环境和设置解算子*/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+
69、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);66/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
70、 - - - - - - - - - - - - - -检测误差并释放内存空间检测误差并释放内存空间- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */*误差检测*/ierr = VecAXPY(&neg_one,u,x);CHKERRQ(ierr);ierr = VecNorm(x,NORM_2,&norm); CHKERRQ(ierr);/*打印收敛信息*/ierr = PetscPrintf(PETSC_COMM_
71、WORLD,Norm of error %A iterations %dn,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;67结束语结束语谢谢大家聆听！谢谢大家聆听！68
