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,若∈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，若∈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
(2)for i-0 to N-1 do ifi≠ s then dist(i)=∞ (4)for=0 to N-l do (4.1)从ⅥL中找一个顶点u,使得dist(u)最小 (4.2)for(每个顶点veVL)and(∈E)d if dist(u)+w(u, v)<dist(v) then dist( v)=dist(u)+w(u, v) endif 5)VL=VL-ful endfor 132最短路径并行算法 在上一小节串行算法的基础上,这里将其并行化。思路如下 (1)上述算法的(1)(2)两步中,每个处理器分派n=N/p个节点进行初始化 (2)第(4.1)步可以并行化为:首先每个处理器分派n个节点分别求局部最小值,然后 再p个处理器合作求全局最小值,最后再将这一最小值广播出去。p个处理器合作方法如下: 当p为偶数时,后p/2个处理器将自己的局部最小值分别发送到对应的前p/2个处理器中, 由前p/2个处理器比较出2个局部最小值中相对较小者并保留。当p为奇数时,设p=2h+1, 则后h个处理器的值分别发送到前h个处理器中,比较并保留小值。一次这样的循环过后 p个最小值减少了一半,两次后,变为1/4,如此一层一层的比较,lgP次循环后,就可以 得出唯一的全局最小值 (3)第(4.2)步可以每个处理器分配n个顶点,然后独立进行更新dst的操作。 根据以上思路,最短路径的并行算法见算法15.5,使用了p个处理器,其时间复杂度 为O(N2/p+Ngp)。这里的实现算法和对应的源程序中,处理器0只进行输入输出工作,不 参与任何其它计算;因此实际参加运算的处理器为p1个,所以有n=Np-1);另外,布尔 数组bdst用来标记各顶点是否已从ⅥL中取出 算法155 Dijkstra并行算法 输入:边权邻接矩阵W(约定顶点i,j之间无边连接时w(ij)=∞,且w(i)=0)、 待计算顶点的标号s 输出:dis(0:N-1),其中ds(表示顶点s到顶点i的最短路径(1≤i≤N edure Dijkstra-parallel Begin /*以下对应源程序中 ReadMatrixO函数* (1)处理器0读入边权邻接矩阵W /*以下初始化dst和 bdist数组,对应源程序中IntO函数* 2)for每个顶点 i par do if i=s then bdist(i= true dist(i)=w(i, s)
(2) for i=0 to N-1 do if i≠s then dist(i) = ∞ endfor (3) VL = V (4) for i=0 to N-1 do (4.1) 从 VL 中找一个顶点 u，使得 dist(u)最小 (4.2) for (每个顶点 v∈VL) and (∈E) do if dist(u)+w(u,v)<dist(v) then dist(v) = dist(u)+w(u,v) endif endfor (5) VL = VL-{u} endfor End 1.3.2 最短路径并行算法 在上一小节串行算法的基础上，这里将其并行化。思路如下： ⑴ 上述算法的(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()函数 */ (1) 处理器 0 读入边权邻接矩阵 W /* 以下初始化 dist 和 bdist 数组，对应源程序中 Init()函数 */ (2) for 每个顶点 i par_do if i=s then dist(i) = 0 bdist(i) = TRUE else dist(i) = w(i,s)
bdist(i= false endif endor /*以下是算法主循环,对应源程序中 FindMin Wayo函数* (3)for i=l to N-l do /*各处理器找出自己负责范围内未搜索节点中最小的dst* (31)for每个顶点 j par_ do index=min(j bdistO=FALSE) num=dist(index) endfor *以下各处理器协作对p-1个 index求最小* caninum=p for j=l to log(p-1)p if calum mod2=0then∥本次循环参加比较的数据个数为偶数 (i)calum =calum/2 (ii)序号大于caum的处理器将 index和num发送给序号比自己小 calum的处理器 (i)l序号小于 calum的处理器(不包含0号)在自己原有的和收到 的两个num之间,比较并保留较小者,同时记录对应的 index 本次循环参加比较的数据个数为奇数 (i)calum=(calum+1)/2 (i)序号大于 canon的处理器将 index和num发送给序号比自己小 calum的处理器 (i)序号小于 cannon的处理器(不包含0号)在自己原有的和收到 的两个num之间,比较并保留较小者,同时记录对应的ind endif endor (3.3)处理器1广播通知其它处理器自己的num和 index /*以下并行更新* 34)for每个顶点 j par_ do if bdist(=FALSE then dist()=min( dist(), num+wG, index)) endif endfor (3.5)顶点 index对应处理器将对应的 bdist( index)更改为TRUE MPI源程序请参见所附光盘。 14最小生成树 一个无向连通图G的生成树是指满足如下条件的G的子图T:①T和G顶点数相同 ②T有足够的边使得所有顶点连通,同时不出现回路。如果对G的每条边指定一个权值, 那么,边权总和最小的生成树称为最小代价生成树,记为MCST( Minimum Cost Spanning
bdist(i) = FALSE endif endfor /* 以下是算法主循环，对应源程序中 FindMinWay()函数 */ (3) 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 源程序请参见所附光盘。 1.4 最小生成树 一个无向连通图 G 的生成树是指满足如下条件的 G 的子图 T：①T 和 G 顶点数相同； ②T 有足够的边使得所有顶点连通，同时不出现回路。如果对 G 的每条边指定一个权值， 那么，边权总和最小的生成树称为最小代价生成树，记为 MCST（Minimum Cost Spanning
Tree),常简称为最小生成树(记为MST)。最小生成树问题就是给定G,找出G的一个最 小生成树T的问题。 本节将介绍用于求解MST问题的 Sollin算法,并将其实现。 141最小生成树串行算法 Sollin算法 (Sollin Algorithm)是众多用于求解MST问题的算法之一,其原理为:算法开 始时,图中的N个顶点视为一片森林,每个顶点视为一棵树:算法主循环的流程中,森林 里每棵树同时决定其连向其它树的最小邻边,并将这些边加入森林中,实现树的合并:循环 到森林中只剩一棵树时终止。由于森林中树的数目每次至少减少一半,所以只要lgN次循 环就可以找出MST。 算法中以D(i)为顶点i所在的树的编号,edge(i)和 closet(i)为树i连向其它树的最小邻 边和其长度,则 Sollin算法的形式描述见算法15.6,其时间复杂度为O(N2begN)。 算法156 Sollin mst算法 输入:无向图G的边权矩阵W 输出:G的最小生成树的边集合T procedure Sollin-MST /*以下所有顶点初始化为一棵孤立的树* (1)for i=0 to N-l do (1) endfor (2)T初始化为空集 3)while T<N-1 do /*以下为各树寻找连向其它树的最小边权* (3.1)for每棵树ido closet()=∝ endfor (32)for图G中每条边(v,u)do if D(v)*D(u) and w(v, u)<closet(D(v))then (i)closet(D(v))=w(v,u) endif endfor /*以下是树的合并* (33)for每棵树ido (i)ifD(v)≠Du)then T=T+{(v,u)} 树D(v)和D(u合并 endor endwhile
Tree），常简称为最小生成树（记为 MST）。最小生成树问题就是给定 G，找出 G 的一个最 小生成树 T 的问题。 本节将介绍用于求解 MST 问题的 Sollin 算法，并将其实现。 1.4.1 最小生成树串行算法 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 /* 以下所有顶点初始化为一棵孤立的树 */ (1) for i=0 to N-1 do D(i) = i endfor (2) T 初始化为空集 (3) 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
142最小生成树并行算法 在上一小节 Sollin算法的基础上,本节将其并行化。基本思路是由多个处理器同时负责 为多个树查找最短边。为了方便并行处理,各树寻找外连最短边的的任务分两步进行:①所 有处理器独立并行为每个顶点找连向其它树的最短边,并将这些数据保存:②各处理器独立 并行检索每棵树的各顶点,为整棵树找出外连最短边。注意在以上两步中,处理器间分别按 照顶点和树进行分配,对象有了变化 为了配合算法运行,设置了4个一维数组D、C、J、W:①Di)为顶点i所在的树的编 号,也就是该树中最初的顶点号,初始化为i②C(0为离顶点i最近的其它树的编号:③J(i) 为对应的顶点号:④W(i)为对应的边权,即距离。以上变量约定对MST并行算法的源程序 同样适用。运行期间每个处理器处理n=N/p个顶点或n个连通分量。 各树找出外连最短边后,接下来分3步进行树的合并:①让所有的DiC(i);②对所 有i,进行C(i=C(C(i)运算,并循环吧N次:③对所有i,更新D(为C(和DC()中较小 者。以上合并过程类似顶点倒塌算法中超顶点的合并 设MST为输出的最小生成树的邻接矩阵,运行期间由0号处理器维护,则并行算法的 描述见算法157,使用了p个处理器,其时间复杂度为O(N/plgN) 算法157并行 Sollin MSt算法 输入:原始图G的邻接矩阵A 输出:所求最小生成树的邻接矩阵MST procedure Sollin-MST-parallel /*以下初始化将图初始化为N棵树* (1)for每个顶点 i par do endfor (2)whle图G未连通do /*以下各处理器为所负责的顶点找出距离最近的树 对应源程序中 d to CO函数* (2.1for每个顶点ipar_do (i)W(=MAX (i)for每个顶点jdo if a(1,j>0 and DO*D(i and a(1,j<w(i) then C(1=DO) W()=a(1j) J()=J endif (ni)if w(i)=MAX then C()=D(1) endif endfor /*以下各处理器为所负责的树找出外连的最短边 对应源程序中 C to CO函数/ (2.2)for每棵树ipar_do
1.4.2 最小生成树并行算法 在上一小节 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 棵树 */ (1) for 每个顶点 i par_do D(i) = i endfor (2) 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
MAX (i)for每个顶点jdo j)= i and C()≠ i and w() th C()=C() temp=wO) dif (iv) if temp<N and J(temp)<n then 通知处理器0将边 temp, J(temp)加入MST数组 (2. 3)for i=0 to N-1 (2. 4)for i=l to log N /*以下对应源程序中 CC to CO函数* for j=0 to N-l par do C()=C(C() /*以下对应源程序中 Cd to DO函数*/ (2. 5)for i=0 to N-l par do D()=min{C(1),D(C(U)} endwhile MPI源程序请参见所附光盘。 1.5小结 本章针对传递闭包、连通分量、单源最短路径、最小生成树等4个经典图论问题,分 别介绍了它们的一种典型算法,以及在MPI并行环境下的算法实现。其中除连通分量外, 其它3个问题的并行算法实现均由串行算法并行化而得到;通过这3个问题,较为详细的介 绍了图论问题串行解法并行化的一般思路:均匀数据划分,独立计算。如何划分才能尽量做 到计算的独立,是划分的关键。在最小生成树 Sollin算法的并行化中,甚至在程序运行中按 照不同的对象对数据进行了划分 除了常见的串行算法并行化之外,本章还通过连通分量问题的解法,接触了根据并行机 特点,直接为并行环境设计的图论问题解法。这种情况在非数值算法问题中并非特例 本章的编写主要参考了文献。其中,并行传递闭包和连通分量算法也可以参见文献21, 单源最短路径 Dijkstra算法也可以参见文献,而最小生成树Soli算法也可以参见文献凵
(ⅱ) 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 源程序请参见所附光盘。 1.5 小结 本章针对传递闭包、连通分量、单源最短路径、最小生成树等 4 个经典图论问题，分 别介绍了它们的一种典型算法，以及在 MPI 并行环境下的算法实现。其中除连通分量外， 其它 3 个问题的并行算法实现均由串行算法并行化而得到；通过这 3 个问题，较为详细的介 绍了图论问题串行解法并行化的一般思路：均匀数据划分，独立计算。如何划分才能尽量做 到计算的独立，是划分的关键。在最小生成树 Sollin 算法的并行化中，甚至在程序运行中按 照不同的对象对数据进行了划分。 除了常见的串行算法并行化之外，本章还通过连通分量问题的解法，接触了根据并行机 特点，直接为并行环境设计的图论问题解法。这种情况在非数值算法问题中并非特例。 本章的编写主要参考了文献[1]。其中，并行传递闭包和连通分量算法也可以参见文献[2]， 单源最短路径 Dijkstra 算法也可以参见文献[3]，而最小生成树 Sollin 算法也可以参见文献[4]