3图论 图论在计算机科学、信息科学、人工智能、网络理论、系统工程、控制论、运筹学和经 济管理等领域有着广泛的应用。但很多图论问题虽易表达,却难以求解,其中有相当多的图 论问题均属NP完全问题。本章主要介绍工程实用简单图论问题的并行算法及其MPI编程 实现,包括传递闭包、连通分量、最短路径和最小生成树等。 1.1传递闭包 设A是一个含N个顶点的有向图G的布尔邻接矩阵( Boolean Adjacent Matrix),即元素 a1=1当且仅当从顶点i到j有一条有向边。所谓A的传递闭包( Transitive Closure),记为 A,是一个N×N的布尔矩阵,其元素b=1当且仅当:①j=j;或②从i出发存在有向路径 到j,又称为顶点i到j可达。从G的布尔邻接矩阵A求出G的传递闭包,就称为传递闭包 问题 传递闭包问题有很强的应用背景,在科学计算中广泛存在。传递闭包问题的经典解法之 就是利用布尔矩阵的乘法来求解。本节将把这一算法分别在串行和并行机上实现。 111传递闭包串行算法 利用布尔矩阵相乘求解传递闭包问题的原理为:对于布尔矩阵(A+)中的任一元素b b产=1表示从i到j存在长度小于等于k的可达路径,否则,b=0。显然对于k=1,(A+)中 b=1当且仅当从i到j路径长度为0(i=j)或为1(从i到j存在有向边);(A+D2中,b=1 当且仅当从i到j路径长度小于等于2;(A+)2)2中,b=1当且仅当从i到j路径长度小于等 于4,等等。因为任意两点间如果存在可达路径,长度最多为N-1,所以k≥N1时,(A+D)k 就是所求的传递闭包A+。于是(A+)矩阵的bN次自乘之后所得的矩阵就是所求的传递闭包 根据前面的叙述,很自然的有下面的传递闭包串行算法15.1,其时间复杂度为O(N3k 算法151传递闭包串行算法 输入:图G的布尔邻接矩阵A 输出:A的传递闭包M edure closure Begi (1)读入矩阵A /*以下作A=A+I的运算 (2)for 0 to N-1 do a(1,1)= endfor /*以下是A矩阵的lgN次自乘,结果放入M矩阵;每次乘后,结果写回A矩阵* (3)for k-l to log n do (3. 1)for i=0 to N-1 do for j=0 to N-l do
1. 3 图论 图论在计算机科学、信息科学、人工智能、网络理论、系统工程、控制论、运筹学和经 济管理等领域有着广泛的应用。但很多图论问题虽易表达，却难以求解，其中有相当多的图 论问题均属 NP 完全问题。本章主要介绍工程实用简单图论问题的并行算法及其 MPI 编程 实现，包括传递闭包、连通分量、最短路径和最小生成树等。 1.1 传递闭包 设 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 的传递闭包，就称为传递闭包 问题。 传递闭包问题有很强的应用背景，在科学计算中广泛存在。传递闭包问题的经典解法之 一就是利用布尔矩阵的乘法来求解。本节将把这一算法分别在串行和并行机上实现。 1.1.1 传递闭包串行算法 利用布尔矩阵相乘求解传递闭包问题的原理为：对于布尔矩阵(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 (1) 读入矩阵 A /* 以下作 A = A+I 的运算 */ (2) for i=0 to N-1 do a(i, i) = 1 endfor /* 以下是 A 矩阵的㏒ N 次自乘，结果放入 M 矩阵；每次乘后，结果写回 A 矩阵 */ (3) for k=1 to ㏒ N do (3.1)for i=0 to N-1 do for j=0 to N-1 do s=0
hile(s<N)and (a(i, s)=0 or a(s,j=0)do endwhile ifs<N then m(1,j=l else m(i,j=0 /*计算结果从M写回A* (3. 2)for i=0 to N for i=0 to N-l do enamor 112传递闭包并行算法 本小节将把上一小节里的算法并行化。在图论问题的并行化求解中,经常使用将N个 顶点(或连通分量)平均分配给p个处理器并行处理的基本思想。其中每个处理器分配到n 个顶点,即n=Np。无法整除时,一般的策略是在尽量均分的前提下,给最后一个处理器分 配较少的顶点。为了使算法简明,在本章中,仅在算法的MPI实现中才考虑不能整除的情 况。在以后的几节中,N、p、n都具有上面约定的意义,不再多说。 为了并行处理,这里将矩阵(A+1)进行划分,分别交给p个处理器。在每次矩阵乘法的 计算中,将NXN的结果矩阵(A+)2均匀划分成p×p个子块,每块大小为nxn。处理器 负责计算位于第i子块行上的p个子块。对整个子块行的计算由p次循环完成,每次计算 个子块。每个处理器为了计算一个n×n大小的子块,需要用到源矩阵(A+D中对应的连续n 行(局部数据a)和连续n列的数据(局部数据b)。计算完成后,各处理器循环交换局部数 据b,就可以进行下一个子块的计算了。 于是,总体算法由2层循环构成,外层是矩阵M=A+1的ogN次自乘,每次首先完成矩 阵M的一次自乘,然后将结果写回M;内层是p个子块的计算,每次首先完成各自独立的 计算,然后处理器间循环交换局部数据b。算法运行期间,矩阵M的数据作为全局数据由 处理器0维护 根据以上思想,并行算法描述见下面的算法15.2,使用了p个处理器,其时间复杂度为 O(NpgN)。其中myid是处理器编号,下同。 算法152传递闭包并行算法 输入:图G的布尔邻接矩阵A 输出:A的传递闭包M procedure closure-parallel 体*由处理器0读入矩阵A到M中,并进行M=M+运算 对应源程序中 readmatrixo函数*/ (1) if myid=0 then (1.1)读入矩阵A,放入M中 (1. 2)for i =0 to N-l do
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 1.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 的传递闭包 M procedure closure-parallel Begin /* 由处理器 0 读入矩阵 A 到 M 中，并进行 M=M+I 运算 对应源程序中 readmatrix()函数 */ (1) if myid=0 then (1.1) 读入矩阵 A，放入 M 中 (1.2) for i=0 to N-1 do m(i,i)=1
(2)各处理器变量初始化 *以下是主循环,矩阵M的ogN次自乘* (3)for i=l to log N par do /*以下向各处理器发送对应子块行和列数据 对应源程序中 sendmatrixo函数* (3. 1)for l to p-l de i)处理器0发送第j个子块行的数据给处理器j,成为j的局部数据a (ⅱi)处理器0发送第j个子块列的数据给处理器j,成为j的局部数据b endor /*以下是各处理器接收接收和初始化局部数据a和b 对应源程序中 getmatrixO函数* (32)处理器0更新自己的局部数据a和b (3.3)处理器1到p1从处理器0接受数据,作为局部数据a和b /*以下是乘法运算过程,对应源程序中 paramus函数* (3.4)for j=0 to p-l do (i)处理器k计算结果矩阵的子块(k,(k+j)modp) (ⅱi)各处理器将数据b发送给前一个处理器 i)各处理器从后一个处理器接收数据b /*以下是各处理器将局部计算结果写回M数组 对应源程序中 writeback(函数* (3.5 )if myid=0 then (i)计算结果直接写入矩阵M (ⅱi)接受其它处理器发送来的计算结果并写入M el se 发送计算结果的myid子块行数据给处理器0 endor MPI源程序请参见所附光盘。 12连通分量 图G的一个连通分量 Connected Component)是G的一个最大连通子图,该子图中每对 顶点间均有一条路径。根据图G,如何找出其所有连通分量的问题称为连通分量问题。解决 该问题常用的方法有3种:①使用某种形式的图的搜索技术:②通过图的布尔邻接矩阵计算 传递闭包:;③顶点倒塌算法。本节将介绍如何在并行环境下实现顶点倒塌算法。 121顶点倒塌法算法原理描述 顶点倒塌( Vertex Collapse)算法中,一开始图中的N个顶点看作N个孤立的超顶点(S Vertex),算法运行中,有边连通的超顶点相继合并,直到形成最后的整个连通分量。每个顶
endfor endif (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 的局部数据 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 源程序请参见所附光盘。 1.2 连通分量 图 G 的一个连通分量(Connected Component)是 G 的一个最大连通子图，该子图中每对 顶点间均有一条路径。根据图 G，如何找出其所有连通分量的问题称为连通分量问题。解决 该问题常用的方法有 3 种：①使用某种形式的图的搜索技术；②通过图的布尔邻接矩阵计算 传递闭包；③顶点倒塌算法。本节将介绍如何在并行环境下实现顶点倒塌算法。 1.2.1 顶点倒塌法算法原理描述 顶点倒塌（Vertex Collapse）算法中，一开始图中的N个顶点看作N个孤立的超顶点(Super Vertex)，算法运行中，有边连通的超顶点相继合并，直到形成最后的整个连通分量。每个顶
点属于且仅属于一个超顶点,超顶点中标号最小者称为该超顶点的根 该算法的流程由一系列循环组成。每次循环分为三步:①发现每个顶点的最小标号邻接 超顶点:②把每个超顶点的根连到最小标号邻接超顶点的根上:③所有在第2步连接在一起 的超顶点倒塌合并成为一个较大的超顶点 图G的顶点总数为N,因为超顶点的个数每次循环后至少减少一半,所以把每个连通 分量倒塌成单个超顶点至多lgN次循环即可。顶点i所属的超顶点的根记为D(i),则一开 始时有D(i)=i,算法结束后,所有处于同一连通分量中的顶点具有相同的D(i)。 122连通分量并行算法 算法中为顶点设置数组变量D和C,其中D(i)为顶点i所在的超顶点号,C(i)为和顶点 i或超顶点ⅰ相连的最小超顶点号等,根据程序运行的阶段不同,意义也有变化。算法的主 循环由5个步骤组成:①各处理器并行为每个顶点找出对应的C();②各处理器并行为每个 超顶点找出最小邻接超顶点,编号放入C(i)中;③修改所有D(i)=C(i);④修改所有 C(i)=C(C(i),运行lgN次:⑤修改所有D)为C(和D(C(i)中较小者。其中最后3步对应 意义为超顶点的合并。 顶点倒塌算法是专为并行程序设计的,多个顶点的处理具有很强的独立性,很适合分配 给多个处理器并行处理。这里让p个处理器分管N个顶点。则顶点倒塌算法的具体描述见 算法15.3,使用了p个处理器,其时间复杂度为ON/gN,其中步骤(2)为主循环,内含 5个子步骤对应上面的描述 算法153顶点倒塌算法 输入:图G的邻接矩阵A 输出:向量D(0:N-1),其中D(i)表示顶点i的标识 procedure collapse-vertices /*初始化*/ (1)for每个顶点 i par d D()=1 (2)for k=l to log n do /*以下并行为每个顶点找邻接超顶点中最小的 对应源程序中DtoC函数*/ (2.1)for每个i,j:0≤i,j≤N- I par do C(i)=min( DO I a(i,j=I and D(i)DO j if没有满足要求的Dj)then C()=D(1) endif endfor /*以下并行求每个超顶点的最小邻接超顶点 对应源程序中 C to CO函数* (2.2)for每个i,j:0≤i,j≤N- I par do C(i=min CO DO=i and CO*i) if没有满足要求的C()the C()=D()
点属于且仅属于一个超顶点，超顶点中标号最小者称为该超顶点的根。 该算法的流程由一系列循环组成。每次循环分为三步：①发现每个顶点的最小标号邻接 超顶点；②把每个超顶点的根连到最小标号邻接超顶点的根上；③所有在第 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-vertices Begin /* 初始化 */ (1) for 每个顶点 i par_do D(i) = i endfor (2) 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- I par do endfor (2. 4) for i=l to log n do /*以下对应源程序中 CC to CO函数* for每个j:0≤j≤N-1 par do C()=C(C() endfor /*以下对应源程序中 cd to d)函数* (2.5)for每个i:0≤i≤N- I par do D()=min{C(),D(C()} endion endfor MPⅠI源程序请参见章末附录 13单源最短路径 单源最短路径 Single source shortest path)问题是指求从一个指定顶点s到其它所有顶点 之间的距离,因为是单一顶点到其它顶点的距离,所以称为单源。设图G(VE)是一个有向 加权网络,其中V和E分别为顶点集合和边集合,其边权邻接矩阵为W,边上权值w(ij)> 设dis(i)为最短的路径长度,即距离,其中s∈Ⅴ且i≠s。这里采用著名的 Dijkstra算 法,并将其并行化 131最短路径串行算法 Dijkstra算法( Dijkstra Algorithm)是单源最短路径问题的经典解法之一,基本思想如 下 假定有一个待搜索的顶点表ⅥL,初始化时做:dst(s)←0,dist(i)=∞(i≠s),ⅥL=V 每次从VL(非空集)中选取这样的一个顶点u,它的dsu)最小。将选出的u点作为搜索 顶点,对于其它还在L内的顶点v,若<u>∈E而且dst(u)+w(uv)<dst(v),则更新dis(yv) 为disu)+wuv),同时从VL中将u删除,直到VL成为空集时算法终止 算法154中给出了最短路径问题的 Dijkstra串行算法,其时间复杂度为ON2) 算法154 Dijkstra串行算法 输入:边权邻接矩阵W(约定顶点i,j之间无边连接时w(ij)=∞,且w(i,i)=0) 待计算顶点的标号s 输出:dst(0:N-1),其中dst(表示顶点s到顶点i的最短路径(1≤i≤N procedure Dijkstra Begin 1)dist(s)=0
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 源程序请参见章末附录。 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}。 设 dist( i )为最短的路径长度，即距离，其中 s∈V 且 i≠s。这里采用著名的 Dijkstra 算 法，并将其并行化。 1.3.1 最短路径串行算法 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 (1) dist(s) = 0