首页 文章

如何对存储为“压缩稀疏行”的矩阵进行稀疏矩阵索引?

提问于
浏览
1

我使用 Intel MKL 将我的大型稀疏对称矩阵存储为压缩稀疏行(CSR) . 为了举例,让我们假设我的对称稀疏矩阵是 5x5

A =
    1    -1     0    -3     0
   -1     5     0     0     0
    0     0     4     6     4
   -3     0     6     7     0
    0     0     4     0    -5

values   = {1,  -1,  -3,   5,   4,   6,   4,   7,  -5}; // symmetric sparse matrix
columns  = {0,   1,   3,   1,   2,   3,   4,   3,   4}; // zero-based
rowIndex = {0,   3,   4,   7,   8,   9}; // zero-based

我试图在给定行和列的情况下找到 A 的子矩阵,例如 A(1:3, 2:4)

A(1:3,2:4) =
   0     0     0
   4     6     4
   6     7     0

values   = {4,   6,   4,   6,   7}; // General sparse matrix (sub-matrix is not necessarily symmetric)
columns  = {0,   1,   2,   0,   1}; // zero-based
rowIndex = {0,   0,   3,   5}; // zero-based

我很高兴知道矩阵索引是如何完成的 . 我能想到的一种方法是将 CSR 转换为坐标格式 COO 并应用矩阵索引,然后将其转换回 CSR ,我认为这不是一种有效的方法 .

有人能让我知道稀疏矩阵索引的有效或常用方法吗?

1 回答

  • 1

    诀窍是通过输出列(它们的行)在下三角形中查找值 . 您可以为每行的数据保留一个索引,因为您按行顺序访问输入的行顺序时的条目 .

    与博览会类型

    struct CSR {  // sometimes implicitly symmetric
      std::vector<...> vals;
      std::vector<int> cols,rowStart;
    };
    

    我们有

    // Return the [r0,r1) by [c0,c1) submatrix, never
    // using any symmetry it might have.
    CSR submatrix(const CSR &sym,int r0,int r1,int c0,int c1) {
      const int m=r1-r0,n=c1-c0;
      std::vector<int> finger(sym.rowStart.begin()+c0,sym.rowStart.begin()+c1);
      CSR ret;
      ret.rowStart.reserve(m+1);
      ret.rowStart.push_back(0);
      for(int r=0,rs=r0;r<m;++r,++rs) {
        // (Strictly) lower triangle:
        for(int cs=c0,c=0;cs<rs;++cs,++c)
          for(int &f=finger[c],f1=sym.rowStart[cs+1];f<f1;++f) {
            const int cf=sym.cols[f];
            if(cf>rs) break;
            if(cf==rs) {
              ret.vals.push_back(sym.vals[f]);
              ret.cols.push_back(c);
            }
          }
        // Copy the relevant subsequence of the upper triangle:
        for(int f=sym.rowStart[rs],f1=sym.rowStart[rs+1];f<f1;++f) {
          const int c=sym.cols[f]-c0;
          if(c<0) continue;
          if(c>=n) break;
          ret.vals.push_back(sym.vals[f]);
          ret.cols.push_back(c);
        }
        ret.rowStart.push_back(ret.vals.size());
      }
      return ret;
    }
    

    对于大型矩阵,可以通过使用二进制搜索来优化上三角形循环以找到 f 的相关范围 .

相关问题