求圈的CSR x CSR矩阵乘法

Question

我试图找出无向图中包含图中每个顶点u的顶点u的指定长度(K)的圈数。为此，我试图找出邻接矩阵的k‘’th幂。我从边列表中创建了图形的CSR表示。它工作得真快。但CSR乘法部分非常慢(输入大小为500 kx500k矩阵，需要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

我解决了我的问题，并想与那些也缺乏必要的“术语”的人分享，找出关于这个主题的大量资源。当你搜索“稀疏矩阵乘法”时，很难找到稀疏矩阵x稀疏矩阵。这被称为SpGEMM。关于这一过程，有许多信息丰富的论文。

我使用的算法的伪码：通用SpGEMM算法

我对算法做了一点修改，以产生CSR输出。与此相关的挑战似乎是对结果数组的分配，以容纳csr数组(值、index_array等)。解决这一问题的方法有多种，例如：

1. 分配与上限大小相同的数组。如果你的矩阵太大，这可能是个问题。如果您决定这样做，您可以查看：https://math.stackexchange.com/questions/1042096/bounds-of-sparse-matrix-multiplication。
2. 在为结果分配任何内存之前，可以执行乘法操作以确定结果中非零的数量。由于在这个空间中不存在内存写操作，所以结果很快就出来了。因此，结果数组所需的内存可以在“虚拟运行”之后分配。
3. 分配预先确定的金额，当分配一个新数组时，将内容复制到新的、较大的数组。

我实现了CPU (使用OpenMP)和GPU (使用CUDA)的功能。在OpenMP方法中，我使用了一个类似于我列出的选项3的方法。我对每一行的结果使用了不同的向量。比我加的结果向量还要多。向量方法可能比手动执行重新分配操作要慢，但它更容易，所以我选择了这种方法，而且速度足够快(测试矩阵有500 k行和500 k列，乘法操作大约需要1.3秒，在我的测试机器上使用60个线程)。对于GPU方法，我使用了选项2。首先，我计算了所需的数量，然后实际操作发生。

编辑：这个方法也能找到“行走”而不是路径。所以可能会有重复的顶点。