ParMetis 网格分割笔记

Table of Contents

• ParMetis
  • ParMetis简介
  • ParMetis算法
  • ParMetis存储格式
  • ParMetis的API

ParMetis

ParMetis简介

ParMetis是一个MPI并行库. 它内置了许多非结构图和网格分割(partionting)和细分(reparationing)的算法. ParMetis特别适合处理大规模数值模拟中的非结构网格划分问题. 计算流体力学就是一个典型的例子, 流体问题的模拟往往需要成千上万个的集群节点同时计算. 在计算前首先要进行计算网格进行区域分解, 计算中不同区域间需要信息交流. 因此网格区域分解对系统负载均衡和进程间通信较大的影响. PartMetis能使分割的区域网格大致相同并且交界面个数尽可能最少, 使交换的数据量最小化, 进而显著的减少进程间通信的时间.
ParMetis主要有以下功能:

1. 非结构图和网格分割;
2. 自适应网格的再分割;
3. 多相和多物理网格的分割;
4. 提高已有分割的质量;
5. 构造网格的对偶图(Dual Graph).

ParMetis使用和Metis中相同的分割算法, 可以看做是后者的并行版本. 除了Metis的基本功能外, ParMetis增加了并行的模块, 使其特别适合并行计算, 尤其是大规模数值计算.

ParMetis算法

ParMetis提供以一系列分割的模块, 用户可以根据任务, 选择合适的模块进行调用. 下图是一个ParMetis的主要模块.
在本文中主要介绍最常用的两个: 图划分和网格划分.

ParMetis存储格式

ParMetis提供C/C++和Fortran的调用接口, 用户只需要按照接口的要求调用函数就可以实现图和网格的分割功能.

头文件

任何掉用ParMetis的程序必须首先包含 parmetis.h 头文件. 这个头文件提供了所有API, 常数和数据类型的定义.

在ParMetis中, 有两个重要的数据类型, 分别是 idx_t 和 real_t . 这两个类型分别是用来存储整型和浮点型变量. indx_t 可以是32位或64位有符号整数, ~real_t~可以是单精度或双精度浮点数. 所有API的输入都是这两个类型的数组.

图和网格的存储

下面介绍如何ParMetis中图和网格的存储, 图和网格的存储方式类似, 这里主要介绍图存储的CSR格式.

图(Graph)

对于图来说, 要保存如下的结构. 首先需要一个数组存储所有的点, 还需要另一个数组用来放点之间的边. 如果点和边有权重的话, 令需要数列存储这些数据.
在Metis中, 图的结构使用CSR(Compressed Sparse Graphs)格式存储, 在ParMetis中将其拓展到分布式存储结构中. 下面先介绍串行的CSR, 再介绍相应的并行版本.

串行的CSR格式

CSR格式是一个广泛使用的存储稀疏图的格式. 在这个格式中, 图的相邻结构通过两个数组 xadj 和 adjncy 存储. 如果节点和边存在权重的话, 可以使用 vwget 和 adjwgt 表示. 例如, 一个图有 个节点 个边. 那么用如下的几个数组表示这个图:

xadj[n+1], vwgt[n], adjncy[2m], adjwgt[2m]

这里需要注意 adjncy 和 adjwgt 都有 个元素, 这是因为每个边都被统计了两次(比如 和 ). 如果点和边都没有权重, 数组 vwgt 和 adjwgt 都为空.

图的相邻结构按照如下的方式存储. 首先假设节点的从0开始计数. 第 个节点的相邻节点存储在 adjncy 数组中, 从第 adj[i] 个下标到 xadj[i+1] (但步包括 xadj[i+1] ). 也就是说,从 adjncy[xadj[i]] 到 adjncy[adj[i+1]-1] . 因此每个节点的相邻节点都存在 adjncy 中. xadj 数组的作用是指出节点列的起始和结束位置.下图是一个15个节点的图, 下面的表格是对象CSR存储结构.

  xadj  0 2 5 8 11 13 16 20 24 28 32 33 36 39 42 44
adjncy  1 5 0 2 6 1 2 7 2 4 8 3 9 0 6 10  1 5 7 11 2 6 8 12 3 7 9 13 4 8 14 5 11 6 10 12 7 11 13 8 12 14 9 13

同时可以发现:
1. adjncy中是和每个点相邻的点.
2. 这些点是顺时针排列的.

并行的CSR格式

ParMetis对CSR格式进行了拓展,使其能够存储分布在不同处理器上的相邻结构.基本原理是这样的,对于 个节点, 个进程,每个进程都存储 个点. 每个进程都有相应的 xadj[n_i+1] , vwgt[n_i] , adjncy[m_i] 和 adjwgt[m_i] . 如果图是没有权重的,最简单的方法是给每个进程分配 个节点. 除了以上的数组, 还需要另外一个数组- vtxdisp[p+1] . 每个进程都有这个数组,进程 储存的节点从 vtxdist[i] 到 vtxdist[i]-1 .
下图表示一个三个进程的系统. 图(a)中的15个节点被平均分配个3个进程,每个节点有5个节点. 进程0获得了0-4号节点,进程1或得到5-9;进程3获得了10-14号.图中分别给出了每个进程上的各类数组.

Processor 0:  xadj  0 2 5 8 11 13
            adjncy  1 5 0 2 6 1 3 7 2 4 8 3 9
           vtxdist  0 5 10 15

Processor 1:  xadj  0 3 7 11 15 18
            adjncy  0 6 10 1 5 7 11 2 6 8 12 3 7 9 13 4 8 14 
           vtxdist  0 5 10 15

Processor 2:  xadj  0 2 5 8 11 13
            adjncy  5 11 6 10 12 7 10 12 7 11 13 8 12 14 9 13
           vtxdist  0 5 10 15

网格(Mesh)

在网格存储格式中,有三个重要的数列, eptr , eind 和 elmdist ,它们的作用和图(graph)中的 xadj , adjncy 和 vtxdist 类似. 如果有 个单元的话,则有 eptr 有 个元素. 通过 epr 可以确定单元的点在 eind 中的分布. 比如为了获得第 个单元的的顶点列表. 首先获得 eptr[i] 和 eptr[i+1] 这两个元素, 然后节点列表就是 eind[eptr[i]] 到 eint[eptr[i+1]-1] 的点. 这个和图的存储格式类似.

ParMetis的API

有了上面的知识, 下面介绍一个具体调用函数. 这里使用的是 ParMETIS_V3_PartMeshKway 函数. 先看该函数的完整接口.

    int ParMETIS_V3_PartMeshKway(
            indx_t *elmdist, indx_t *eptr, indx_t *eind, indx_t *elmwgt, indx_t *wgtflag, indx_t numflag,
            indx_t *ncon, indx_t *ncommomnodes, indx_t *nparts, real_t *tpwgts, real_t *ubvec,
            intx_t * otions, indx_t *edgcut, indx_t *part, MPI_Comm *comm
            )

1. 函数描述
   将一个在 个进程上的网格划分成 个区域. 网格可以包含不同的类型.
2. 主要参数
   
   • elmdist
     描述网格单元是如何在进程上分布的. 每个进程包含相同的elmdist.
   • eptr, eind
     见上节网格的描述.
   • elmwgt
     每个单元的权重.
   • wgtflag
     wgtflag有两个值, 0表示没有权重, 2表示只给节点权重.
   • numflag
     数组的计算格式, 0表示C语言格式, 也就是从0开始; 1表示FORTRAN的格式, 从1开始.
   • ncon
     描述每个节点的权重.
   • ncommonnodes
     每个节点在dual graph的联通度.
   • nparts
     网格分成多少个区域.
   • tpwgts
     一个 的数组, 用来确定每个节点权重的分数. 如果每个区域单元的个数一样. 那么数组每个元素都为.
   • part
     这个数组用来标记划分成功后的节点.
   • comm
     这是MPI的通信子.
3. 返回值
   
   • METIS_OK
     函数正常返回.
   • METIS_ERROR
     出现错误.