《并行计算》课程辅导：图论

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外，其余顶点构成另一个连通分量