MPI並行編程系列之五：圖的單源最短路徑演算法

    圖的單源最短路徑問題是指求一個從指定的頂點s到其他所有頂點i之間的最短距離。因為是從一點到其他頂點的距離，所以稱謂單源。本文將介紹dijkstra演算法：按路徑長度遞增次序的最短路徑演算法，並給出了對該演算法MPI並行化演算法。

一、演算法思想 

    該演算法的思想是：引入一個輔助向量D，它的每個分量D[i]為源頂點v到其他頂點v[i]的路徑長度。初始態為：如有從v到vi的路徑，則D[i]為弧[v,vi]的權值；否則D[i]為無窮大。顯然

                  D[j] = min{D[i]}

為頂點v出發到其他頂點的一條最短路徑的長度，其路徑為（v，vj）。

    那麼下一條最短路徑長度如何計算呢？其實很簡單，下一條最短路徑長度要麼是源頂點v直接到某一頂點vk的長度，即{v，vk}。要麼是源頂點v經過頂點vj到某一頂點的長度，即{v，vj，vk}。

   我們假設S為已經求得最短路徑的頂點的集合，那麼下一條最短路徑（設其終點為x），要麼是弧{v， vx}，要麼為中間只經過S中頂點而最後到達終點X的路徑。

    在一般情況下，下一條最短路徑的長度為：

                 D[j] = min{D[i] | vi 屬於 V-S}

    其中V為圖頂點的集合， D[i]為弧{v， vi}的權值，或者為D[k]和弧{vk， vi}權值之和。

二、演算法描述

     根據以上思想，我們得到演算法描述如下：

輸入：圖G的鄰接矩陣m[0...n-1][0...n-1]，源頂點 v

輸出：最短路徑值向量D[0...n-1], 最短路徑矩陣p[0...n-1][0...n-2].其中D[i]為源頂點v到頂點vi的最短路徑長度，向量p[i]為源頂點v到頂點vi的最短路徑

1）初始化D[i]

     D[i] = m[v][vi] == 無窮大 ? 無窮大 : m[v][vi]

2)計算當前最短路徑值

     min = min{D[i]}

     final[i] = 1   //標記頂點i已經取得最短路徑

3）更新最短路徑值及最短路徑

     for(i = 0; i < n; i++)

         if(!final[i])

             if(D[i] > min + m[vk][vi])

                 D[i] = min + m[vk][vi];

             end if

          endif

         for( j = 0; j < n-1 ; j++)

               if(p[vk][j] != 無窮大)

                   p[vi][j] = p[vk][j];

              end if

         end for

         p[vi][j] = vi

     end for

4) 輸出最短路徑和最短路徑值

三、演算法實現

     這裡只給出演算法的核心代碼， 如下：

   1:  void short_path_function(

   2:          int **matrix,

   3:          int vertex_num,

   4:          int vertex_source){

   5:   

   6:      int *short_path_value;

   7:      int *final;

   8:   

   9:      int *short_path_storage;

  10:      int **short_path;

  11:   

  12:      int min;

  13:      int vertex;

  14:   

  15:      int i,j,k;

  16:   

  17:      //分配空間

  18:      short_path_value = my_malloc(sizeof(int) * vertex_num);

  19:      final = my_malloc(sizeof(int) * vertex_num);

  20:      dynamic_allocate_matrix((void *)&short_path, (void *)&short_path_storage, vertex_num,

  21:              vertex_num-1, sizeof(int));

  22:   

  23:      //初始化

  24:      for(i = 0; i < vertex_num; i++){

  25:          final[i] = 0;

  26:          short_path_value[i] = matrix[vertex_source][i];

  27:   

  28:          //設定空路徑

  29:          for(j = 0; j < vertex_num-1; j++)

  30:              short_path[i][j] = -1;

  31:   

  32:          //初始化路徑

  33:          if(short_path_value[i] < MAX_VAULE)

  34:              short_path[i][0] = i;

  35:      }

  36:   

  37:      final[vertex_source] = 1;

  38:      for(i = 1; i < vertex_num; i++){

  39:          //找出從源頂點出發的一條最短路徑長度頂點

  40:          min =MAX_VAULE ;

  41:          for(j = 0; j < vertex_num; j++)

  42:              if(!final[j])

  43:                  if(short_path_value[j] < min){

  44:                      min = short_path_value[j];

  45:                      vertex = j;

  46:                  }

  47:   

  48:          final[vertex] = 1;

  49:   

  50:          //跟新最短路徑

  51:          for(j = 0; j < vertex_num; j++){

  52:              if(!final[j])

  53:                  if(min + matrix[vertex][j] < short_path_value[j]){

  54:                      //跟新最短路徑長度

  55:                      short_path_value[j] = min + matrix[vertex][j];

  56:   

  57:                      //更新路徑

  58:                      k = 0;

  59:                      while(short_path[vertex][k] != -1){

  60:                          short_path[j][k] = short_path[vertex][k];

  61:                          k++;

  62:                      }

  63:   

  64:                      short_path[j][k] = j;

  65:                  }

  66:          }

  67:      }

  68:   

  69:      //列印結果

  70:      array_int_print(vertex_num, short_path_value);

  71:      print_int_matrix(short_path, vertex_num-1, vertex_num);

  72:  }

四、並行演算法描述

      我們對這個演算法進行並行化分析，顯然初始化向量D（對應於演算法的第一步），更新最短路徑值和最短路徑（對應於演算法的第三步）是可以並行化實現的，因為只要有了當前最短路徑和最短路徑值，各個頂點的最短路徑的演算法是相互獨立的。在本並行化演算法中，如何求得當前的最短路徑和最短路徑值成為關鍵。

      我們假設一共用P個進程，圖右n個頂點，我們讓每個進程負責n/p個頂點，每個進程都有自己的向量D和最短路徑p，我們如何求得當前的最短路徑長度和最短路徑呢？首先，我們可以計算出各個進程的當前最短路徑，並將局部最短路徑發往進程0，進程0對這些進程的最短路徑進行比較，取得當前的全域最短路徑，並將這一結果廣播到所有的進程。

      其演算法描述如下：

      我們假設總共有p個進程

      輸入：圖G的鄰接矩陣m[0...n-1][0...n-1]，源頂點 v

      輸出：最短路徑值向量D[0...n-1], 最短路徑矩陣p[0...n-1][0...n-2].

    1） 進程0讀取鄰接矩陣m和源節點v，並將m，v廣播到其他所有的進程

    2） 各進程並行初始化各自的局部D和P

    3）求最短路徑

         1）各進程並行計算出各自的局部最短路徑值，並將其發送0號進程

         2）0號進程求出全域最短路徑值和對應的進程號，並將其廣播到其他所有的進程。

         3）擁有全域最短路徑值的進程將其對應的最短路徑廣播到其他進程（用於更新各個進程的最短路徑）

    4）各進程並行跟新各自的D和P

五、 並行演算法實現

    我們這裡只給出演算法中第三步的具體實現：

   1:  int get_short_path_value(

   2:          int *final,

   3:          int *vertex,

   4:          int **short_path,

   5:          int *short_path_value,

   6:          int *short_path_copy,

   7:          int *shortest,

   8:          int vertex_num,

   9:          int local_vertex_num,

  10:          int process_id,

  11:          int process_size,

  12:          MPI_Comm comm){

  13:   

  14:      int min;

  15:      int local_vertex;

  16:   

  17:      int shortest_process_id;

  18:   

  19:      MPI_Status status;

  20:   

  21:      int j;

  22:      

  23:      //取到每個進程的路徑的最小值

  24:      min = MAX_VAULE;

  25:      for(j = 0; j < local_vertex_num; j++)

  26:          if(!final[j])

  27:              if(short_path_value[j] < min){

  28:                  min = short_path_value[j];

  29:                  local_vertex = j;

  30:              }

  31:   

  32:      //各進程將自己的最小路徑值發往0號進程

  33:      if(process_id)

  34:           MPI_Send(&min, 1, MPI_INT, 0, DATA_MESSAGE, comm); 

  35:      else{

  36:          //0號進程收集所有進程的最小值

  37:          shortest[0] = min;

  38:          for(j = 1; j < process_size; j++)

  39:              MPI_Recv(shortest+j, 1, MPI_INT, j, DATA_MESSAGE, comm, &status);

  40:   

  41:           //0號進程計算出全域最小值，並廣播到其他各個進程

  42:            min = shortest[0];

  43:            shortest_process_id = 0;

  44:            for(j = 1; j < process_size; j++)

  45:               if(shortest[j] < min){

  46:                   min = shortest[j];

  47:                   shortest_process_id = j;

  48:                }

  49:      }

  50:      //0號進程將最小值廣播到所有進程,如果最小值為無窮大，則返回

  51:      MPI_Bcast(&min, 1, MPI_INT, 0, comm);

  52:      if(min == MAX_VAULE)

  53:          return min;

  54:   

  55:      //0號進程將最小值所屬的進程廣播到所有的進程

  56:      MPI_Bcast(&shortest_process_id, 1, MPI_INT, 0, comm);

  57:   

  58:      //擁有最小值的進程，將其對應的標誌位標1，並計算該座標的全域座標

  59:      if(process_id == shortest_process_id){

  60:          final[local_vertex] = 1;

  61:          *vertex = BLOCK_LOW(shortest_process_id, process_size, vertex_num) + local_vertex;

  62:   

  63:          //複製最短路徑

  64:          for(j = 0; j < vertex_num -1; j++)

  65:              short_path_copy[j] = short_path[local_vertex][j];

  66:      }

  67:   

  68:      //廣播最短路徑

  69:      MPI_Bcast(short_path_copy, vertex_num-1, MPI_INT, shortest_process_id,comm);

  70:      MPI_Bcast(vertex, 1, MPI_INT, shortest_process_id, comm);

  71:   

  72:      return min;

  73:  }

六、MPI函數說明

     本並行演算法用到的MPI函數主要為點對點通訊函數MPI_Send和MPI_Recv，全域通訊函數MPI_Bcast。MPI名字起的真好，訊息傳遞編程，我覺得ＭＰＩ演算法的核心就是任務的劃分和任務間的通訊。

　　下篇將介紹圖的最下產生樹及其並行演算法。