BFS：Breadth-First Search，广度优先搜索
SpMV：Sparse-Matrix (Dense-)Vector Multiplication，稀疏矩阵向量乘法
SpMSpV：Sparse-Matrix Sparse-Vector Multiplication，稀疏矩阵稀疏向量乘法
基于CSR的SpMV计算方法
CSR（Compressed Sparse Row，压缩行存储）是存储稀疏矩阵的一种有效方式，避免了使用二维数组方式时存储大量0值的情况。事实上，CSR对矩阵的稀疏性没有要求，是一种适用于任何矩阵的通用存储方法，在诸稀疏矩阵的存储方式中也不见得是最高效的。
Example 1: 图1是来自Parallel Programming for FPGAs书中的举例，a)和b)分别展示了矩阵M的二维数组表示和CSR表示。CSR作为一种存储方式，也作为一种数据结构，由三个数组构成（所有的index默认从0开始）：
values：存储M中的非零元素的值；
columnIndex：存储非零元素在二维数组中水平方向的位置，即columnIndex[i]记录了values[i]元素的列序号；
rowPtr：该数组的长度为矩阵的行数加1，其中rowPtr[0]=0，其后的rowPtr[i]记录了矩阵前i行（即序号0至i-1的行）中包含非零元素的数量。
Reference: Kastner, R., Matai, J., & Neuendorffer, S. (2018). Parallel Programming for FPGAs. ArXiv, abs/1805.03648.
我们可以来模拟一下矩阵M（CSR表示）和向量的乘法的计算过程：
首先要用M的第一行向量和x相乘，由rowPtr[1]-rowPtr[0]=2知第一行共有2个非零元素；
从columnIndex中找到这两个元素所处的列，即columnIndex[0]=0和columnIndex[1]=1，这就意味着只需要考虑和分别与该行前两列元素的乘积（因为其余列元素都是0）；
再从values中找到这两个元素的值values[0]=3和values[1]=4，分别计算values[0]×和values[1]×并将两个乘积加起来即得到的值。
后面的计算都是先从rowPtr[i+1]-rowPtr[i]=k开始，然后查找columnIndex[rowPtr[i]]至columnIndex[rowPtr[i+k-1]]这k个列序号值，再然后根据列序号值在values中查找对应位置的k个元素值，分别与至
相乘后将k个乘积相加得到的值。
在程序实现中，上述操作需要2个嵌套循环，外层循环访问每一行，内层循环访问当前行中每一列元素。具体的实现代码从略。
尽管上文并没有刻意给出一个十分“稀疏”矩阵的例子，但也容易分析知：相比于二维数组，采用CSR存储方式避免了大量乘零冗余运算，的确能提高SpMV的实现效率。
由SpV切入的SpMSpV计算方法
在前一部分中，我们讨论了基于CSR的矩阵向量乘法，其核心思想是从矩阵的压缩着手，想办法避免由矩阵中的零元素造成的冗余运算。由（前一部分）Mx的计算过程可以看出，矩阵要先“主动宣布”自己在当前行有哪些列上的元素是非零的，然后向量再“被动派出”与这些非零元素对应的分量（“对应”是指非零元素的列序号与分量的行序号相等）与之相乘。当向量也十分稀疏时，它所“派出”的分量也有很大概率为零，这同样会导致相当多的冗余运算。既然由SpM切入可以使SpMV的计算更加高效，那么由SpV切入能否改善SpMSpV的计算效率呢？
Example 2：
由上例的表达式可以看出，对于仅含两个非零分量和的稀疏向量x，想要计算出M和x的乘积y，需要进行的乘法运算仅有以及，也就是说，参与Mx计算的M中的元素仅有列序号为r和s（和、的行序号对应相等）的两列。这样一来，我们只需要根据向量非零分量的行序号，选出矩阵对应列上的全部非零元素，例如我们假设列r的中仅有和这两个非零元素，那么只需要将二者分别与相乘，将两个乘积分别暂时保存在两个列表中，这两个列表分别存储所有形如和的结果，这是为了方便最后计算和（直接求对应列表中全部元素的和即得）。可以看到，此时SpV成为了“主动”的一方。
从上面所描述的计算过程可以看出，这里我们需要一种存储（或表示）方法，使得按列获取矩阵的非零元素更为便捷。第一部分中提到的CSR是按行优先访问非零元素的，它对行进行压缩；我们这里需要按列优先访问非零元素，那当然也可以对列进行压缩，这便是与CSR对应的CSC（Compressed Sparse Column，压缩列存储）。CSC与CSR的原理和结构极为类似，只不过是将“行”和“列”互换了而已，因此这里就不再详述了。
我们再来考虑一个更为简化的问题：将M和x中的非零元素都置为1，计算Mx。具体如下：
Example 3: Mx=y
如果我们只关注结果向量有哪些分量（行）非零，那么就根据向量x的第2、4两个分量非零，去矩阵M中查看第2、4两列，发现两列中的非零元素分布于第3、5、6、8四行（不重复列举），于是结果向量中第3、5、6、8四个分量就是非零的。如果要求出这四个分量具体的值，也很简单，第6行第2、4两列都是1，而其余没行都只有一列是1，所以结果向量第6个分量是2，第3、5、8三个分量都是1。可结合上面右边的图示来思考这一计算过程，这张图在后文会再次出现。
两种层同步并行BFS算法
层同步并行BFS有两种算法，分别是自顶向下（top-down）和自底向上（bottom-up）方法。通俗来讲，我们熟知的BFS树扩展方向是自顶向下，也就是parent-to-child的；而还有一种child-to-parent的BFS树扩展方向，就是自底向上。本部分主要介绍这两种算法的操作过程，以及引入自底向上方法的原因。
需要声明的一点是，本文并不对什么是“层同步并行”做解释（笔者本人也不懂这方面的知识），这里只需要了解这两种算法的执行过程，看懂其中的抽象操作即可（对这两种算法的分析也不需要有关“层同步并行”的知识）。下文所写多为简要概括性的内容，如有兴趣了解更多，可以参考Beamer, S., Asanović, K., & Patterson, D. (2013). Direction-optimizing breadth-first search. Scientific Programming, 21(3-4).这篇论文。
层同步并行BFS算法（Level Synchronized BFS）是最常用的并行图遍历方法：
1. 从任意随机顶点开始，作为当前层活跃顶点（Current Frontier）；
2. 由Current Frontier向未访问的邻接顶点（Neighbor）扩展；
3. 新遍历的顶点集合成为下一层遍历迭代的活跃顶点（Next Frontier）。
图2
Reference: https://www.ccita.net/wp-content/uploads/dlm_uploads/2021/01/4.2.4-面向高通量计算机的图算法优化技术-（中科院计算所-曹华伟）.pdf
Top-Down方法：
1. 扫描Current Frontier的所有顶点；
2. 检查其所有的Neighbor；
3. 添加未访问的Neighbor到Next Frontier。
Remark: Current Frontier活跃顶点数量较多时，遍历过程存在大量的冗余检查；并行执行时存在更新冲突。
Bottom-Up方法：
1. 扫描未访问的顶点；
2. 检查其Neighbor是否在Current Frontier，若有，则终止并检查下一顶点；
3. 将该顶点添加到Next Frontier。
Remark: 从未访问的顶点开始遍历，发现Neighbor位于Current Frontier即跳出循环。
Pseudocode Reference: 深入浅出BFS(1)_delibk的博客-CSDN博客
C++ Implementation: 深入浅出BFS(2)：自顶向下(Top-down method) & 自底向上(Bottom-up method)_自底向上bfs_delibk的博客-CSDN博客
BFS与矩阵向量乘法的联系
对于笔者这样的初学者来说，能够发现两个看似毫不相干概念之间的某种内在联系，无疑是一件令人兴奋的事情。让我们回头看看前文的Example 3，例子中给出的矩阵M实际上正是下面这个（有向无权）图G的邻接矩阵（Adjacency Matrix）：
为方便表示，设置该图的顶点表为： .
构造这个邻接矩阵的规则是：
注意，这里列下标代表的是有向边的起点，行下标代表的是有向边的终点。至于为什么要这么设定（通常都应该是行代表起点、列代表终点的，这里反过来了），看到后面自然就会明白。
回到Example 3，我们给向量x变个形，变形规则是：
经这样的变形后，Example 3中的矩阵向量乘法运算式就变成了下面这样：
Example 4: Mx‘=y‘
为什么要对向量做这么一个变形呢？有何用意？我们回顾一下第三部分所讲的层同步并行BFS算法：如果在图G中选择从顶点A开始层同步并行BFS过程，那么下一层遍历迭代的顶点就是B和D；接下来，以B和D为Current Frontier，进行下一轮扩展，其未访问的Neighbor有C、E、F、H（其中E、F是B的Neighbor，C、F、H是D的Neighbor，二者共有的Neighbor是F），于是C、E、F、H四个顶点成为Next Frontier。再看看Example 4中结果向量y'的各分量与一旁标注出来的顶点之间的对应关系，对照着想想前面层同步并行BFS树扩展的过程，是否就有所发现了？
图G以及Example 3、Example 4都是来源于GraphLily: Accelerating Graph Linear Algebra on HBM-Equipped FPGAs这篇论文中的举例，如图3所示。
Reference: Hu, Y., Du, Y., Ustun, E., & Zhang, Z. (2021). GraphLily: Accelerating Graph Linear Algebra on HBM-Equipped FPGAs. 2021 IEEE/ACM International Conference On Computer Aided Design (ICCAD), (1–9).
这张图已经把层同步并行BFS和矩阵向量乘法之间的联系表现得足够清晰，前文的分析也已经有足够的提示性，因此，为了省时省力，笔者就不在这里给出这种联系的具体的数学语言表达了。
经进一步分析我们还可以发现，层同步并行BFS的两种方法（top-down和bottom-up）还分别与SpMSpV和SpMV具有对应关系。图3(a)右边表现的就是第一部分所讲的按行优先遍历的SpMV计算过程，左边就是这个SpMV所对应的BFS过程，可以看出是按bottom-up方法进行的（其实看右边矩阵中各行的红色箭头所指的方向，可以看作是“由‘行’指向‘列’”，而“行”代表的是终点dst vertex，“列”代表的是起点src vertex，“由终点指向起点”，正是child-to-parent的方向，也就判断出是按bottom-up方法遍历的了；后面对top-down方法也可以像这样去判断）；而图3(b)右边表现的是按列优先遍历的SpMSpV计算过程，左边就是与之对应的按top-down方法进行的BFS过程。这样对照着看下来，本文标题所说的“两种层同步并行BFS算法与SpMV和SpMSpV之间的对应关系”，便十分明晰了。
GraphLily: Accelerating Graph Linear Algebra on HBM-Equipped FPGAs中还提到：
BFS typically applies SpMSpV at the first few iterations when the frontier is small (i.e., the input vector has a high sparsity), and switches to SpMV as the frontier grows (i.e., the input vector becomes denser).
BFS通常在前几次迭代中当Frontier的规模较小（即输入向量较为稀疏）时应用SpMSpV算法，并在当Frontier的规模增长（即输入向量变得更为稠密）时切换到SpMV算法。
至此，我们知道了不仅两种层同步并行BFS算法与SpMV和SpMSpV之间存在着两个对应关系，它们之间还能随着遍历迭代中由量变引起的质变而相互转化。我们可以将四者之间的联系表现在下面的这张图中，并以此作为本文的结尾：
Author: @Deng ZY文章来源：https://www.toymoban.com/news/detail-772660.html
Time: Mar 28, 2023文章来源地址https://www.toymoban.com/news/detail-772660.html
到了这里，关于【BFS应用】两种层同步并行BFS算法与SpMV和SpMSpV之间的对应关系的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
