3 图论
图论在计算机科学、信息科学、人工智能、网络理论、系统工程、控制论、运筹学和经济管理等领域有着广泛的应用。但很多图论问题虽易表达,却难以求解,其中有相当多的图论问题均属NP完全问题。本章主要介绍工程实用简单图论问题的并行算法及其MPI编程实现,包括传递闭包、连通分量、最短路径和最小生成树等。
传递闭包
设A是一个含N个顶点的有向图G的布尔邻接矩阵(Boolean Adjacent Matrix),即元素aij=1当且仅当从顶点i到j有一条有向边。所谓A的传递闭包(Transitive Closure),记为A+,是一个N×N的布尔矩阵,其元素bij=1当且仅当:①i=j;或②从i出发存在有向路径到j,又称为顶点i到j可达。从G的布尔邻接矩阵A求出G的传递闭包,就称为传递闭包问题。
传递闭包问题有很强的应用背景,在科学计算中广泛存在。传递闭包问题的经典解法之一就是利用布尔矩阵的乘法来求解。本节将把这一算法分别在串行和并行机上实现。
传递闭包串行算法
利用布尔矩阵相乘求解传递闭包问题的原理为:对于布尔矩阵(A+I)k中的任一元素bij,bij=1表示从i到j存在长度小于等于k的可达路径,否则,bij=0。显然对于k=1,(A+I)1中bij=1当且仅当从i到j路径长度为0(i=j)或为1(从i到j存在有向边);(A+I)2中,bij=1当且仅当从i到j路径长度小于等于2;((A+I)2) 2中,bij=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的传递闭包M
procedure closure
Begin
读入矩阵A
/* 以下作A = A+I的运算 */
for i=0 to N-1 do
a(i, i) = 1
endfor
/* 以下是A矩阵的㏒N次自乘,结果放入M矩阵;每次乘后,结果写回A矩阵 */
for k=1 to ㏒N do
(3.1)for i=0 to N-1 do
for j=0 to N-1 do
s=0
while (s<N) and (a(i,s)=0 or a(s,j)=0) do
s = s+1
endwhile
if s<N then m(i,j)=1 else m(i,j)=0
endfor
endfor
/* 计算结果从M写回A */
(3.2)for i=0 to N-1 do
for j=0 to N-1 do
a(i, j) = m(i, j)
endfor
endfor
endfor
End
传递闭包并行算法
本小节将把上一小节里的算法并行化。在图论问题的并行化求解中,经常使用将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的传递闭包M
procedure closure-parallel
Begin
/* 由处理器0读入矩阵A到M中,并进行M=M+I运算
对应源程序中readmatrix()函数 */
if myid=0 then
(1.1) 读入矩阵A,放入M中
(1.2) for i=0 to N-1 do
m(i,i)=1
endfor
endif
各处理器变量初始化
/* 以下是主循环,矩阵M的㏒N次自乘 */
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的局部数据b
endfor
/* 以下是各处理器接收接收和初始化局部数据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发送给前一个处理器
(ⅲ) 各处理器从后一个处理器接收数据b
endfor
/* 以下是各处理器将局部计算结果写回M数组
对应源程序中writeback()函数 */
(3.5)if myid=0 then
(ⅰ) 计算结果直接写入矩阵M
(ⅱ) 接受其它处理器发送来的计算结果并写入M
else
发送计算结果的myid子块行数据给处理器0
endif
endfor
End
MPI源程序请参见所附光盘。
连通分量
图G的一个连通分量(Connected Component)是G的一个最大连通子图,该子图中每对顶点间均有一条路径。根据图G,如何找出其所有连通分量的问题称为连通分量问题。解决该问题常用的方法有3种:①使用某种形式的图的搜索技术;②通过图的布尔邻接矩阵计算传递闭包;③顶点倒塌算法。本节将介绍如何在并行环境下实现顶点倒塌算法。
顶点倒塌法算法原理描述
顶点倒塌(Vertex Collapse)算法中,一开始图中的N个顶点看作N个孤立的超顶点(Super Vertex),算法运行中,有边连通的超顶点相继合并,直到形成最后的整个连通分量。每个顶点属于且仅属于一个超顶点,超顶点中标号最小者称为该超顶点的根。
该算法的流程由一系列循环组成。每次循环分为三步:①发现每个顶点的最小标号邻接超顶点;②把每个超顶点的根连到最小标号邻接超顶点的根上;③所有在第2步连接在一起的超顶点倒塌合并成为一个较大的超顶点。
图G的顶点总数为N,因为超顶点的个数每次循环后至少减少一半,所以把每个连通分量倒塌成单个超顶点至多㏒N次循环即可。顶点i所属的超顶点的根记为D(i),则一开始时有D(i)=i,算法结束后,所有处于同一连通分量中的顶点具有相同的D(i)。
连通分量并行算法
算法中为顶点设置数组变量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-vertices
Begin
/* 初始化 */
for 每个顶点i par_do
D(i) = i
endfor
for k=1 to ㏒N do
/* 以下并行为每个顶点找邻接超顶点中最小的
对应源程序中D_to_C()函数 */
(2.1) for 每个i, j : 0 ≤ i, j ≤ N-1 par_do
C(i) = min{ D(j) | a(i,j)=1 and D(i)≠D(j) }
if 没有满足要求的D(j) then
C(i) = D(i)
endif
endfor
/* 以下并行求每个超顶点的最小邻接超顶点
对应源程序中C_to_C()函数 */
(2.2) for 每个i, j : 0 ≤ i, j ≤ N-1 par_do
C(i) = min{ C(j) | D(j)=i and C(j)≠i }
if 没有满足要求的C(j) then
C(i) = D(i)
endif
endfor
(2.3) for 每个i : 0 ≤ i ≤ N-1 par_do
D(i) = C(i)
endfor
(2.4) for i=1 to ㏒N do
/* 以下对应源程序中CC_to_C()函数 */
for 每个j : 0 ≤ j ≤ N-1 par_do
C(j) = C(C(j))
endfor
endfor
/* 以下对应源程序中CD_to_D()函数 */
(2.5) for 每个i : 0 ≤ i ≤ N-1 par_do
D(i) = min{ C(i), D(C(i)) }
endfor
endfor
End
MPI源程序请参见章末附录。
单源最短路径
单源最短路径(Single Source Shortest Path)问题是指求从一个指定顶点s到其它所有顶点i之间的距离,因为是单一顶点到其它顶点的距离,所以称为单源。设图G(V,E)是一个有向加权网络,其中V和E分别为顶点集合和边集合,其边权邻接矩阵为W,边上权值w(i,j) > 0,i,j∈V,V={0,1,…,N-1}。
设dist( i )为最短的路径长度,即距离,其中s∈V且i≠s。这里采用著名的Dijkstra算法,并将其并行化。
最短路径串行算法
Dijkstra算法(Dijkstra Algorithm)是单源最短路径问题的经典解法之一,基本思想如下。
假定有一个待搜索的顶点表VL,初始化时做: dist(s)←0,dist(i)=∞(i≠s),VL=V。每次从VL(非空集)中选取这样的一个顶点u,它的dist(u)最小。将选出的u点作为搜索顶点,对于其它还在VL内的顶点v,若<u,v>∈E,而且dist(u)+w(u,v)<dist(v),则更新dist(v)为dist(u)+w(u,v),同时从VL中将u删除,直到VL 成为空集时算法终止。
算法15.4中给出了最短路径问题的Dijkstra串行算法,其时间复杂度为O(N2)。
算法 15.4Dijkstra串行算法
输入:边权邻接矩阵W(约定顶点i,j之间无边连接时w(i,j)=∞,且w(i,i) = 0)、
待计算顶点的标号s
输出:dist(0 : N-1),其中dist(i)表示顶点s到顶点i的最短路径(1≤i≤N)
procedure Dijkstra
Begin
dist(s) = 0
for i=0 to N-1 do
if i≠s then dist(i) = ∞
endfor
VL = V
for i=0 to N-1 do
(4.1) 从VL中找一个顶点u,使得dist(u)最小
(4.2) for (每个顶点v∈VL) and (<u,v>∈E) do
if dist(u)+w(u,v)<dist(v) then dist(v) = dist(u)+w(u,v)
endif
endfor
VL = VL-{u}
endfor
End
最短路径并行算法
在上一小节串行算法的基础上,这里将其并行化。思路如下:
⑴ 上述算法的(1)(2)两步中,每个处理器分派n=N/p个节点进行初始化。
⑵ 第(4.1)步可以并行化为:首先每个处理器分派 n个节点分别求局部最小值,然后再p个处理器合作求全局最小值,最后再将这一最小值广播出去。p个处理器合作方法如下:当p为偶数时,后p/2个处理器将自己的局部最小值分别发送到对应的前p/2个处理器中,由前p/2个处理器比较出2个局部最小值中相对较小者并保留。当p为奇数时,设p=2h+1,则后h个处理器的值分别发送到前h个处理器中,比较并保留小值。一次这样的循环过后,p个最小值减少了一半,两次后,变为1/4,如此一层一层的比较,㏒P次循环后,就可以得出唯一的全局最小值。
⑶ 第(4.2)步可以每个处理器分配n个顶点,然后独立进行更新dist的操作。
根据以上思路,最短路径的并行算法见算法15.5,使用了p个处理器,其时间复杂度为O(N2/p+N㏒p)。这里的实现算法和对应的源程序中,处理器0只进行输入输出工作,不参与任何其它计算;因此实际参加运算的处理器为p-1个,所以有n=N/(p-1);另外,布尔数组bdist用来标记各顶点是否已从VL中取出。
算法 15.5 Dijkstra并行算法
输入:边权邻接矩阵W(约定顶点i,j之间无边连接时w(i,j)=∞,且w(i,i) = 0)、
待计算顶点的标号s
输出:dist(0 : N-1),其中dist(i)表示顶点s到顶点i的最短路径(1≤i≤N)
procedure Dijkstra-parallel
Begin
/* 以下对应源程序中ReadMatrix()函数 */
处理器0读入边权邻接矩阵W
/* 以下初始化dist和bdist数组,对应源程序中Init()函数 */
for 每个顶点i par_do
if i=s then
dist(i) = 0
bdist(i) = TRUE
else
dist(i) = w(i,s)
bdist(i) = FALSE
endif
endfor
/* 以下是算法主循环,对应源程序中FindMinWay()函数 */
for i=1 to N-1 do
/* 各处理器找出自己负责范围内未搜索节点中最小的dist */
(3.1) for 每个顶点j par_do
index = min{ j | bdist(j)=FALSE }
num = dist(index)
endfor
/* 以下各处理器协作对p-1个index求最小 */
(3.2) calnum = p-1
for j=1 to ㏒(p-1) par_do
if calnum mod 2=0 then //本次循环参加比较的数据个数为偶数
(ⅰ) calnum = calnum/2
(ⅱ) 序号大于calnum的处理器将index和num发送给序号比自己小
calnum的处理器
(ⅲ) 序号小于calnum的处理器(不包含0号)在自己原有的和收到
的两个num之间,比较并保留较小者,同时记录对应的index
else //本次循环参加比较的数据个数为奇数
(ⅰ) calnum = (calnum+1)/2
(ⅱ) 序号大于calnum的处理器将index和num发送给序号比自己小
calnum的处理器
(ⅲ) 序号小于calnum的处理器(不包含0号)在自己原有的和收到
的两个num之间,比较并保留较小者,同时记录对应的index
endif
endfor
(3.3) 处理器1广播通知其它处理器自己的num和index
/* 以下并行更新 */
(3.4) for 每个顶点j par_do
if bdist(j)=FALSE then dist(j) = min{ dist(j), num+w(j,index) }
endif
endfor
(3.5) 顶点index对应处理器将对应的bdist(index)更改为TRUE
endfor
End
MPI源程序请参见所附光盘。
最小生成树
一个无向连通图G的生成树是指满足如下条件的G的子图T:①T和G顶点数相同;②T有足够的边使得所有顶点连通,同时不出现回路。如果对G的每条边指定一个权值,那么,边权总和最小的生成树称为最小代价生成树,记为MCST(Minimum Cost Spanning Tree),常简称为最小生成树(记为MST)。最小生成树问题就是给定G,找出G的一个最小生成树T的问题。
本节将介绍用于求解MST问题的Sollin算法,并将其实现。
最小生成树串行算法
Sollin算法(Sollin Algorithm)是众多用于求解MST问题的算法之一,其原理为:算法开始时,图中的N个顶点视为一片森林,每个顶点视为一棵树;算法主循环的流程中,森林里每棵树同时决定其连向其它树的最小邻边,并将这些边加入森林中,实现树的合并;循环到森林中只剩一棵树时终止。由于森林中树的数目每次至少减少一半,所以只要㏒N次循环就可以找出MST。
算法中以D(i)为顶点i所在的树的编号, edge(i) 和closet(i)为树i连向其它树的最小邻边和其长度,则Sollin算法的形式描述见算法15.6,其时间复杂度为O(N2㏒N)。
算法 15.6 Sollin MST算法
输入:无向图G的边权矩阵W
输出:G的最小生成树的边集合T
procedure Sollin-MST
Begin
/* 以下所有顶点初始化为一棵孤立的树 */
for i=0 to N-1 do
D(i) = i
endfor
T初始化为空集
while |T| < N-1 do
/* 以下为各树寻找连向其它树的最小边权 */
(3.1) for 每棵树i do
closet(i) = ∞
endfor
(3.2) for 图G中每条边(v,u) do
if D(v)≠D(u) and w(v,u)<closet(D(v)) then
(ⅰ) closet(D(v)) = w(v,u)
(ⅱ) edge(D(v)) = (v,u)
endif
endfor
/* 以下是树的合并 */
(3.3) for 每棵树i do
(ⅰ) (v,u) ← edge(i)
(ⅱ) if D(v)≠D(u) then
T = T + { (v,u) }
树D(v)和D(u)合并
endif
endfor
endwhile
End
最小生成树并行算法
在上一小节Sollin算法的基础上,本节将其并行化。基本思路是由多个处理器同时负责为多个树查找最短边。为了方便并行处理,各树寻找外连最短边的的任务分两步进行:①所有处理器独立并行为每个顶点找连向其它树的最短边,并将这些数据保存;②各处理器独立并行检索每棵树的各顶点,为整棵树找出外连最短边。注意在以上两步中,处理器间分别按照顶点和树进行分配,对象有了变化。
为了配合算法运行,设置了4个一维数组D、C、J、W:① D(i)为顶点i所在的树的编号,也就是该树中最初的顶点号,初始化为i;② C(i)为离顶点i最近的其它树的编号;③ J(i)为对应的顶点号;④ W(i)为对应的边权,即距离。以上变量约定对MST并行算法的源程序同样适用。运行期间每个处理器处理n=N/p个顶点或n个连通分量。
各树找出外连最短边后,接下来分3步进行树的合并:①让所有的D(i)=C(i);②对所有i,进行C(i)=C(C(i))运算,并循环㏒N次;③对所有i,更新D(i)为C(i)和D(C(i))中较小者。以上合并过程类似顶点倒塌算法中超顶点的合并。
设MST为输出的最小生成树的邻接矩阵,运行期间由0号处理器维护,则并行算法的描述见算法15.7,使用了p个处理器,其时间复杂度为O(N2/p㏒N)。
算法 15.7 并行Sollin MST算法
输入:原始图G的邻接矩阵A
输出:所求最小生成树的邻接矩阵MST
procedure Sollin-MST-parallel
Begin
/* 以下初始化将图初始化为N棵树 */
for 每个顶点i par_do
D(i) = i
endfor
while 图G未连通 do
/* 以下各处理器为所负责的顶点找出距离最近的树
对应源程序中D_to_C()函数 */
(2.1) for 每个顶点i par_do
(ⅰ) W(i) = MAX
(ⅱ) for 每个顶点j do
if a(i, j)>0 and D(j)≠D(i) and a(i, j)<W(i) then
C(i) = D(j)
W(i) = a(i,j)
J(i) = j
endif
endfor
(ⅲ) if W(i)=MAX then
C(i) = D(i)
endif
endfor
/* 以下各处理器为所负责的树找出外连的最短边
对应源程序中C_to_C()函数 */
(2.2) for 每棵树i par_do
(ⅰ) tempj = N+1
(ⅱ) tempw = MAX
(ⅲ) for 每个顶点j do
if D(j)=i and C(j)≠i and W(j)<tempw then
C(i) = C(j)
tempw = W(j)
tempj = j
endif
endfor
(ⅳ) if tempj<N and J(tempj)<N then
通知处理器0将边(tempj, J(tempj))加入MST数组
endif
endfor
(2.3) for i=0 to N-1 par_do
D(i) = C(i)
endfor
(2.4) for i=1 to ㏒N do
/* 以下对应源程序中CC_to_C()函数 */
for j=0 to N-1 par_do
C(j) = C(C(j))
endfor
endfor
/* 以下对应源程序中CD_to_D()函数 */
(2.5) for i=0 to N-1 par_do
D(i) = min{ C(i), D(C(i)) }
endfor
endwhile
End
MPI源程序请参见所附光盘。
小结
本章针对传递闭包、连通分量、单源最短路径、最小生成树等4个经典图论问题,分别介绍了它们的一种典型算法,以及在MPI并行环境下的算法实现。其中除连通分量外,其它3个问题的并行算法实现均由串行算法并行化而得到;通过这3个问题,较为详细的介绍了图论问题串行解法并行化的一般思路:均匀数据划分,独立计算。如何划分才能尽量做到计算的独立,是划分的关键。在最小生成树Sollin算法的并行化中,甚至在程序运行中按照不同的对象对数据进行了划分。
除了常见的串行算法并行化之外,本章还通过连通分量问题的解法,接触了根据并行机特点,直接为并行环境设计的图论问题解法。这种情况在非数值算法问题中并非特例。
本章的编写主要参考了文献[1]。其中,并行传递闭包和连通分量算法也可以参见文献[2],单源最短路径Dijkstra算法也可以参见文献[3],而最小生成树Sollin算法也可以参见文献[4]。
参考文献
陈国良 编著.并行算法的设计与分析(修订版).高等教育出版社,2002.11
Hirschberg D S. Parallel Algorithms for Transitive Closure and the Connected Component Problem. Proc. 8th Annu. ACM STOC, New York, 1976,55-57
Dijkstra E. A Note on Two Problems in Connetion with Graphs. Numerische Mathematic, 1959,1:269-271
Goodman S E and Hedetnimi S T ed. Introduction to the Design and Analysis of Algorithms(in An Algorithm Attributed to Sollin). McGraw-Hill, New York, 1977
附录 连通分量并行算法的MPI源程序
1. 源程序connect.c
/*本算法中处理器数目须小于图的顶点数*/
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <math.h>
#include <mpi.h>
/*使用动态分配的内存存储邻接矩阵,以下为宏定
义*/
#define A(i,j) A[i*N+j]
/*以下为
N:顶点数
n:各处理器分配的顶点数
p:处理器数*/
int N;
int n;
int p;
int *D,*C;
int *A;
int temp;
int myid;
MPI_Status status;
/*输出必要信息*/
void print(int *P)
{
int i;
if(myid==0)
{
for(i=0;i<N;i++)
printf("%d ",P[i]);
printf("\n");
}
}
/*读入邻接矩阵*/
void readA()
{
char *filename;
int i,j;
printf("\n");
printf("Input the vertex num:\n");
scanf("%d",&N);
n=N/p;
if(N%p!=0) n++;
A=(int*)malloc(sizeof(int)*(n*p)*N);
if(A==NULL){
printf("Error when allocating memory\n");
exit(0);
}
printf("Input the adjacent matrix:\n");
for(i=0;i<N;i++)
for(j=0;j<N;j++)
scanf("%d",&A(i,j));
for(i=N;i<n*p;i++)
for(j=0;j<N;j++)
A(i,j)=0;
}
/*处理器0广播特定数据*/
void bcast(int *P)
{
MPI_Bcast(P,N,MPI_INT,0,
MPI_COMM_WORLD);
}
/*两者中取最小的数学函数*/
int min(int a,int b)
{
return(a<b?a:b);
}
/*为各顶点找最小的邻接超顶点,对应算法步骤
(2.1)*/
void D_to_C()
{
int i,j;
for(i=0;i<n;i++){
C[n*myid+i]=N+1;
for(j=0;j<N;j++)
if((A(i,j)==1)
&&(D[j]!=D[n*myid+i])
&&(D[j]<C[n*myid+i])){
C[n*myid+i]=D[j];
}
if(C[n*myid+i]==N+1)
C[n*myid+i]=D[n*myid+i];
}
}
/*为各超顶点找最小邻接超顶点,对应算法步骤(2.2)*/
void C_to_C()
{
int i,j;
for(i=0;i<n;i++){
temp=N+1;
for(j=0;j<N;j++)
if((D[j]==n*myid+i)
&&(C[j]!=n*myid+i)
&&(C[j]<temp)){
temp=C[j];
}
if(temp==N+1) temp=D[n*myid+i];
C[myid*n+i]=temp;
}
}
/*调整超顶点标识*/
void CC_to_C()
{
int i;
for(i=0;i<n;i++)
C[myid*n+i]=C[C[myid*n+i]];
}
/*产生新的超顶点,对应算法步骤(2.5)*/
void CD_to_D()
{
int i;
for(i=0;i<n;i++)
D[myid*n+i]=
min(C[myid*n+i],D[C[myid*n+i]]);
}
/*释放动态内存*/
void freeall()
{
free(A);
free(D);
free(C);
}
/*主函数*/
void main(int argc,char *argv)
{
int i,j,k;
double l;
int group_size;
/*以下变量用来记录运行时间*/
double starttime,endtime;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,
&group_size);
MPI_Comm_rank(MPI_COMM_WORLD,
&myid);
p=group_size;
MPI_Barrier(MPI_COMM_WORLD);
if(myid==0)
starttime=MPI_Wtime();
/*处理器0读邻接矩阵*/
if(myid==0) readA();
MPI_Barrier(MPI_COMM_WORLD);
MPI_Bcast(&N,1,MPI_INT,0,
MPI_COMM_WORLD);
if(myid!=0){
n=N/p;
if(N%p!=0) n++;
}
D=(int*)malloc(sizeof(int)*(n*p));
C=(int*)malloc(sizeof(int)*(n*p));
if(myid!=0)
A=(int*)malloc(sizeof(int)*n*N);
/*初始化数组D,步骤(1)*/
for(i=0;i<n;i++) D[myid*n+i]=myid*n+i;
MPI_Barrier(MPI_COMM_WORLD);
MPI_Gather(&D[myid*n],n,MPI_INT,D,n,
MPI_INT,0,MPI_COMM_WORLD);
bcast(D);
MPI_Barrier(MPI_COMM_WORLD);
/*处理器0向其它处理器发送必要数据*/
if(myid==0)
for(i=1;i<p;i++)
MPI_Send(&A(i*n,0),n*N,
MPI_INT,i,i,
MPI_COMM_WORLD);
else
MPI_Recv(A,n*N,MPI_INT,
0,myid,MPI_COMM_WORLD,
&status);
MPI_Barrier(MPI_COMM_WORLD);
l=log(N)/log(2);
/*主循环体:算法步骤(2)*/
for(i=0;i<l;i++){
if(myid==0) printf("Stage %d:\n",i+1);
/*算法步骤(2.1)*/
D_to_C();
MPI_Barrier(MPI_COMM_WORLD);
MPI_Gather(&C[n*myid],n,MPI_INT,C,
n,MPI_INT,0,
MPI_COMM_WORLD);
print(C);
bcast(C);
MPI_Barrier(MPI_COMM_WORLD);
/*算法步骤(2.2)*/
C_to_C();
print(C);
MPI_Barrier(MPI_COMM_WORLD);
MPI_Gather(&C[n*myid],n,MPI_INT,C,
n,MPI_INT,0,
MPI_COMM_WORLD);
MPI_Gather(&C[n*myid],n,MPI_INT,D,
n,MPI_INT,0,
MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
/*算法步骤(2.3)*/
if(myid==0)
for(j=0;j<n;j++)
D[j]=C[j];
/*算法步骤(2.4)*/
for(k=0;k<l;k++){
bcast(C);
CC_to_C();
MPI_Gather(&C[n*myid],n,
MPI_INT,C,n,MPI_INT,0,
MPI_COMM_WORLD);
}
bcast(C);
bcast(D);
/*算法步骤(2.5)*/
CD_to_D();
MPI_Gather(&D[n*myid],n,
MPI_INT,D,n,MPI_INT,0,
MPI_COMM_WORLD);
print(D);
bcast(D);
}
if(myid==0) printf("Result: \n");
print(D);
if(myid==0){
endtime=MPI_Wtime();
printf("The running time is %d\n",
endtime-starttime);
}
freeall();
MPI_Finalize();
}
2. 运行实例
编译:mpicc connect.c
运行:本实例中使用3个处理器。
mpirun –np 3 a.out
运行结果:
Input the vertex num:
8
Input the adjacent matrix:
0 0 0 0 0 0 0 1
0 0 0 0 0 1 1 0
0 0 0 0 0 0 0 0
0 0 0 0 0 1 0 1
0 0 0 0 0 0 1 1
0 1 0 1 0 0 0 0
0 1 0 0 1 0 0 0
1 0 0 1 1 0 0 0
最终输出结果:
Result:
0 0 2 0 0 0 0 0
说明:使用了和上一小节传递闭包时相同的输入矩阵,得到的结果也相同。除顶点2外,其余顶点构成另一个连通分量