CSR x CSR 矩阵乘法找出周期

Question

我正在尝试找出具有指定长度 (k) 的无向图中的循环数，其中包含图中每个顶点 u 的顶点 u。为此，我试图找出邻接矩阵的第 k 次方。我从边缘列表创建了图形的 CSR 表示。它工作得非常快。但是 CSR x CSR 乘法部分真的很慢，（输入大小为 500k x 500k 矩阵似乎需要 50 分钟）。我很好奇更好的解决方案。由于这是邻接矩阵，是否有更有效的方法？或者有没有更好的 CSRxCSR 矩阵乘法可供我查看？我找不到任何 CSR X CSR 矩阵乘法示例作为算法或 C++ 实现。

    void multiply_matrix(std::vector<int> &adj, std::vector<int> &xadj, std::vector<int> &values, std::vector<int> &adj2, std::vector<int> &xadj2, std::vector<int> &values2, int size)
  {
          std::vector<int> result_adj;
          std::vector<int> result_xadj(size+1,0);
          std::vector<int> result_value(values.size(),0);
          for(int i = 0; i<size; i++)
          {
                  for(int j = 0; j<size; j++)
                  {
                          int result = 0;
                          int startIndex = xadj[i];
                          int endIndex = xadj[i+1];
                          for(int index = startIndex; index<endIndex; index++)
                          {
                                  int currentValRow = values[adj[index]];
                                  bool shouldContinue = false;
                                  for(int colIndex = xadj2[j]; colIndex<xadj2[j+1]; colIndex++)
                                  {
                                          if(adj[index] == adj2[colIndex])
                                          {
                                                  shouldContinue = true;
                                                  break;
                                          }
                                  }
                                  if(!shouldContinue)
                                          continue;
                                  int currentValCol = values2[adj2[index]];
                                  result += currentValCol*currentValRow;
                          }
                          if(result != 0)
                          {
                                  result_xadj[i+1]++;
                                  if(i+2 < result_xadj.size())
                                          result_xadj[i+2] = result_xadj[i+1];
                                  result_adj.push_back(j);
                                  result_value[j] = result;
                          }
                  }
          }
  }

Answer 1

我解决了我的问题，并想与那些同样缺乏所需“术语”的人分享，以找到有关该主题的大量资源。当您 google “稀疏矩阵乘法”时，很难找到稀疏矩阵 x 稀疏矩阵。这称为 SpGEMM。有很多关于该过程的信息性论文。

我使用的算法伪代码： General SpGEMM algorithm

我稍微修改了算法以生成 CSR 输出。挑战似乎是结果数组的分配来保存 csr 数组（值，index_array，等等）。有多种方法可用于解决该问题，例如：

分配与上限一样大的数组。如果您的矩阵太大，这可能是个问题。如果您决定走这条路，您可以查看：https://math.stackexchange.com/questions/1042096/bounds-of-sparse-matrix-multiplication.
在为结果分配任何内存之前，可以进行乘法运算以确定结果中非零的数量。由于此 space 中不存在内存写入操作，因此结果非常快。所以结果数组所需的内存可以在这个“dummy 运行”之后分配。
分配一个预先确定的数量，当它不够时分配一个新数组并将内容复制到新的、更大的数组。

我为 CPU（使用 OpenMP）和 GPU（使用 CUDA）实现了该功能。在 OpenMP 方法中，我使用了一种类似于我列出的选项 3 的方法。我为每一行的结果使用了单独的向量。比我添加结果向量。向量方法可能比手动进行重新分配操作慢，但它更容易，所以我选择了这种方式并且它足够快（测试矩阵有 500k 行和 500k 列，乘法操作大约需要 1.3 秒，使用我的 60 个线程试机）。对于 GPU 方法，我使用了选项 2。首先我计算了所需的数量，然后才开始实际操作。

编辑： 此外，此方法找出“步行”而不是路径。所以可能会有重复的顶点。

CSR x CSR 矩阵乘法找出周期

CSR x CSR Matrix Multiplication for Finding Out Cycles

graph

matrix

sparse-matrix

matrix-multiplication

adjacency-matrix