燕山大學(xué)并行處理與體系結(jié)構(gòu)—算法實(shí)踐部分---單源最短路徑Dijkstra專業(yè):班級(jí):姓名:學(xué)號(hào):一．求解問(wèn)題描述：?jiǎn)卧醋疃搪窂絾?wèn)題即指：已知一個(gè)n結(jié)點(diǎn)有向圖G=(V,E)和邊的權(quán)函數(shù)c(e)，求由G中某指定結(jié)點(diǎn)v0到其他各個(gè)結(jié)點(diǎn)的最短路徑。這里還假定所有的權(quán)值都是正的。1.單源最短路徑的串行算法為了制定產(chǎn)生最短路徑的貪心算法的基礎(chǔ)算法，對(duì)于這個(gè)問(wèn)題應(yīng)該想出一個(gè)多級(jí)解決辦法和一個(gè)最優(yōu)的度量標(biāo)準(zhǔn)。方法之一是逐條構(gòu)造這些最短路徑，可以使用迄今已生成的所有路徑長(zhǎng)度之和作為一種度量，為了使這一度量達(dá)到最小，其單獨(dú)的每一條路徑都必須有最小長(zhǎng)度。使用這一度量標(biāo)準(zhǔn)的時(shí)候，假定已經(jīng)構(gòu)造了I條最短路徑，則下面要構(gòu)造的路徑應(yīng)該使下一條最短的最小長(zhǎng)度路徑。生成從v0到所有其他結(jié)點(diǎn)的最短路徑的貪心算法就是按照路徑長(zhǎng)度的非降次序生成這些路徑。首先，生成一條到最近結(jié)點(diǎn)的最短路徑。然后生成一條到第二近結(jié)點(diǎn)的最短路徑等等。為了按這樣的次序生成最短路徑，就需要確定與其生成最短路徑的下一個(gè)結(jié)點(diǎn)，以及到達(dá)這個(gè)結(jié)點(diǎn)的最短路徑。設(shè)S表示對(duì)其已經(jīng)生成了最短路徑的那些結(jié)點(diǎn)（包括v0）的集合。對(duì)于不在S中的w，設(shè)DIST(w)是從v0開(kāi)始只經(jīng)過(guò)S中的結(jié)點(diǎn)而在w結(jié)束的那條最短路徑的長(zhǎng)度，則有：如果下一條最短路徑是到結(jié)點(diǎn)u，則這條路徑是從v0處開(kāi)始而在u處終止，并且只通過(guò)那些在S中的結(jié)點(diǎn)。所生成的下一條路徑的終點(diǎn)u必須是所有不在S內(nèi)的結(jié)點(diǎn)中具有最小距離DIST(u)的結(jié)點(diǎn)。在像2)中那樣選出了結(jié)點(diǎn)u并生成了從v0到u的最短路徑之后，結(jié)點(diǎn)u就成為S中的一個(gè)成員。這可由DIST的定義和上述1)而得到。在存在幾個(gè)不在S中而又有相同DIST的結(jié)點(diǎn)的情況下，則可以選擇這些結(jié)點(diǎn)中的任何一個(gè)。在像2)中那樣選出了結(jié)點(diǎn)u并生成了從v0到u的最短路徑之后，結(jié)點(diǎn)u就成為了S中的一個(gè)成員。此時(shí)，那些從v0開(kāi)始，只通過(guò)S中的結(jié)點(diǎn)并且在S外的結(jié)點(diǎn)w處結(jié)束的最短路徑可能會(huì)減小。即DIST(w)的值可能改變。因此可以得出結(jié)論，如果DIST(w)會(huì)減少，那是由于有一條從v0經(jīng)u到w更短路徑得緣故，其中從v0到u的路徑是這樣一條最短路徑,而從u到v的路徑是邊<u,w>。這條路徑的長(zhǎng)度是DIST(u)+c(u,w)。由上述可知單源最短路徑問(wèn)題存在一個(gè)簡(jiǎn)單算法。這個(gè)算法（通稱Dijkstra算法）實(shí)際上只求出從v0到G中的所有其他結(jié)點(diǎn)的最短路徑長(zhǎng)度。下面給出一個(gè)Dijkstra算法。輸入：邊權(quán)鄰接矩陣W(約定頂點(diǎn)i,j之間無(wú)邊連接時(shí)w(i,j)的值為無(wú)窮大，且w(i,j)＝0，待計(jì)算頂點(diǎn)的標(biāo)號(hào)s輸出：dist(0:N-1),其中dist(i)表示頂點(diǎn)s到頂點(diǎn)i的最短路徑（1<=i<=N）procedureDijkstraBegindist(s)=0fori=0toN-1doif≠sthendist(i)=∞endforVL=Vfori=0toN-1do(4.1)從VL中找一個(gè)頂點(diǎn)u,使得dist(u)最小(4.2)for(每個(gè)頂點(diǎn)v∈VL)and(<u,v>∈E)doifdist(u)+w(u,v)<dist(v)thendist(u)+w(u,v)endifendforVL=VL-{u}endforEnd2.單源最短路徑并行算法Dijkstra算法1．Dijkstra算法基本原理：假定有一個(gè)待搜索的頂點(diǎn)表VL，初始化時(shí)做：dist(s)0,dist(i)=∞（i≠s）,VL=V。每次從VL（非空集）中選取這樣的一個(gè)頂點(diǎn)u,它的dist(u)最小。將選出的u點(diǎn)作為搜索頂點(diǎn)，對(duì)于其他還在VL內(nèi)的頂點(diǎn)v,若〈u，v〉∈E，而且dist(u)+w<u,v><dist(v),則更新dist(v)為dist(u)+w<u,v>，同時(shí)從VL中將u刪除，直到VL成為空集時(shí)算法終止。2．Djkstra并行算法：1）每個(gè)處理器分配n=N/p各節(jié)點(diǎn)進(jìn)行初始化（N為節(jié)點(diǎn)數(shù)，p為處理器個(gè)數(shù)）。2）首先每個(gè)處理器分配n個(gè)節(jié)點(diǎn)分別求局部最小值，然后再p個(gè)處理器和作求全局最小值，最后再將這一最小值廣播出去。P個(gè)處理器合作方法：當(dāng)p為偶數(shù)時(shí)，后p/2個(gè)處理器將自己的局部最小值分別發(fā)送到對(duì)應(yīng)的前p/2個(gè)處理器中，由前p/2個(gè)處理器比較出2個(gè)局部最小值中相對(duì)較小者并保留。當(dāng)p為奇數(shù)時(shí)，設(shè)p=2h+1,則后h個(gè)處理器的值分別發(fā)送到前h個(gè)處理器中，比較并保留最小值。一次這樣的循環(huán)過(guò)程后，p個(gè)最小值減少了一半，兩次后，變?yōu)?/4，如此一層一層的比較，logp次循環(huán)后，就可以得出唯一的全局最小值，然后將其廣播出去。3）每個(gè)處理器分配n個(gè)頂點(diǎn)，然后獨(dú)立進(jìn)行更新dist的操作。其中廣播實(shí)現(xiàn)步驟如下：輸入：數(shù)據(jù)u(存放在單元B(1)中)；輸出：將u廣播到數(shù)組B的所有單元中去。BeginB(1)←u；（2）fori←0to-1do(3) foreachj:2i+1≤j≤2i+1pardoB(j)←B(j-2j)EndforEndforEnd.3．Djkstra算法復(fù)雜性分析：根據(jù)以上步驟可得并行算法使用了p個(gè)處理器，其時(shí)間復(fù)雜度為O(N2/p+Nlogp)。這里應(yīng)該注意的是處理器0只進(jìn)行輸入/輸出工作，不參與任何其他的計(jì)算，因此實(shí)際參加運(yùn)算的處理器為p-1個(gè)，所以有n=N/(p-1)。二．并行MPI源代碼#include<stdio.h>#include<stdlib.h>#include<malloc.h>#include<math.h>#include<mpi.h>#defineA(i,j)A[i*N+j]#defineMST(i,j)MST[i*n*p+j]#defineMAX1000intN;intn;intp;inti,j,k;doublel;int*D,*C,*W,*J;int*A;int*MST;intmyid;MPI_Statusstatus;/**輸出結(jié)果MST**/voidprintM(){ inti,j; if(myid==0){ printf("theMSTis:\n"); for(i=0;i<N;i++) { for(j=0;j<N;j++) printf("%d",MST[i*n*p+j]); printf("\n"); } }}/**輸出數(shù)組D內(nèi)容**/voidprintD(){ inti; if(myid==0){ printf("thearrayDis:\n"); for(i=0;i<N;i++) { printf("%d",D[i]); } printf("\n"); }}/**讀入鄰接矩陣**/intreadA(){ inti,j; intp1; p1=p; printf("Inputthesizeofmatrix:"); scanf("%d",&N); n=N/p1; if(N%p1!=0)n++; A=(int*)malloc(sizeof(int)*(n*p1)*N); if(A==NULL){ printf("Errorwhenallocatingmemory\n"); exit(0); } for(i=0;i<N;i++) for(j=0;j<N;j++) scanf("%d",&(A(i,j))); for(i=N;i<n*p1;i++) for(j=0;j<N;j++) A(i,j)=MAX-1; p=p1; return(0);}/**廣播特定數(shù)組**/voidbcast(int*P){MPI_Bcast(P,N,MPI_INT,0,MPI_COMM_WORLD);}/**求最小的數(shù)學(xué)函數(shù)**/intmin(inta,intb){return(a<b?a:b);}/**檢查MST是否已完成**/intconnected(){ inti; intflag; flag=1; for(i=1;i<N;i++) if(D[i]!=D[0]) { flag=0; break; } return(flag);}/**為各頂點(diǎn)找出距離最近的樹(shù)**/voidD_to_C(){ inti,j; for(i=0;i<n;i++){ W[n*myid+i]=MAX; for(j=N-1;j>=0;j--) if((A(i,j)>0)&&(D[j]!=D[n*myid+i])&&(A(i,j)<=W[n*myid+i])) { C[n*myid+i]=D[j]; W[n*myid+i]=A(i,j); J[n*myid+i]=j; } if(W[n*myid+i]==MAX)C[n*myid+i]=D[n*myid+i];}}/**為各樹(shù)找出外連最短邊**/voidC_to_C(){ inttempw,tempj; for(i=0;i<n;i++){ tempj=N+1; tempw=MAX; for(j=N-1;j>=0;j--) if((D[j]==n*myid+i)&&(C[j]!=n*myid+i)&&(W[j]<=tempw)) { C[myid*n+i]=C[j]; tempw=W[j]; tempj=j; } if(myid==0) { if((tempj<N)&&(J[tempj]<N)) MST(tempj,J[tempj])=MST(J[tempj],tempj)=tempw; for(j=1;j<p;j++) { MPI_Recv(&tempj,1,MPI_INT,j,j,MPI_COMM_WORLD,&status); MPI_Recv(&tempw,1,MPI_INT,j,0,MPI_COMM_WORLD,&status); if((tempj<N)&&(tempw<N))MST(tempj,tempw)=MST(tempw,tempj)=A(tempj,tempw); } } else { MPI_Send(&tempj,1,MPI_INT,0,myid,MPI_COMM_WORLD); MPI_Send(&J[tempj],1,MPI_INT,0,0,MPI_COMM_WORLD); } MPI_Barrier(MPI_COMM_WORLD); }}/**調(diào)整數(shù)組頂點(diǎn)的標(biāo)識(shí)**/voidCC_to_C(){ for(i=0;i<n;i++) C[myid*n+i]=C[C[myid*n+i]];}/**產(chǎn)生新樹(shù)**/voidCD_to_D(){ for(i=0;i<n;i++) D[myid*n+i]=min(C[myid*n+i],D[C[myid*n+i]]);}/**釋放動(dòng)態(tài)內(nèi)存**/voidfreeall(){free(A);free(D);free(C);}l/**主函數(shù)**/intmain(intargc,char**argv){ inti,j,k; intgroup_size; MPI_Init(&argc,&argv); MPI_Comm_size(MPI_COMM_WORLD,&group_size); MPI_Comm_rank(MPI_COMM_WORLD,&myid); p=group_size; /**處理器0讀入鄰接矩陣**/ if(myid==0) {readA(); MST=(int*)malloc(sizeof(int)*n*p*n*p); } 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)); W=(int*)malloc(sizeof(int)*(n*p)); J=(int*)malloc(sizeof(int)*(n*p)); if(myid!=0) A=(int*)malloc(sizeof(int)*n*N); /**數(shù)組D初始化，步驟(1)**/ for(i=0;i<n;i++)D[myid*n+i]=myid*n+i; MPI_Gather(&D[myid*n],n,MPI_INT,D,n,MPI_INT,0,MPI_COMM_WORLD); bcast(D); /**處理器0向其它處理器發(fā)送必要信息**/ 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); i=1; /**主循環(huán)，步驟(2)**/ while(!connected()){ if(myid==0)printf("loop%d\n",i); printD(); printM(); /**步驟(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); MPI_Gather(&W[n*myid],n,MPI_INT,W,n,MPI_INT,0,MPI_COMM_WORLD); bcast(C); bcast(W); MPI_Barrier(MPI_COMM_WORLD); /**步驟(2.2)**/ C_to_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); bcast(D); i++; } if(myid==0)printf("loop%d\n",i); printD(); printM(); freeall(); MPI_Finalize(); return(0);}三.運(yùn)行結(jié)果：本實(shí)例中使用了5個(gè)處理器。mpirun–np5mst運(yùn)行結(jié)果：Inputthesizeofmatrix:9045862573508316495720641389456057216458702347414750784789321067421435809532897620最終輸出結(jié)果：theMSTis:000002000000011000000001012000000010010002000211020100000001000001100000002000000說(shuō)明：輸出結(jié)果為最小生成樹(shù)的邊。四．Dijkstra的單源最短路徑并行算法參數(shù)n存儲(chǔ)圖中的頂點(diǎn)總數(shù)，參數(shù)source存儲(chǔ)計(jì)算單元最短路徑的源頂點(diǎn)，參數(shù)wgt指向圖的加權(quán)鄰接矩陣的本地存儲(chǔ)部分。參數(shù)lengths指向一個(gè)向量，該向量將用來(lái)存儲(chǔ)從源頂點(diǎn)到本地存儲(chǔ)頂點(diǎn)的最短路徑。最后，參數(shù)comm是將被MPI理性程序使用的通信器。這個(gè)例行程序假設(shè)頂點(diǎn)數(shù)目是處理器數(shù)目的整數(shù)倍。SingleSource(intn,intsource,int*wgt,int*lengths,MPI_Commcomm){inti,j;intnlocal;int*marker;intfirstvtx;intlastvtx;intu,udist;intlminpair[2],gminpair[2];intnpes,myrank;MPI_Statsstatus;MPI_Comm_size(comm.,&npes);MPI_Comm_rank(comm.,&myrank);nlocal=n/npes;firstvtx=myrank*nlocal;lastvtx=firstvtx+nlocal-1;for(j=0;j<nlocal;j++)lengths[j]=wgt[source*nlocal+j];marker=(int*)malloc(nlocal*sizeof(sizeof(int));for(j=0;j<nlocal;j++);marker[j]=1;if(source>=firstvtx&&sorce<=lastvtx)marker[source-firstvtx]=0;for(i=1;i<n;i++){lminpair[0]=MAXINT;lminpair[1]=-1;for(j=0;j<nlocal;j++{})if(marker[j]&&lengths[j]<lminpair[0]{lminpair[0]=lengths[j];lminpair[1]=firstvtx+j;})}MPI_Allreduce(lminpair,gminpair,1,MPI_2INT,MPI_MINLOC,comm);udist=gminpair[0];u=gminpair[1];if(u==lminpair[1])marker[u-firstvtx]=0;for(j=0;j<nlocal;j++){if(marker[j]&&udist+wgt[u*nlocal+j]<lengths[j])lengths[j]=udist+wgt[u*nlocal+j];}}free(marker);}五．算法性能及改進(jìn)Dijkstra的并行單源最短路