单源最短路径问题并行算法分析

单源最短路径问题并行算法分析实验报告
一、实验名称
单源最短路径问题并行算法分析。

二、实验目的
分析单源最短路径Dijkstra并行算法和MPI源程序，并分析比较Dijkstra并行算法和Moore并行算法的性能。

三、实验内容
1、分析单源最短路径Dijkstra并行算法和MPI源程序。

2、分析单源最短路径问题的Moore并行算法，比较两种并行算法的性能。

四、实验步骤
1、问题描述
单源最短路径问题即指：已知一个n结点有向图G=(V,E)和边的权函数c(e)，求由G中某指定结点v0到其他各个结点的最短路径。

这里还假定所有的权值都是正的。

2、比较串行Dijkstra算法和Moore算法
2.1、Dijkstra算法基本思想
假定有一个待搜索顶点表VL，初始化时做：dist(s)←0；dist(i)←∞(i≠s)；VL←V。

算法执行时，每次从VL（≠Φ）中选取这样一个顶点u，它的dist(u)值最小。

将选出的u作为搜索顶点，若<u，v>∈E，而且dist(u)+w(u,v)<dist(v)，则更新dist(v)为dist(u)+w(u,v)，直到VL=Φ时算法终止。

算法描述如下：输入：加权邻接矩阵W，约定i，j之间无边连接时w(i,j)=∞，且w(i,i)=∞；
输出：dist(1:n)，其中，dist(i)表示顶点s到顶点i的最短路径（1≤i≤n）。

begin
/*初始化*/
(1)dist(s)←0；
(2)for i←1 to n do
if i≠s then dist(i)←∞endif
endfor；
(3)VL←V；
(4)for i←1 to n do /*找最短距离*/
(5)find a vertex u∈VL,such that dist(u) is minimal；
(6)for each(<u，v>∈E) ∧(v∈VL) do
if dist(u)+w(u,v)<dist(v) then
dist(v)←dist(u)+w(u,v)
endif
endfor；
(7)VL←VL-{u}
endfor
end.
2.2、Moore算法的基本思想
设源点为s∈V,从s到其它各顶点的最短路径长度用一个一
维数组dist存储。

首先置dist(s)=0，dist(v)←∞，v≠s，v∈V。

Moore算法同Dijkstra算法不同之处是：将Dijkstra算法中的待搜索顶点表替换成待搜索队列。

算法开始执行时，队列仅含源点s。

以后每次只要待搜索队列非空，则将队头的顶点从队列中移出作为本次搜索的顶点，然后检查u的所有射出边<u，v>∈E。

若dist(u)+w(u,v)<dist(v)，则此时找出了一条从s到v的更短路径，置dist(v)等于dist(u)+w(u,v)。

若v不在搜索队列中，则把v加入到队列的队尾。

如此重复进行，直到整个待搜索队列空时，算法终止。

单处理器上的Moore算法如下：
输入：有向图G(V,E)的加权邻接矩阵W={w ij},i,j∈V；
输出：从源s到所有其它顶点i(i≠s)的最短路径dist(i)，i∈V。

begin
(1)for i←1 to n do /*初始化*/
call INITIALIZE(i)
endfor;
(2) insert s into the queue; /*插入s到队列中*/
(3) while the queue is not empty do
call SEARCH
endwhile
end.
Procedure SEARCH;
begin
(3.1) dequeue vertex u; /*从队列中删去顶点u*/
(3.2) for each edge (u,v)∈E do
(ⅰ) new_dist←dist(u)+w(u,v);
(ⅱ) if new_dist<dist(v)
then dist(v)←new_dist
endif;
(ⅲ) if v is not in the queue then
enqueue vertex v
endif
endfor
end.
Procedure INITIALIZE
begin
(1.1)dist(s)←0
(1.2)dist(v)←∞(v≠s,v∈V)
(1.3)queue←0
end
2.3、两种算法的性能比较
Dijkstra算法是由E.W.Dijkstra于1959年提出的一个适用于所有弧度权均为非负的最短路径算法，也是目前公认的高效经典算法之一。

其时间复杂度为O (n2)，其中n为结点个数。

Moore算法分别由Bellman，Ford和Moore在五、六十年代提出。

其时间复杂度是O (nm)，期中m是边/弧数。

目前这样的时间复杂度在所有带有负权弧顶最短路径算法中是最好的。

但其实际运算效果却往往不及Dijkstra算法。

3、串行算法并行化
3.1、Dijkstra并行算法
1）每个处理器分配n=N/p各节点进行初始化（N为节点数，p为处理器个数）。

2）首先每个处理器分配n个节点分别求局部最小值，然后再p个处理器和作求全局最小值，最后再将这一最小值广播出去。

P 个处理器合作方法：当p为偶数时，后p/2 个处理器将自己的局部最小值分别发送到对应的前p/2个处理器中，由前p/2个处理器比较出2个局部最小值中相对较小者并保留。

当p为奇数时，设p=2h+1,则后h个处理器的值分别发送到前h个处理器中，比较并保留最小值。

一次这样的循环过程后，p个最小值减少了一半，两次后，变为1/4，如此一层一层的比较，logp次循环后，就可以得出唯一的全局最小值，然后将其广播出去。

3）每个处理器分配n个顶点，然后独立进行更新dist的操作。

其中广播实现步骤如下：
输入：数据u(存放在单元B(1)中)；
输出：将u广播到数组B的所有单元中去。

Begin
(1) B(1)←u；
log-1 do
(2) for i←0 to ⎡⎤p
(3) for each j:2i+1≤j≤2i+1 pardo
B(j)←B(j-2j)
endfor
endfor
End.
3.2、Moore并行算法
首先，队列用源点初始化。

然后创建了许多异步进程，每个进程都从队列中删除一个顶点，检查其射出边，将已发现有更短路径的顶点加入到队列。

算法的第（1）步采用预调（Prescheduling）方法很容易并行化。

而第（3）步的while循环需作适当的修改，为此引入两个用于同步的变量：数组变量waiting，它记住哪一个进程正处于等待状态；布尔变量halt，仅当队列为空和所有进程处于等待状态时为真。

INITIALIZE过程置数组waiting中的第一个元素为假。

SEARCH过程亦必须作适当的修改。

因为对队列的插入、删除操作不是原子操作，所以执行上述操作时必须给队列上锁；其次，在一个进程将刚找到的v路径new_dist与当前的最短路径dist(v)比较之前，变量dist(v)也必须上锁，否则两个进程有可能同时修改它；最后，若一个进程发现队列为空时，则置waiting中的相应元素为真。

若第一个进程处于等待状态，则它要检查每个进程是否都处在等待状态，如果是，则halt置为真，而在第一个进程检查每个进程是否都处于等待状态时，队列亦必须上锁。

算法描述：
MIMD紧耦合多处理机上的单源最短路径算法
输入：有向加权图G的权矩阵W；
输出：从源s到其余顶点i（i<>s）的最短路径dist(i), i属于V。

BEGIN
(1)for each i : 1 ≤i ≤p par-do /*p为进程数*/
for j = i to n step p do /*初始化*/
INITIALIZE(J)
End for
End for
(2)enqueue s /*s 入队*/
(3)halt ← false
(4)for each i : 1 ≤i ≤p par-do
while not halt do SEARCH(i) end while
end for
END
Procedure SEARCH(i)
Begin
(5)lock the queue /*队列上锁*/
(6)if queue is empty then /*队列空时，等待进程为真*/
(6.1) waiting (i) ← true
(6.2) if i = 1 then /*进程1等待时，其他进程均须等待*/
halt ← waiting(2)∧waiting(3)∧…∧waiting(p) end if
(6.3) unlock the queue /*队列开锁*/
else
(6.4)dequeue u /*从队列中删除u*/
(6.5)waiting(i)←false
(6.6)unlock the queue /*队列开锁*/
(6.7)for every edge(u,v)in graph do /*检查每条射出边*/
(ⅰ) new_dist←dist(u)+w(u,v);
(ⅱ) lock(dist(v)) /*dist(v)上锁*/
(ⅲ) if new_dist<dist(v) then
(a)dist(v)←new_dist /*更新new_dist*/
(b)unlock(dist(v)) /*dist(v)开锁*/
(c)if v !∈queue then
①lock the queue /*队列上锁*/
②enqueue v /*v入队*/
③unlock the queue /*队列开锁*/
endif
else unlock(dist(v)) /*dist(v)开锁*/
end if
end for
end if
end
值得指出的是：虽然创建更多的进程会缩短各个算法的执行时间（因为可以同时检查若干个顶点的射出边），但每个进程对队列的插入和删除要进行互斥控制，所以最大加速比最终还是受限制的。

3.3、两种并行算法的性能比较
Djkstra算法复杂性分析：根据以上步骤可得并行算法使用了p个处理器，其时间复杂度为O(N2/p+nlogp)。

这里应该注意的是处理器0只进行输入/输出工作，不参与任何其他的计算，因此实际参加运算的处理器为p-1个，所以有n=N/(p-1)。

在图搜索中，若使用Dijkstra算法则要求首先检查最近的顶点，因而Dijkstra算法在这个阶段中进行图的搜索时则必须按照顺序
处理的方式进行搜索。

这样Dijkstra算法在图的搜索阶段的并行性能也就大大降低了。

这无形中给多处理机系统的应用带来了极大的负面影响。

这时在Moore算法中，不必以特定的次序检查这些边，所以，多选用Moore算法。

4、并行算法实现
4.1、Dijkstra并行算法实现
以下是Dijkstra并行算法的核心代码：
辅助数据结构：
FILE *fp; /*存储图的邻接矩阵数据，第一个数值为矩阵的维数（即图中点的个数），剩余数据为按行存储的边权值，且如果两点之间无边相连，则权值由M或者m表示*/
Double *W; /*按行存储的图的邻接矩阵*/
int S = 0; /*源点s的序号*/
辅助函数原型说明：
4.1.1、v oid fatal(char* err)
函数功能：打印出错信息，并退出程序；
输入参数：出错信息字符串；
输出参数：无；
函数返回：无。

4.1.2、void GetChar(char &ch, FILE* fp)
函数功能：从文件fp中读出一个字节数据，存入变量ch；
输入参数：文件指针fp；
输出参数：字符引用&ch；
函数返回：无；
附加说明：文件中的合法字符为：数字、M、m，该函数要求过滤掉非法字符。

4.1.3、double GetNextNum()
函数功能：逐个读出文件fp中的数值；
输入参数：无；
输出参数：无；
函数返回：double型数值；
附加说明：该函数调用GetChar()函数，并用于读出邻接矩阵的维数以及边权值。

4.1.4、void ReadMatrix()
函数功能：从文件fp读入邻接矩阵W并分配存储空间，从终端读入源点的序号S；
输入参数：无；
输出参数：无；
函数返回：无。

核心函数说明及其实现：
4.1.5、int main(int argc,char** argv)
函数功能：算法主函数；
输入参数：int argc,char** argv；
输出参数：无；
函数返回：int型函数执行结果；
函数实现：
int main(int argc,char** argv)
{
int group_size = 0; /*所有参加运算的进程个数*/
int my_rank = 0; /*当前正在运行的进程的标识*/
MPI_Status status;
int i,j;
int ep;
int mynum;
MPI_Init(&argc,&argv); /*MPI begin*/
MPI_Comm_size(MPI_COMM_WORLD,&group_size);
MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
/*每个进程分配n=N/p个节点进行初始化（N为节点数，p为计算进程（不包括0号进程）个数）*/
/*0号进程读入邻接矩阵，并把顶点个数和源点标识广播给其他各个计算进程*/
if(my_rank == 0)
{
ReadMatrix();
for(i=1;i<group_size;i++)
{
MPI_Send(&nodenum,1,MPI_INT,i,i,
MPI_COMM_WORLD);
MPI_Send(&S,1,MPI_INT,i,i,
MPI_COMM_WORLD);
}
}
else
{/*非0号计算进程接收0号进程广播的顶点个数和源点标识*/
MPI_Recv(&nodenum,1,MPI_INT,0,my_rank,
MPI_COMM_WORLD,&status);
MPI_Recv(&S,1,MPI_INT,0,my_rank,
MPI_COMM_WORLD,&status);
}
/*0号进程不参与运算，所以计算的进程只group_size-1个*/
/*保证epnodenum/(group_size-1)向上取整*/
ep = nodenum/(group_size-1);
if(ep*group_size-ep < nodenum) ep++;
if(ep*my_rank <= nodenum) /*每个进程分到节点数ep*/
{
mynum = ep;
}
else if(ep*my_rank < nodenum+ep) /*最后一个进程分配的节点数可能小于平均值ep*/
{
mynum = nodenum - ep*(my_rank-1);
}
else mynum = 0;
if (my_rank == 0) mynum = 0; /*0号进程不参与计算*/
/*初始化各进程数据*/
Init(my_rank,group_size,ep);
/*输出图的邻接矩阵*/
OutPutMatrix(my_rank, group_size, ep, mynum);
/*单源最短路径主算法*/
FindMinWay(my_rank,group_size,ep,mynum);
OutPutResult(my_rank,group_size,ep,mynum);
MPI_Finalize();
free(W);
free(dist);
free(bdist);
return 0;
}
4.1.6、void Init(int my_rank,int group_size,int ep)
函数功能：初始化各个进程数据
输入参数：
int my_rank : 本进程在该通信域中的标识；
int group_size : 该通信域中并行的进程个数；
int ep : 每个计算进程处理的平均节点个数。

输出参数：无
函数返回：无
函数实现：
void Init(int my_rank,int group_size,int ep)
{
int i = 0;
MPI_Status status;
/*0号进程发送节点数据到其他各个计算进程，发送的数据为：按行划分的邻接矩阵数值，ep为该计算进程需要处理的节点个数，nodenum为邻接矩阵的维数，则发送给该计算进程的数值个数应为ep*nodenum，且该段数据的缓存起始地址为
&W(ep*(i-1),0)，该消息的标识即为接收进程的标识*/
if(my_rank == 0)
{
for(i=1;i<group_size;i++)
{
MPI_Send(&W(ep*(i-1),0),ep*nodenum,
MPI_DOUBLE,i,i, MPI_COMM_WORLD);
}
}
/*计算进程为得到的数据进行空间分配*/
else
{
/*从源节点到各个节点的最短路径*/
dist = (double*)malloc(sizeof(double)*ep);
/*是否找到了最短路径*/
bdist = (int*) malloc(sizeof(BOOL)*ep);
W = (double*)malloc(sizeof(double)*ep*nodenum);
if(W == NULL || dist == NULL || bdist == NULL)
{
fatal("Dynamic allocate space for matrix fail!");
}
/*从0号进程接收需要处理的节点数据*/
MPI_Recv(W,ep*nodenum,MPI_DOUBLE,0,my_rank,
MPI_COMM_WORLD,&status);
/*轮询分配到的节点：如果分配到的节点i是源节点S，则初始化dist[i]=0, bdist[i]=TRUE；否则，dist[i] = W(I,S);bdist[i] = FALSE；*/
for(i=0;i<ep; i++)
{
if(i+(my_rank-1)*ep == S)
{
dist[i] = 0;
bdist[i] = TRUE;
}
else
{
dist[i] = W(i,S);
bdist[i] = FALSE;
}
}
}
}
4.1.7、void OutPutMatrix(int my_rank,int group_size,int
ep,int mynum)
函数功能：输出邻接矩阵
输入参数：
int my_rank : 本进程在该通信域中的标识；
int group_size : 该通信域中并行的进程个数；
int mynum: 本进程需要处理的节点个数。

输出参数：无
函数返回：无
函数实现：
邻接矩阵的输出格式为：
Processor1: X X X X
Processor1: X X X X
void OutPutMatrix(int my_rank,int group_size,int ep,
int mynum)
{
int i,j;
if(my_rank != 0)
{
for(i=0;i<mynum;i++)
{
printf("Processor %d:\t",my_rank);
for(j=0;j<nodenum;j++)
{
if(W(i,j) > 1000000) printf("M\t");
else printf("%d\t",(int)W(i,j));
}
printf("\n");
}
}
}
4.1.8、void OutPutResult(int my_rank,int group_size,int
ep,int mynum)
函数功能：输出结果
输入参数：
int my_rank : 本进程在该通信域中的标识；
int group_size : 该通信域中并行的进程个数；
int mynum: 本进程需要处理的节点个数。

输出参数：无
函数返回：无
函数实现：
最短路径结果格式如下：
node i: X
node j: Y
void OutPutResult(int my_rank,int group_size,int ep,
int mynum)
{
int i,j;
if(my_rank != 0)
{
for(i=0;i<mynum;i++)
{
printf("node %d:\t%d\n",
(my_rank-1)*ep+i,(int)dist[i]);
}
}
}
4.1.9、void FindMinWay(int my_rank,int group_size,int
ep,int mynum)
函数功能：算法主循环
输入参数：
int my_rank : 本进程在该通信域中的标识；
int group_size : 该通信域中并行的进程个数；
int mynum: 本进程需要处理的节点个数。

输出参数：无
函数返回：无
函数实现：
void FindMinWay(int my_rank,int group_size,int ep,int mynum) {
int i,j;
int index,index2;
double num,num2;
int calnum;
MPI_Status status;
int p = group_size;
for(i=0; i<nodenum;i++)
{
index = 0;
num = M;
/*每个进程计算局部最短路径*/
for(j=0;j<mynum;j++)
{
if(dist[j] < num && bdist[j]==FALSE)
{
num = dist[j];
index = ep*(my_rank-1)+j; /*j在图中的序号*/ }
}
MPI_Barrier(MPI_COMM_WORLD); /*同步路障*/
calnum = group_size-1;
/*while循环对各个计算进程局部最小值进行规约*/
while(calnum > 1)
{
/**节点数目为偶数时**/
if(calnum % 2 == 0)
{
calnum = calnum/2;
/*后p/2 个进程将自己的局部最小值发到对应的前p/2个进程*/
if(my_rank > calnum)
{
MPI_Send(&index,1,MPI_INT,
my_rank-calnum,
my_rank-calnum,
MPI_COMM_WORLD);
MPI_Send(&num,1,MPI_DOUBLE,
my_rank-calnum,
my_rank-calnum,
MPI_COMM_WORLD);
}
/*前p/2 个进程接收后p/2个进程发过来的局部最小值，并比较得到更优的最小值*/
else if(my_rank!=0)
{
MPI_Recv(&index2,1,MPI_INT,
my_rank+calnum,
my_rank,
MPI_COMM_WORLD,&status);
MPI_Recv(&num2,1,MPI_DOUBLE,
my_rank+calnum,
my_rank,
MPI_COMM_WORLD,&status);
if(num2 < num)
{
num = num2;
index = index2;
}
}
}
/*节点数目为奇数时，同理*/
else
{
calnum = (calnum+1)/2;
if(my_rank > calnum)
{
MPI_Send(&index,1,MPI_INT,
my_rank-calnum,
my_rank-calnum,
MPI_COMM_WORLD);
MPI_Send(&num,1,MPI_DOUBLE,
my_rank-calnum,
my_rank-calnum,
MPI_COMM_WORLD);
}
else if(my_rank!=0 && my_rank < calnum)
{
MPI_Recv(&index2,1,MPI_INT,
my_rank+calnum,
my_rank,
MPI_COMM_WORLD,&status);
MPI_Recv(&num2,1,MPI_DOUBLE,
my_rank+calnum,
my_rank,
MPI_COMM_WORLD,&status);
if(num2 < num)
{
num = num2;
index = index2;
}
}
}
/*while一次循环之后都要进行同步*/
MPI_Barrier(MPI_COMM_WORLD);
}
/*广播最小值*/
MPI_Bcast(&index,1,MPI_INT,1,
MPI_COMM_WORLD);
MPI_Bcast(&num,1,MPI_DOUBLE,1,
MPI_COMM_WORLD);
/*各个进程根据广播的当前最小值num，index更新
自己的局部值*/
for(j=0;j<mynum;j++)
{
if((bdist[j]==FALSE)
&&(num + W(j,index) < dist[j])
)
dist[j] =num + W(j,index);
}
/*各个进程更新自己的bdist[]数组*/
if(my_rank == index/ep+1)
{
bdist[index%ep] = TRUE;
}
/*一次for循环同步一次*/
MPI_Barrier(MPI_COMM_WORLD);
}
}
4.2、并行算法的运行结果
实验平台：
以下CPU信息截图不全，但是从中可以看出2个cpu：
MPI环境搭建：
下载MPI包：mpich2-1.2.1p1.tar.tar 解压：tar –xvf mpich2-1.2.1p1.tar.tar
安装：./configure –prefix=/opt/mpi/
make
make install
编译源程序：
编译结果：
运行：./mpirun –np 1 shortest
运行出错，根据提示信息编辑配置文件：
Secretsword=****
再次运行，出错提示如下：
结果：未果。