3 图论图论在计算机科学、信息科学、人工智能、网络理论、系统工程、控制论、运筹学和经济管理等领域有着广泛的应用。

但很多图论问题虽易表达，却难以求解，其中有相当多的图论问题均属NP完全问题。

本章主要介绍工程实用简单图论问题的并行算法及其MPI编程实现，包括传递闭包、连通分量、最短路径和最小生成树等。

1.1 传递闭包设A是一个含N个顶点的有向图G的布尔邻接矩阵(Boolean Adjacent Matrix)，即元素a ij=1当且仅当从顶点i到j有一条有向边。

所谓A的传递闭包（Transitive Closure），记为A+，是一个N×N的布尔矩阵，其元素b ij=1当且仅当：①i=j；或②从i出发存在有向路径到j，又称为顶点i到j可达。

从G的布尔邻接矩阵A求出G的传递闭包，就称为传递闭包问题。

传递闭包问题有很强的应用背景，在科学计算中广泛存在。

传递闭包问题的经典解法之一就是利用布尔矩阵的乘法来求解。

本节将把这一算法分别在串行和并行机上实现。

1.1.1 传递闭包串行算法利用布尔矩阵相乘求解传递闭包问题的原理为：对于布尔矩阵(A+I)k中的任一元素b ij，b ij=1表示从i到j存在长度小于等于k的可达路径，否则，b ij=0。

显然对于k=1，(A+I)1中b ij=1当且仅当从i到j路径长度为0（i=j）或为1（从i到j存在有向边）；(A+I)2中，b ij=1当且仅当从i到j路径长度小于等于2；((A+I)2) 2中，b ij=1当且仅当从i到j路径长度小于等于4，等等。

因为任意两点间如果存在可达路径，长度最多为N-1，所以k≥N-1时，(A+I)k 就是所求的传递闭包A+。

于是(A+I)矩阵的㏒N次自乘之后所得的矩阵就是所求的传递闭包。

根据前面的叙述，很自然的有下面的传递闭包串行算法15.1，其时间复杂度为O(N3㏒N)。

算法15.1传递闭包串行算法输入：图G的布尔邻接矩阵A输出：A的传递闭包Mprocedure closureBegin(1)读入矩阵A/* 以下作A = A+I的运算*/(2)for i=0 to N-1 doa(i, i) = 1endfor/* 以下是A矩阵的㏒N次自乘，结果放入M矩阵；每次乘后，结果写回A矩阵*/(3)for k=1 to㏒N do(3.1)for i=0 to N-1 dofor j=0 to N-1 dos=0while (s<N) and (a(i,s)=0 or a(s,j)=0) dos = s+1endwhileif s<N then m(i,j)=1 else m(i,j)=0endforendfor/* 计算结果从M写回A */(3.2)for i=0 to N-1 dofor j=0 to N-1 doa(i, j) = m(i, j)endforendforendforEnd1.1.2 传递闭包并行算法本小节将把上一小节里的算法并行化。

在图论问题的并行化求解中，经常使用将N个顶点（或连通分量）平均分配给p个处理器并行处理的基本思想。

其中每个处理器分配到n 个顶点，即n=N/p。

无法整除时，一般的策略是在尽量均分的前提下，给最后一个处理器分配较少的顶点。

为了使算法简明，在本章中，仅在算法的MPI实现中才考虑不能整除的情况。

在以后的几节中，N、p、n都具有上面约定的意义，不再多说。

为了并行处理，这里将矩阵(A+I)进行划分，分别交给p个处理器。

在每次矩阵乘法的计算中，将N×N的结果矩阵(A+I)2均匀划分成p×p个子块，每块大小为n×n。

处理器i 负责计算位于第i子块行上的p个子块。

对整个子块行的计算由p次循环完成，每次计算一个子块。

每个处理器为了计算一个n×n大小的子块，需要用到源矩阵(A+I)中对应的连续n 行（局部数据a）和连续n列的数据（局部数据b）。

计算完成后，各处理器循环交换局部数据b，就可以进行下一个子块的计算了。

于是，总体算法由2层循环构成，外层是矩阵M=A+I的㏒N次自乘，每次首先完成矩阵M的一次自乘，然后将结果写回M；内层是p个子块的计算，每次首先完成各自独立的计算，然后处理器间循环交换局部数据b。

算法运行期间，矩阵M的数据作为全局数据由处理器0维护。

根据以上思想，并行算法描述见下面的算法15.2，使用了p个处理器，其时间复杂度为O(N3/p㏒N)。

其中myid是处理器编号，下同。

算法15.2传递闭包并行算法输入：图G的布尔邻接矩阵A输出：A的传递闭包Mprocedure closure-parallelBegin/* 由处理器0读入矩阵A到M中，并进行M=M+I运算对应源程序中readmatrix()函数*/(1)if myid=0 then(1.1) 读入矩阵A，放入M中(1.2) for i=0 to N-1 dom(i,i)=1endforendif(2)各处理器变量初始化/* 以下是主循环，矩阵M的㏒N次自乘*/(3)for i=1 to㏒N par_do/* 以下向各处理器发送对应子块行和列数据对应源程序中sendmatrix()函数*/(3.1)for j=1 to p-1 do(ⅰ) 处理器0发送第j个子块行的数据给处理器j，成为j的局部数据a(ⅱ) 处理器0发送第j个子块列的数据给处理器j，成为j的局部数据bendfor/* 以下是各处理器接收接收和初始化局部数据a和b对应源程序中getmatrix()函数*/(3.2)处理器0更新自己的局部数据a和b(3.3)处理器1到p-1从处理器0接受数据，作为局部数据a和b/* 以下是乘法运算过程，对应源程序中paramul()函数*/(3.4)for j=0 to p-1 do(ⅰ) 处理器k计算结果矩阵的子块(k, ((k+j) mod p))(ⅱ) 各处理器将数据b发送给前一个处理器(ⅲ) 各处理器从后一个处理器接收数据bendfor/* 以下是各处理器将局部计算结果写回M数组对应源程序中writeback()函数*/(3.5)if myid=0 then(ⅰ) 计算结果直接写入矩阵M(ⅱ) 接受其它处理器发送来的计算结果并写入Melse发送计算结果的myid子块行数据给处理器0endifendforEndMPI源程序请参见所附光盘。

1.2 连通分量图G的一个连通分量(Connected Component)是G的一个最大连通子图，该子图中每对顶点间均有一条路径。

根据图G，如何找出其所有连通分量的问题称为连通分量问题。

解决该问题常用的方法有3种：①使用某种形式的图的搜索技术；②通过图的布尔邻接矩阵计算传递闭包；③顶点倒塌算法。

本节将介绍如何在并行环境下实现顶点倒塌算法。

1.2.1 顶点倒塌法算法原理描述顶点倒塌（Vertex Collapse）算法中，一开始图中的N个顶点看作N个孤立的超顶点(Super Vertex)，算法运行中，有边连通的超顶点相继合并，直到形成最后的整个连通分量。

每个顶点属于且仅属于一个超顶点，超顶点中标号最小者称为该超顶点的根。

该算法的流程由一系列循环组成。

每次循环分为三步：①发现每个顶点的最小标号邻接超顶点；②把每个超顶点的根连到最小标号邻接超顶点的根上；③所有在第2步连接在一起的超顶点倒塌合并成为一个较大的超顶点。

图G的顶点总数为N，因为超顶点的个数每次循环后至少减少一半，所以把每个连通分量倒塌成单个超顶点至多㏒N次循环即可。

顶点i所属的超顶点的根记为D(i)，则一开始时有D(i)=i，算法结束后，所有处于同一连通分量中的顶点具有相同的D(i)。

1.2.2 连通分量并行算法算法中为顶点设置数组变量D和C，其中D(i)为顶点i所在的超顶点号，C(i)为和顶点i或超顶点i相连的最小超顶点号等，根据程序运行的阶段不同，意义也有变化。

算法的主循环由5个步骤组成：①各处理器并行为每个顶点找出对应的C(i)；②各处理器并行为每个超顶点找出最小邻接超顶点，编号放入C(i)中；③修改所有D(i)=C(i)；④修改所有C(i)=C(C(i))，运行㏒N次；⑤修改所有D(i)为C(i)和D(C(i))中较小者。

其中最后3步对应意义为超顶点的合并。

顶点倒塌算法是专为并行程序设计的，多个顶点的处理具有很强的独立性，很适合分配给多个处理器并行处理。

这里让p个处理器分管N个顶点。

则顶点倒塌算法的具体描述见算法15.3，使用了p个处理器，其时间复杂度为O(N2/p㏒N)，其中步骤(2)为主循环，内含5个子步骤对应上面的描述。

算法15.3 顶点倒塌算法输入：图G的邻接矩阵A输出：向量D( 0 ：N-1 )，其中D(i)表示顶点i的标识procedure collapse-verticesBegin/* 初始化 */(1)for每个顶点i par_doD(i) = iendfor(2)for k=1 to㏒N do/* 以下并行为每个顶点找邻接超顶点中最小的对应源程序中D_to_C()函数 */(2.1) for每个i, j : 0 ≤ i, j ≤ N-1 par_doC(i) = min{ D(j) | a(i,j)=1 and D(i)≠D(j) }if没有满足要求的D(j) thenC(i) = D(i)endifendfor/* 以下并行求每个超顶点的最小邻接超顶点对应源程序中C_to_C()函数 */(2.2) for每个i, j : 0 ≤ i, j ≤ N-1 par_doC(i) = min{ C(j) | D(j)=i and C(j)≠i }if没有满足要求的C(j) thenC(i) = D(i)endifendfor(2.3)for每个i : 0 ≤ i ≤ N-1 par_doD(i) = C(i)endfor(2.4)for i=1 to㏒N do/* 以下对应源程序中CC_to_C()函数 */for每个j : 0 ≤ j ≤ N-1 par_doC(j) = C(C(j))endforendfor/* 以下对应源程序中CD_to_D()函数 */(2.5)for每个i : 0 ≤ i ≤ N-1 par_doD(i) = min{ C(i), D(C(i)) }endforendforEndMPI源程序请参见章末附录。

1.3 单源最短路径单源最短路径(Single Source Shortest Path)问题是指求从一个指定顶点s到其它所有顶点i之间的距离，因为是单一顶点到其它顶点的距离，所以称为单源。

设图G(V,E)是一个有向加权网络，其中V和E分别为顶点集合和边集合，其边权邻接矩阵为W，边上权值w(i,j) > 0，i,j∈V，V={0,1,…,N-1}。