Как можно выполнить разреженную индексацию матрицы для матрицы, хранящейся как «Compressed Sparse Row»?

Question

Моя большая разреженная симметричная матрица хранится как сжатая разреженная строка (CSR) с использованием Intel MKL.Для примера, давайте предположим, что моя симметричная разреженная матрица имеет вид 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, что я не считаю эффективным способом.

Может ли кто-нибудь сообщить мне эффективный или распространенный способ разреженной индексации матриц?

Answer 1

Хитрость заключается в поиске значений в нижнем треугольнике по выходному столбцу (который является их строкой).Вы можете сохранить индекс данных для каждой строки, поскольку вы посещаете записи в порядке столбцов по мере продвижения в порядке строк для вывода.

С типом экспозиции

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

мыhave

// 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.