并行的Dijkstra算法
要解决的问题
代码段
1.生成矩阵模块
2.主程序模块
3.输入模块
实验分析

要解决的问题

代码段
1.生成矩阵模块
#include "stdafx.h"
#include <stdlib.h>
#include <time.h>
#include <fstream.h>
#define Size 2000

int aMatrix[Size][Size];

int createMatrix()
{
	int i,j;
	srand(time(NULL));
	fstream fs;
	fs.open("Matrix2000.txt",ios::out);
	for(i=0;i<Size;i++)
	{
		for(j=0;j<Size;j++)
		{
			if(i==j)
			{
				aMatrix[i][j]=0;
			}
			else
			{
				aMatrix[i][j]=30+rand()%71;//1+rand()%100;
			}
		}
	}
	for(i=0;i<Size;i++)
	{
		for(j=0;j<i;j++)
		{
			aMatrix[i][j]=aMatrix[j][i];
		}
	}
	for(i=0;i<Size;i++)
	{
		for(j=0;j<Size;j++)
		{
			fs<<aMatrix[i][j]<<' ';
		}
	}
	fs.close();
	return 1;
}

int main(int argc, char* argv[])
{
	createMatrix();
	printf("OK!\n");	
	return 0;
}



2.主程序模块
#include <fstream.h>
#include <windows.h>
#include <stdio.h>
#include <process.h>
#include <string.h>
#include "stdafx.h"
#include <stdlib.h>
#include <time.h>


#define Size 2000
#define Threadnum 10//×î´óÏß³ÌÊý ÓÃÓÚ³õʼ»¯
#define Inf 1000000 

int xcs;//Óû§ÊäÈëµÄÏß³ÌÊý
int c1,c2;//c1,c2״̬

HANDLE hMutex;
struct ProcParam {
	int qidian;//ÿ´ÎÕÒµÄÆðµã
	int zhongdian;//ÿ´ÎÕÒµÄ×îºóÒ»¸öÊý
	int threadid;//Ïß³Ìid
	int *qbook;//´«bookÊý×éµÄÖ¸Õë
	int *qdst;//´«dstÊý×éµÄÖ¸Õë
	int *qmin;//´«jubuminµÄÖ¸Õë
	int pos;//¼Ç¼×îСֵµÄ×ø±ê
	int *qroute;//´«Â·¾¶routeÊý×éµÄÖ¸Õë
	int condition;//
	int *qtm;//
};


struct ProcParam2 {
	int start;
	int end;
	int sumtime;
	int zt;
};

int **aMatrix=new int *[Size];//´´½¨¾ØÕó
int readMatrix()
{
	int i,j;
	// ´´½¨¶¯Ì¬¶þάÊý×éaMatrix
	for(i=0;i<Size;i++)
	{
		aMatrix[i]=new int [Size];
	}
	fstream fs;
	fs.open("Matrix2000.txt",ios::in); //´ò¿ªÎļþmatrix.txt£¬¶ÁÈëÊý¾Ý
	for(i=0;i<Size;i++)
	{
		for(j=0;j<Size;j++)
		{
			fs>>aMatrix[i][j];
		}
	}
	fs.close();
 	return 1;
}

void time_Qiu(int *qiutime)//Çóʱ¼ä
{
	int t=0;
	for(int i=1;i<=xcs;i++)
	{
		printf("Ïß³Ì%d:¿ªÊ¼ÓÚ%d£¬½áÊøÓÚ%d,¹²ÓÃʱ¼äΪ%dms\n",i,qiutime[(i-1)*2+0],qiutime[(i-1)*2+1],qiutime[(i-1)*2+1]-qiutime[(i-1)*2+0]);
		qiutime[(i-1)*2+1]-qiutime[(i-1)*2+0]>t?(t=qiutime[(i-1)*2+1]-qiutime[(i-1)*2+0]):t;
	}
	printf("T=max{");
	for(i=1;i<=xcs;i++)
	{
		printf("t%d ",i);
	}
	printf("}=%dms\n",t);

}

unsigned __stdcall biJiao(void * arg)
{
	ProcParam *param;
	param = (ProcParam*)arg;
	if(param->condition==1) param->qtm[(param->threadid*2)+0]=GetTickCount();
	int min=Inf+1;//Èç¹û³öÏÖ¶¼ÊÇÎÞÇî¡£¡£´ý¶¨
	int pos;//ÊÇ·ñÒª¸³³õÖµ -1ÔõôÑù£¿
	for(int j=param->qidian;j<param->zhongdian;j++)
	{
		if(param->qbook[j] == 0 && param->qdst[j] < min)
		{  
			min = param->qdst[j];  
			pos=j;
		}  
	}
	param->qmin[param->threadid]=min;
	param->pos=pos;
	param->qtm[(param->threadid*2)+1]=GetTickCount();

    return 0;
}


int minFind(int *book,int *dst,int dyc,int *qiutime)//dycµÚÒ»´Î
{
	//int xcs=Threadnum;// ·ÖÅäµÄÏ̸߳öÊý
	int eachfind=Size/xcs;//ÿ¸öÏß³ÌËù×ö¸öÊý
	HANDLE hThread[Threadnum]; //HANDLE³ÆΪ¾ä±ú£¬±¾ÖÊÊÇÒ»¸ö×ÊÔ´µÄID±êʶ·û
	unsigned threadID[Threadnum];
	ProcParam *param=new ProcParam[Threadnum];
	int jubuMin[Threadnum];//¾Ö²¿×îСֵÊý×é
	for(int i=0;i<xcs;i++)//Êý×é³õʼ»¯Îª×î´óÖµ¡£
	{
		jubuMin[i]=Inf;
	}
	
	for(i=0;i<xcs;i++)
	{
		param[i].qtm=qiutime;
		param[i].condition=dyc;//ÅжÏÊÇ·ñÊǵÚÒ»´Î
		param[i].qmin=jubuMin;
		param[i].qbook=book;
		param[i].qdst=dst;
		param[i].qidian=0+eachfind*i;//Æðµã
		if(i!=xcs-1)
		{
			param[i].zhongdian=param[i].qidian+eachfind;//ÖÕµã
		}
		else
		{
			param[i].zhongdian=Size;
		}
		param[i].threadid=i;
		hThread[i] = (HANDLE)_beginthreadex(NULL, 0, biJiao, (void*)&param[i], 0, &threadID[i]);
	}
	WaitForMultipleObjects(xcs,hThread,TRUE,INFINITE);//µÈ´ýËùÓоä±ú¶¼Íê³ÉÔٱȽÏÆäÖеÄ×îСֵ  ÆäÖÐÏß³ÌÊýÒªÉ趨
	for(i=0;i<xcs;i++)
	{
		CloseHandle(hThread[i]);// ¹Ø±ÕÕâ¸öÏ̶߳ÔÓ¦µÄ¾ä±ú£¬ÊÍ·Å×ÊÔ´
	}
	int pos;//×îСֵËùÔÚµÄλÖÃ
	int minReal=Inf;//ËùÓоֲ¿×îСֵÖеÄ×îСֵ
	for(int j=0;j<xcs;j++)
	{
		if(jubuMin[j]<minReal)
		{
			minReal=jubuMin[j];
			pos=param[j].pos;
		}
	}
	return pos;
}

unsigned __stdcall refresh(void * arg)
{
	ProcParam *param;
	param = (ProcParam*)arg;
	for(int k=param->qidian;k<(param->zhongdian);k++)
	{
		if((param->qdst[k] > param->qdst[param->pos] + aMatrix[param->pos][k]) && param->qbook[k] == 0)
		{ 
			param->qroute[k]=param->pos;
			param->qdst[k] = param->qdst[param->pos] + aMatrix[param->pos][k];
			
		}    
	}
    return 0;
}
int printroute(int *route,int end,int start,int dst)
{
	int newroute[Size];//ÕýÐòµÄ·¾¶
	char str[Size]="";//¼Ç¼·¾¶µÄ×Ö·û´®
	int i=0;
	printf("v%d->v%dµÄ×î¶Ì·¾¶ÊÇ:\n",end,start);
	fstream fs;
	fs.open("print.txt",ios::app);//ÒÑ×·¼Ó·½Ê½´ò¿ª£¬ÈôÊÇios::out¾Í»áÖØдtxt
	fs<<"v"<<end<<"->"<<"v"<<start<<"×î¶Ì·¾¶ÊÇ£º"<<endl;
	
	while(route[end]!=-1)
	{
		newroute[i]=end;
		end=route[end];
		i++;
	}
	newroute[i]=end;
	if(end!=start)  
	{
		i++;
		newroute[i]=start;
	}
	

	for(int k=0;k<=i;k++)
	{
		if(k==i)
		{
			printf("v%d\n",newroute[k]);
			fs<<"v"<<newroute[k]<<endl;
		}
		else
		{
			printf("v%d->",newroute[k]);
			fs<<"v"<<newroute[k]<<"->";
		}
	}
	printf("¾àÀëΪ£º%d\n",dst);
	fs<<"¾àÀëΪ£º"<<dst<<endl;
	fs<<"      "<<endl;
	fs.close();
	return 1;
}

int reshuffle(int *book,int *dst,int &pos,int start,int end,int *route)//ºÃÏñÕâÀï²»ÓüÓ&
{
	//int xcs=Threadnum;
	int eachfind=Size/xcs;//ÿ¸öÏß³ÌËù×öµÄ¸öÊý

	HANDLE hThread[Threadnum]; //HANDLE³ÆΪ¾ä±ú£¬±¾ÖÊÊÇÒ»¸ö×ÊÔ´µÄID±êʶ·û
	unsigned threadID[Threadnum];
	ProcParam *param=new ProcParam[Threadnum];
	for(int i=0;i<xcs;i++)
	{
		param[i].qroute=route;
		param[i].qbook=book;
		param[i].qdst=dst;
		param[i].qidian=0+eachfind*i;//Æðµã
		if(i!=xcs-1)
		{
			param[i].zhongdian=param[i].qidian+eachfind;//ÖÕµã
		}
		else
		{
			param[i].zhongdian=Size;
		}
		param[i].pos=pos;
		param[i].threadid=i;
		hThread[i] = (HANDLE)_beginthreadex(NULL, 0, refresh, (void*)&param[i], 0, &threadID[i]);
	}
	WaitForMultipleObjects(xcs,hThread,TRUE,INFINITE);//µÈ´ýËùÓоä±ú¶¼Íê³ÉÔٱȽÏÆäÖеÄ×îСֵ  ÆäÖÐÏß³ÌÊýÒªÉ趨
	for(i=0;i<xcs;i++)
	{
		CloseHandle(hThread[i]);// ¹Ø±ÕÕâ¸öÏ̶߳ÔÓ¦µÄ¾ä±ú£¬ÊÍ·Å×ÊÔ´
	}
	return 0;
}



void dijkstra(int *book,int *dst,int start,int end,int &time,int option)//option 0´òӡ·¾¶  option1²»´òӡ·¾¶
{
	int qiutime[1000];//ÇóÒ»ÂÖÖеĸ÷¸ö²¿·Öʱ¼ä
	DWORD st,et;//starttime endtime
	st=GetTickCount();
	int i;
	for(i=0;i<Size;i++)
	{
		book[i]=0;
		dst[i]=aMatrix[start][i];
	}
	printf("\n");
	book[start]=1;
	int route[Size];
	for(i=0;i<Size;i++)
	{
		route[i]=-1;   //ÒòΪ³õʼËùÓнڵ㶼ԴÓÚstart£¬²»ÄÜÖ»ÊÇroute[start]=-1
	}

	for(i = 1; i < Size; i++)
	{
		int t=minFind(book,dst,i,qiutime);//×îСֵ×ø±ê
		book[t]=1;
		if(t==end)break;//²»ÓÃÈ«²¿±éÀúÕÒµ½¼´¿É
		reshuffle(book,dst,t,start,end,route);//startÆðµãendÖÕµãÊÇΪÁËÇ󷾶  route¼Ç¼µÚn¸ö½ÚµãµÄÇ°Ò»¸ö½ÚµãÊÇʲô
	}
	et=GetTickCount();
	time=et-st;
	if(!option)
	{
		printroute(route,end,start,dst[end]);//´òӡ·¾¶ºÍ·¾¶³¤¶È
		time_Qiu(qiutime);
	}

}

int jiexi(char *a,int size,int &start,int &end)//½âÎö×Ö·û´®£¬¶ÁÈ¡ÆäÖеÄÆðµãºÍÖÕµã
{
	int temp=-1,d=1,t=0;//temp=-1ÓÐv0
	for(int i=size-1;i>=0;i--)
	{
		if(48<=a[i]&&a[i]<=57)
		{
			if(48<=a[i+1]&&a[i]<=57)
			{
				d*=10;
			}
			else
			{
				d=1;
				t++;
				temp=0;
			}
			temp+=(a[i]-48)*d;
		}
		else
		{
			if(temp>-1)
			{
				if(t==1)
				{
					start=temp;
					temp=-1;
				}
				if(t==2)
				{
					end=temp;
				}
			}
		}
	}
	return 1;
}

int input(int *qi_zhong)
{
	FILE* f1=fopen("Input.txt","r");
	if(!f1) return 0;
	char buf[1024]={0};
	int mount=0;//¼Ç¼¹²ÓжàÉÙ¸öÆðÖÕµã
	while(!feof(f1))
	{
		int start,end;
		memset(buf, 0, sizeof(buf));//Çå¿Õ»º³åÇø
		fscanf(f1, "%s", buf);//´ÓÎļþÖжÁÈ¡Ò»¶ÎÊý¾Ý´æÈ뻺³åÇø
		jiexi(buf,strlen(buf),start,end);
		qi_zhong[mount++]=start;
		qi_zhong[mount++]=end;
	}
	fclose(f1);	
	return mount;
}

void thread_Time()
{
	for(xcs=2;xcs<11;xcs++)
	{
		int time,time_sum=0;
		int qi_zhong[100];//¼Ç¼ÆðÖÕµãµÄÊý×é
		int mount=input(qi_zhong);
		if(mount)
		{
			for(int i=0;i<mount/2;i++)
			{		
				int book[Size];
				int dst[Size];
				int start=qi_zhong[2*i];
				int end=qi_zhong[2*i+1];
				dijkstra(book,dst,start,end,time,1);
				time_sum+=time;
			}
		}
		printf("%d¸öỊ̈߳º%dms\n",xcs,time_sum);
		fstream fs;
		fs.open("time.txt",ios::app);//ios::out¾Í»áÖØдtxt
		fs<<xcs<<"¸öỊ̈߳º"<<time_sum<<endl;
		fs.close();
	}
}

unsigned __stdcall bXing(void * arg)
{
	ProcParam2 *param2;
	param2=(ProcParam2*)arg;
	int qi_zhong[100];//¼Ç¼ÆðÖÕµãµÄÊý×é
	input(qi_zhong);
	int time=0;
	for(int i=param2->start;i<param2->end;i++)
	{
		int book[Size];
		int dst[Size];
		dijkstra(book,dst,qi_zhong[2*i],qi_zhong[2*i+1],time,0);
		param2->sumtime+=time;
	}
	return 0;
}
void double_Find()
{
	HANDLE hThread[2];
	unsigned threadID[2];
	ProcParam2 *param2=new ProcParam2[2];
	int qi_zhong[100];//¼Ç¼ÆðÖÕµãµÄÊý×é
	int mount=input(qi_zhong);
	int n=mount/2;//¼¸¶ÔÆðÖÕµã
	int k=n/2;
	if(readMatrix())
	{
		for(int i=0;i<2;i++)
		{
			param2[i].start=0+i*k;
			param2[i].sumtime=0;
			if(i!=1)
			{
				param2[i].end=param2[i].start+k;
			}
			else
			{
				param2[i].end=n;
			}
			hThread[i] = (HANDLE)_beginthreadex(NULL, 0, bXing, (void*)&param2[i], 0, &threadID[i]);
		}
		WaitForMultipleObjects(2,hThread,TRUE,INFINITE);
		CloseHandle(hThread[0]);CloseHandle(hThread[1]);
		printf("Ëã·¨¹²ÓÃ=%dms\n",param2[0].sumtime+param2[1].sumtime);
	}
}

int chaxun()
{
	printf("****       ÇëÊäÈëÊ×Ä©½ÚµãºÅ         ****\n");
	int a,b;
	int start,end;
	scanf("%d %d",&a,&b);
	FILE* f1=fopen("print.txt","r");
	if(!f1) return 0;
	char buf[1024]={0};
	while(!feof(f1))
	{
		memset(buf, 0, sizeof(buf));//Çå¿Õ»º³åÇø
		fscanf(f1, "%s", buf);//´ÓÎļþÖжÁÈ¡Ò»¶ÎÊý¾Ý´æÈ뻺³åÇø
		jiexi(buf,strlen(buf),start,end);
		if(start==b&&end==a)
		{
			for(int i=0;i<strlen(buf);i++)
			{
				printf("%c",buf[i]);
			}
			printf("\n");
			memset(buf, 0, sizeof(buf));//Çå¿Õ»º³åÇø
			fscanf(f1, "%s", buf);
			for(i=0;i<strlen(buf);i++)
			{
				printf("%c",buf[i]);
			}
			printf("\n");
			memset(buf, 0, sizeof(buf));//Çå¿Õ»º³åÇø
			fscanf(f1, "%s", buf);
			for(i=0;i<strlen(buf);i++)
			{
				printf("%c",buf[i]);
			}
			printf("\n");
			return 1;
		}
		
	}
	printf("ÕÒ²»µ½\n");
	fclose(f1);	
}


void main()
{
	while(1)
	{
		printf("***************************************\n");
		printf("****        ²âÊÔÇëÊäÈë1            ****\n");
		printf("****        ͳ¼ÆÇëÊäÈë2            ****\n");
		printf("****        ²éѯÇëÊäÈë3            ****\n");
		printf("***************************************\n");
		int option;
		scanf("%d",&option);
		system("cls");
		switch(option)
		{
		case 1:
			printf("***************************************\n");
			printf("****        ÇëÊäÈëÏß³ÌÊý           ****\n");
			printf("***************************************\n");
			scanf("%d",&xcs);
			double_Find();break;
		case 2:thread_Time();break;
		case 3:chaxun();break;
		default:printf("ÇëÊäÈëÕýÈ·µÄÊý\n");
		}			
	getchar();
	}
}


3.输入模块
//input.txt
v8->v10
v231->v32
v342->v1034
v986->v128
v86->v1234

实验分析

译自：Data-Intensive Text Processing with MapReduce, Chap.5.2

一个有向图，由（V，E）组成，其中V是顶点的集合，E为联结各顶点的边，每条边e可能有相应的权重w。
图的表示方式有两种：邻接矩阵和邻接表。其中对于节点数较少的图，用邻接矩阵表示较为方便，计算时也能充分应用矩阵计算的一些优势。但是当节点数特别大，需要借助map-reduce计算时，用邻接表是更为合适的选择。每一行数据，key为NodeId，值为与这个节点邻接的所有节点的AdjacentList（可能还包含了每一个节点的权重）。

有向图的最短路径，即从源点s出发，到达图上任意结点的最短路径。
图论中最经典的最短路径算法即Dijkstra算法，它采用了贪心的策略，利用宽度优先一层层搜索，直到遍历所有结点或已经找到最短路径。算法伪代码如下：
Dijkstra(G,w, s)

    d[s] ← 0

    for all vertex v ∈ V do

        d[v] ← ∞

    Q ← {V }

    while Q != ∅ do

        u ←ExtractMin(Q)

        for all vertex v ∈ u.AdjacencyList do

            if d[v] > d[u] + w(u, v) then

                d[v] ← d[u] + w(u, v)

以下图的例子来说明这个算法：



n1为起始点。首先标记它到自已的距离为0，然后算法把所有节点(n2~n5)加入到优先队列Q中，并标记n1到这些点的距离为∞。
进入循环：
第一次迭代，从Q中取出最近的扩张点n1，找到n1的邻接点n2, n3，并更新到n2, n3的距离分别为10和5。
第二次迭代，从Q中取出最近的扩张点n3，找到n3的邻接点n2, n5, n4，这时发现n1经n3到n2的距离比直接到n2的距离要短，于是更新到n2的距离为8，到n5的距离为7，到n4的距离为14。
第三次迭代，从Q中取出最近的扩张点n5...
当所有点都搜索过后，算法终止，此时已经找到n1到各点的最短距离。
Dijkstra算法关键的一点是优先队列Q（实际实现可以用数组），它保存了全局的从源点出发最近的结点。而map-reduce则无法做到这一点（其实也不是做不到，就是比较麻烦，需要用distributed cache）。


基于map-reduce的并行算法跟Dijkstra算法有点类似，它也基于Dijkstra的迭代思想，伪代码如下：
class Mapper

    method Map(nid n, node N)

    d ← N.Distance

    Emit(nid n,N)                                  //Pass along graph structure [1]

    for all nodeid m ∈ N.AdjacencyList do

        Emit(nid m, d+w)                          //Emit distances to reachable nodes [2]


class Reducer

    method Reduce(nid m, [d1, d2, . . .])

    dmin←∞

    M ← ∅

    for all d ∈ counts [d1, d2, . . .] do

        if IsNode(d) then

            M ← d                                       //Recover graph structure

        else if d < dmin then                  //Look for shorter distance

            dmin ← d

    M.Distance← dmin                        //Update shortest distance

    Emit(nid m, node M)


它每次迭代执行一个map-reduce job，并且只遍历一个节点。在Map中，它先输出这个节点的完整邻接节点数据，即[1]。然后遍历该节点的邻接节点，并输出该节点ID及权重。在Reduce中，对当前节点m，遍历map的输出权重，若比当前的路径值小，则更新。最后输出该节点的路径值及完整邻接节点数据，作为下一次迭代的输入。
实现上有个细节需要注意的是，map的输出有两种类型的数据：邻接节点数据和权重数据，这可以通过一个包装类，并设置一个dataType变量来实现。
当遍历完所有的节点之后，迭代就终止了。
这个实现还有一个小问题就是，不知道完整的最短路径，这可以在[2]中修改一下Map的输出来实现。
下面是一个例子，邻接节点的数据格式为： nodeId --> distance | adj1: w1, adj2: w2...
原始输入为：
n1 --> 0 | n2: 6, n3: 15, n4:20

n2 --> ∞ | n3: 5, n4:8
n3 --> ∞ | n4:2
n4 -> ∞ |


第一次迭代：
Map：
n1 --> 0 | n2: 6, n3: 15, n4:20
n2 --> ∞ | n3: 5, n4:8

n3 --> ∞ | n4:2
n4 -> ∞ |
n2 --> 6
n3 --> 15
n4 --> 20
Reduce：
n1 --> 0 | n2: 6, n3: 15, n4:20
n2 --> 6 | n3: 5, n4:8
n3 --> 15 | n4:2
n4 --> 20 |
第二次迭代：
Map：

n1 --> 0 | n2: 6, n3: 15, n4:20
n2 --> 6 | n3: 5, n4:8
n3 --> 15 | n4:2
n4 --> 20 |
n3 --> 11
n4 --> 14
Reduce：

n1 --> 0 | n2: 6, n3: 15, n4:20
n2 --> 6 | n3: 5, n4:8
n3 --> 11 | n4:2
n4 --> 14 |
第三次迭代：

Map：

n1 --> 0 | n2: 6, n3: 15, n4:20
n2 --> 6 | n3: 5, n4:8
n3 --> 11 | n4:2
n4 --> 14 |
n4 --> 13
Reduce：

n1 --> 0 | n2: 6, n3: 15, n4:20
n2 --> 6 | n3: 5, n4:8
n3 --> 11 | n4:2
n4 --> 13 |
第四次迭代：
Map：

n1 --> 0 | n2: 6, n3: 15, n4:20
n2 --> 6 | n3: 5, n4:8
n3 --> 11 | n4:2
n4 --> 13 |
Reduce：

n1 --> 0 | n2: 6, n3: 15, n4:20
n2 --> 6 | n3: 5, n4:8
n3 --> 11 | n4:2
n4 --> 13 |
迭代终止

一个有向图，由（V，E）组成，其中V是顶点的集合，E为联结各顶点的边，每条边e可能有相应的权重w。

图的表示方式有两种：邻接矩阵和邻接表。其中对于节点数较少的图，用邻接矩阵表示较为方便，计算时也能充分应用矩阵计算的一些优势。但是当节点数特别大，需要借助map-reduce计算时，用邻接表是更为合适的选择。每一行数据，key为NodeId，值为与这个节点邻接的所有节点的AdjacentList（可能还包含了每一个节点的权重）。


有向图的最短路径，即从源点s出发，到达图上任意结点的最短路径。

图论中最经典的最短路径算法即Dijkstra算法，它采用了贪心的策略，利用宽度优先一层层搜索，直到遍历所有结点或已经找到最短路径。算法伪代码如下：

Dijkstra(G,w, s)
    d[s] ← 0
    for all vertex v ∈ V do
        d[v] ← ∞
    Q ← {V }
    while Q != ∅ do
        u ←ExtractMin(Q)
        for all vertex v ∈ u.AdjacencyList do
            if d[v] > d[u] + w(u, v) then
                d[v] ← d[u] + w(u, v)


以下图的例子来说明这个算法：


n1为起始点。首先标记它到自已的距离为0，然后算法把所有节点(n2~n5)加入到优先队列Q中，并标记n1到这些点的距离为∞。

进入循环：

第一次迭代，从Q中取出最近的扩张点n1，找到n1的邻接点n2, n3，并更新到n2, n3的距离分别为10和5。

第二次迭代，从Q中取出最近的扩张点n3，找到n3的邻接点n2, n5, n4，这时发现n1经n3到n2的距离比直接到n2的距离要短，于是更新到n2的距离为8，到n5的距离为7，到n4的距离为14。

第三次迭代，从Q中取出最近的扩张点n5...

当所有点都搜索过后，算法终止，此时已经找到n1到各点的最短距离。

Dijkstra算法关键的一点是优先队列Q（实际实现可以用数组），它保存了全局的从源点出发最近的结点。而map-reduce则无法做到这一点（其实也不是做不到，就是比较麻烦，需要用distributed cache）。


基于map-reduce的并行算法跟Dijkstra算法有点类似，它也基于Dijkstra的迭代思想，伪代码如下：

class Mapper
    method Map(nid n, node N)
    d ← N.Distance
    Emit(nid n,N)                                  //Pass along graph structure [1]
    for all nodeid m ∈ N.AdjacencyList do
        Emit(nid m, d+w)                          //Emit distances to reachable nodes [2]

class Reducer
    method Reduce(nid m, [d1, d2, . . .])
    dmin←∞
    M ← ∅
    for all d ∈ counts [d1, d2, . . .] do
        if IsNode(d) then
            M ← d                                       //Recover graph structure
        else if d < dmin then                  //Look for shorter distance
            dmin ← d
    M.Distance← dmin                        //Update shortest distance
    Emit(nid m, node M)

它每次迭代执行一个map-reduce job，并且只遍历一个节点。在Map中，它先输出这个节点的完整邻接节点数据，即[1]。然后遍历该节点的邻接节点，并输出该节点ID及权重。在Reduce中，对当前节点m，遍历map的输出权重，若比当前的路径值小，则更新。最后输出该节点的路径值及完整邻接节点数据，作为下一次迭代的输入。

实现上有个细节需要注意的是，map的输出有两种类型的数据：邻接节点数据和权重数据，这可以通过一个包装类，并设置一个dataType变量来实现。

当遍历完所有的节点之后，迭代就终止了。

这个实现还有一个小问题就是，不知道完整的最短路径，这可以在[2]中修改一下Map的输出来实现。

下面是一个例子，邻接节点的数据格式为： nodeId --> distance | adj1: w1, adj2: w2...

原始输入为：

n1 --> 0 | n2: 6, n3: 15, n4:20


n2 --> ∞ | n3: 5, n4:8

n3 --> ∞ | n4:2

n4 -> ∞ |


第一次迭代：

Map：

n1 --> 0 | n2: 6, n3: 15, n4:20

n2 --> ∞ | n3: 5, n4:8


n3 --> ∞ | n4:2

n4 -> ∞ |

n2 --> 6 
n3 --> 15

n4 --> 20

Reduce：

n1 --> 0 | n2: 6, n3: 15, n4:20

n2 --> 6 | n3: 5, n4:8

n3 --> 15 | n4:2

n4 --> 20 |

第二次迭代：

Map：


n1 --> 0 | n2: 6, n3: 15, n4:20

n2 --> 6 | n3: 5, n4:8

n3 --> 15 | n4:2

n4 --> 20 |

n3 --> 11 
n4 --> 14

Reduce：


n1 --> 0 | n2: 6, n3: 15, n4:20

n2 --> 6 | n3: 5, n4:8

n3 --> 11 | n4:2

n4 --> 14 |

第三次迭代：


Map：


n1 --> 0 | n2: 6, n3: 15, n4:20

n2 --> 6 | n3: 5, n4:8

n3 --> 11 | n4:2

n4 --> 14 |

n4 --> 13 
Reduce：


n1 --> 0 | n2: 6, n3: 15, n4:20

n2 --> 6 | n3: 5, n4:8

n3 --> 11 | n4:2

n4 --> 13 |

第四次迭代：

Map：


n1 --> 0 | n2: 6, n3: 15, n4:20

n2 --> 6 | n3: 5, n4:8

n3 --> 11 | n4:2

n4 --> 13 |

Reduce： 

n1 --> 0 | n2: 6, n3: 15, n4:20

n2 --> 6 | n3: 5, n4:8

n3 --> 11 | n4:2

n4 --> 13 |

迭代终止

求最短路径的串行算法在互联网上应该一搜一大堆，也非常简单，几行代码搞定。但Floyd的并行算法却很难搜到，github倒是有一些，但不容易运行成功，这里对这个算法的并行化进行详细的讲解，结合论文以及实际实现。


需要该算法的c++代码实现，请联系~~~

前言：对Floyd的介绍

求所有节点间的最短路径是一个基础的图问题，可以应用在生物信息学、社交网络以及交通规划。一个经典的求解方法就是Floyd-Warshal(FW)算法 5 6，它的空间复杂度是O(n^2)，时间复杂度是O(n^3)。基于块的算法的提出，有效利用了当前处理器内存的分层结构 7 8 9。然后在块算法的思路基础上，出现了很多多块CPU、多块GPU以及CPU+GPU结合的算法，但它们很多都只是用于计算最短路径距离矩阵，并不得到最短路径中间节点矩阵。

下文中的名词解释：

最短路径距离矩阵为：从节点i到节点j的最短距离；

最短路径中间节点矩阵为：从节点i到节点j，经过节点k最近。

1.Floyd的串行算法

贴一下代码，理解请看其他博客。

/***********在CPU中计算最短路径函数：floydWarshallCPUReference***********/
//void FloydWarshall::floydWarshallCPUReference(unsigned int *pathDistanceMatrix, unsigned int *pathMatrix, const unsigned int numNodes)
// * @param pathDistanceMatrix    输入，最短路径距离
// * @param pathMatrix            输入，最短路径中间节点存储矩阵
// * @param numNodes              输入，邻接矩阵节点个数
void floydWarshallCPUReference(float *pathDistanceMatrix, int *pathMatrix, const int numNodes)
{
    float distanceYtoX, distanceYtoK, distanceKtoX, indirectDistance;

    int width = numNodes;
    int yXwidth;

    //Floyd算法计算所有点之间的最短路径
    for (int k = 0; k < numNodes; ++k)
    {
        for (int y = 0; y < numNodes; ++y)
        {
            yXwidth = y * numNodes;
            for (int x = 0; x < numNodes; ++x)
            {
                distanceYtoX = pathDistanceMatrix[yXwidth + x];
                distanceYtoK = pathDistanceMatrix[yXwidth + k];
                distanceKtoX = pathDistanceMatrix[k * width + x];

                indirectDistance = distanceYtoK + distanceKtoX;

                if (indirectDistance < distanceYtoX)
                {
                    pathDistanceMatrix[yXwidth + x] = indirectDistance;
                    pathMatrix[yXwidth + x] = k;
                }
            }
        }
    }
};

2.论文阅读翻译

2.1 Blocked United Algorithm for the All-Pairs Shortest Paths Problem on Hybrid CPU-GPU Systems



论文题目：基于块融合的在混合CPU+GPU系统中求所有节点最短路径搜索算法

发表日期：2012年5月27日

作者：全是来自日本会津大学（Aizu）

第一作者的个人简历网址：https://u-aizu.ac.jp/~kazuya-m/index_en.html

第一作者的研究门简历网址：https://www.researchgate.net/profile/Kazuya_Matsumoto

评价：这篇文章即算最短路径距离矩阵，也算最短路径中间节点矩阵。同时发表时间相对较新，运算效率也非常不错

结果：非常不错，30000多节点，最快速度只需要2分钟！！！



解决的难点：

1.邻接矩阵的矩阵规模大于GPU本身的显存；

2.既算出最短路径距离矩阵，又算出最短路径中间节点矩阵;

3.邻接矩阵的矩阵规模n不是block因子b的整数倍；

算法原理：

1.具有两个核函数。第一个核函数解决block大小的子矩阵的子问题；另一个核函数是解决矩阵间的乘法问题。

2.为了高效利用GPU的多个内核，优化了CPU和GPU之间的通信；

3.将原有的n*n邻接矩阵拆分为b*b的距离矩阵和中间节点矩阵（子矩阵），b是block因子。（为了简化，文中的描述默认n是b的倍数）；

4.作者提出的算法如下：每次最外围的迭代和原串行FW算法的迭代次数相同，每次K迭代被划分为4个阶段。



阶段示例如下：



阶段1：更新b*b的中心线block{Dkk,Ckk}；

阶段2：更新b*b的与中心线block相关的列block{DIk,CIk} (I  !=  k)；

阶段3：更新b*b的与中心线block相关的行block{Dkj,Ckj} (j  !=  k)；

阶段4：更新b*b的其他与中心线block完全无关的block{Dij,Cij}(I  !=  k)且 (j  !=  k)；

5.为了求解阶段1的问题，可以使用传统的串行FW算法或者递归地调用基于Block的FW算法；（本文用了前者的方法，加入此算法的总算法如下图所示），算法4的4-16行运行在GPU中，算法4的合并操作17-20运行在CPU中。



6.阶段2-阶段4可以使用矩阵乘法更新，在本问题中，就是极小加代数。极小加代数的乘法和加法是分离执行的，极小加操作(MINPLUS)是运行在GPU中，矩阵加(MMA)运行在CPU中。这个操作减少了Z，C从CPU到GPU的数据传输，也就允许了CPU和GPU之间的高速通信。



具体的实现：

1.平台运行在单精度下；

2.运行环境为：Ubuntu 10.04   Linux核心2.6.32-37   gcc 4.6.2  -02优化等级

3.第一个kernel计算子APSP问题，划分减少了数据交互时间，因为D和C没有必要送给GPU。

算法4引入了另一个新的block因子bs，block因子的大小较大地影响着kernel的运行效率。为了在阶段1引入串行FW算法，每个外层迭代开始前，GPU线程间的同步是非常必要的。这意味着bs被共享内存的大小限制着。本文将bs=64，给两个block   (Dck,ck),(Cck,ck)需要的大小(32KB=2*64^2*4字节)等于共享内存的大小。bs取的小了，就会造成性能恶化。

4.第二个kernel计算矩阵乘法问题，本文设计了SGEMM核（单精度通用矩阵乘法） 25 26  。在此核函数中，一个基本的乘加操作Zi,j=Zi,j+Xi,k*Yk,j可以在GPU中用单FMA实现。但是在极小加核函数中，4个不同的指令才能实现这个操作。这意味着，对于相同的矩阵大小，SGEMM核函数的效率是MINPUS核函数的4倍。这四个指令是：一个加法指令，一个比较指令以及两个复制指令。

解决问题的技巧：

1.解决邻接矩阵的矩阵规模n不是block因子b的整数倍

           将邻接矩阵补充数据，补为block因子b的整数倍，这样会有性能的下降（进行了无用的计算）。为了避免最坏情况下的补充，即比b的整数倍多一个，需要补充b-1个行和列的数据，采用了最忧block因子的选取方法。

n_pad=[(n+b-1)/b]*b

上式中，n_pad为补充后的邻接矩阵的矩阵规模；n为当前邻接矩阵的矩阵规模；b为block因子；这里的数据都是整数型，因此每次计算都有取整操作。

补充的数据，邻接矩阵的值都要设为无穷大（也就是相对于正常的数据非常大的数）。

这个计算最优block因子b的方法，本文未给出。

2.选择另一个block因子bs的方法

根据共享内存的大小，决定block因子bs的大小。

如果共享内存大小为32kB，每次计算需要两个矩阵（需x2），单精度（需x4），则根据如下的公式

32kB=2*4*n^2

得到矩阵的规模大小，即block因子bs的大小为64.

3.对CPU到GPU的数据传输大小进行估计

对于算法3，Phase1.    

有一个距离矩阵输入，一个距离矩阵和一个节点矩阵输出，即3个矩阵；

单精度计算，因此一个数据4个字节；

矩阵的规模为block因子b。因此数据传输大小为：

4*3b^2

对于Phase2和Phase3.

有一个依赖的中心距离矩阵输入（只需要输入距离矩阵，且对其他矩阵通用）；

每行或者列，共(n/b-1)个矩阵。每个矩阵代表有一个距离矩阵输入，一个距离矩阵和一个节点矩阵输出，即3个矩阵；

单精度计算，因此一个数据4个字节；

矩阵的规模为block因子b。因此数据传输大小为：

4*(b^2+3b^2(n/b-1))

对于Phase4.

有一个依赖的中心距离矩阵输入（只需要输入距离矩阵，且对其他矩阵通用）；

同时有行或者列，共(n/b-1)^2个矩阵。每个矩阵代表有一个距离矩阵输入，一个距离矩阵和一个节点矩阵输出，即3个矩阵；

单精度计算，因此一个数据4个字节；

矩阵的规模为block因子b。因此数据传输大小为：

4*(b^2+3b^2*(n/b-1)^2)

4.邻接矩阵的矩阵规模大于GPU本身的显存；

            原来需要传输4*3n^2个字节的数据，现在只需要4*3b^2个字节的数据，从而解决了矩阵规模大于GPU显卡的问题。同时，b越小，越能解决这个问题。

5.算法性能的评价

             带宽的算法为：

BW[字节/s]=GPU峰值计算能力[Flop/s]/计算密集性[Flops/字节]

GPU峰值计算能力是GPU的固有属性；

计算密集性为总共的浮点数操作与数据传输大小的比值，计算方法如下：

计算密集性=4b^3/(4*3b^2)=b/3

因此，使用相当大的block因子b，可以隐藏数据的传输时间（hidden communication latencies）.但不能无限大的增加

block因子b的大小，以为太大就饱和了，带宽到达极限后就不怎么增加了。对文中的机器，给出的值为b=1536.

 

完结撒花~~~~

有什么不懂的评论区见

需要源码的请联系购买~~~