求最短路径Floyd算法的并行化（解APSP问题）

求最短路径的串行算法在互联网上应该一搜一大堆，也非常简单，几行代码搞定。但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.

 

完结撒花~~~~

有什么不懂的评论区见

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