如何对存储为 "Compressed Sparse Row" 的矩阵进行稀疏矩阵索引?
How sparse matrix-indexing can be done for a matrix stored as "Compressed Sparse Row"?
我使用 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
,我认为这不是一种有效的方法.
有人可以告诉我一种有效或常见的稀疏矩阵索引方式吗?
诀窍是通过输出列(即它们的行)查找下三角中的值。您可以为每一行的数据保留一个索引,因为您在按行顺序进行输出时按列顺序访问条目。
与说明型
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
来优化。
我使用 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
,我认为这不是一种有效的方法.
有人可以告诉我一种有效或常见的稀疏矩阵索引方式吗?
诀窍是通过输出列(即它们的行)查找下三角中的值。您可以为每一行的数据保留一个索引,因为您在按行顺序进行输出时按列顺序访问条目。
与说明型
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
来优化。