如何将块压缩行转换为密集矩阵？
Posted
技术标签:
【中文标题】如何将块压缩行转换为密集矩阵？【英文标题】：how to convert block compressed row to dense matrix? 【发布时间】：2017-11-25 13:42:20 【问题描述】：我很高兴创建一个用于以块压缩稀疏行格式存储稀疏矩阵的类 这种存储方法包括将矩阵细分为大小为sz*sz 的方形块，并将该块存储在一个向量中@987654325 @ , 在这里你可以找到关于link 的大部分信息
基本上矩阵是使用 4 个向量存储的：
BA 包含子矩阵（块）的元素以自上而下从左到右的顺序存储（图片中大小为2x2 的第一个块为11,12,0,22）
AN 包含向量BA 的每个起始块的索引（在图片中，块大小为2x2，因此它包含1,5 ...）
AJ包含blocks矩阵中blocks的列索引（图中较小的那个）
AI行指针向量，它存储了i-th行ai[i+1]-a[i] = number of block in i-th row有多少块@
我正在编写将矩阵从密集格式转换为 BCRS 格式的构造函数：
template <typename data_type, std::size_t SZ = 2 >
class BCSRmatrix
public:
constexpr BCSRmatrix(std::initializer_list<std::vector<data_type>> dense );
auto constexpr validate_block(const std::vector<std::vector<data_type>>& dense,
std::size_t i, std::size_t j) const noexcept ;
auto constexpr insert_block(const std::vector<std::vector<data_type>>& dense,
std::size_t i, std::size_t j) noexcept ;
private:
std::size_t bn ;
std::size_t bSZ ;
std::size_t nnz ;
std::size_t denseRows ;
std::size_t denseCols ;
std::vector<data_type> ba_ ;
std::vector<std::size_t> an_ ;
std::vector<std::size_t> ai_ ;
std::vector<std::size_t> aj_ ;
std::size_t index =0 ;
;
template <typename T, std::size_t SZ>
constexpr BCSRmatrix<T,SZ>::BCSRmatrix(std::initializer_list<std::vector<T>> dense_ )
this->denseRows = dense_.size();
auto it = *(dense_.begin());
this->denseCols = it.size();
if( (denseRows*denseCols) % SZ != 0 )
throw InvalidSizeException("Error block size is not multiple of dense matrix size");
std::vector<std::vector<T>> dense(dense_);
bSZ = SZ*SZ ;
bn = denseRows*denseCols/(SZ*SZ) ;
ai_.resize(denseRows/SZ +1);
ai_[0] = 1;
for(std::size_t i = 0; i < dense.size() / SZ ; i++)
auto rowCount =0;
for(std::size_t j = 0; j < dense[i].size() / SZ ; j++)
if(validate_block(dense,i,j))
aj_.push_back(j+1);
insert_block(dense, i, j);
rowCount ++ ;
ai_[i+1] = ai_[i] + rowCount ;
printBCSR();
template <typename T,std::size_t SZ>
inline auto constexpr BCSRmatrix<T,SZ>::validate_block(const std::vector<std::vector<T>>& dense,
std::size_t i, std::size_t j) const noexcept
bool nonzero = false ;
for(std::size_t m = i * SZ ; m < SZ * (i + 1); ++m)
for(std::size_t n = j * SZ ; n < SZ * (j + 1); ++n)
if(dense[m][n] != 0) nonzero = true;
return nonzero ;
template <typename T,std::size_t SZ>
inline auto constexpr BCSRmatrix<T,SZ>::insert_block(const std::vector<std::vector<T>>& dense,
std::size_t i, std::size_t j) noexcept
//std::size_t value = index;
bool firstElem = true ;
for(std::size_t m = i * SZ ; m < SZ * (i + 1); ++m)
for(std::size_t n = j * SZ ; n < SZ * (j + 1); ++n)
if(firstElem)
an_.push_back(index+1);
firstElem = false ;
ba_.push_back(dense[m][n]);
index ++ ;
template <typename T, std::size_t SZ>
auto constexpr BCSRmatrix<T,SZ>::printBCSR() const noexcept
std::cout << "ba_ : " ;
for(auto &x : ba_ )
std::cout << x << ' ' ;
std::cout << std::endl;
std::cout << "an_ : " ;
for(auto &x : an_ )
std::cout << x << ' ' ;
std::cout << std::endl;
std::cout << "aj_ : " ;
for(auto &x : aj_ )
std::cout << x << ' ' ;
std::cout << std::endl;
std::cout << "ai_ : " ;
for(auto &x : ai_ )
std::cout << x << ' ' ;
std::cout << std::endl;
以及测试类的主要功能：
# include "BCSRmatrix.H"
using namespace std;
int main()
BCSRmatrix<int,2> bbcsr2 = 11,12,0,0,0,0,0,0 ,0,22,0,0,0,0,0,0 ,31,32,33,0,0,0,0,0,
41,42,43,44,0,0,0,0, 0,0,0,0,55,56,0,0,0,0,0,0,0,66,67,0,0,0,0,0,0,0,77,78,0,0,0,0,0,0,87,88;
BCSRmatrix<int,4> bbcsr3 = 11,12,0,0,0,0,0,0 ,0,22,0,0,0,0,0,0 ,31,32,33,0,0,0,0,0,
41,42,43,44,0,0,0,0, 0,0,0,0,55,56,0,0,0,0,0,0,0,66,67,0,0,0,0,0,0,0,77,78,0,0,0,0,0,0,87,88;
return 0;
现在回到问题..我得到了图片中的 4 向量..但是从这个 4 向量支持到密集矩阵呢？ 例如如何打印出整个矩阵？
编辑：我已经找到了用向量 AN 的相对索引在图片中绘制“块矩阵”的方法：
template <typename T,std::size_t SZ>
inline auto constexpr BCSRmatrix<T,SZ>::printBlockMatrix() const noexcept
for(auto i=0 ; i < denseRows / SZ ; i++)
for(auto j=1 ; j <= denseCols / SZ ; j++)
std::cout << findBlockIndex(i,j) << ' ' ;
std::cout << std::endl;
template <typename T, std::size_t SZ>
auto constexpr BCSRmatrix<T,SZ>::findBlockIndex(const std::size_t r, const std::size_t c) const noexcept
for(auto j= ai_.at(r) ; j < ai_.at(r+1) ; j++ )
if( aj_.at(j-1) == c )
return j ;
当我主要打电话时：
bbcsr3.printBlockMatrix();
给我正确的结果：
1 0 0 0
2 3 0 0
0 0 4 5
0 0 0 6
现在只是缺少整个矩阵，我想我可能会错过一些东西..但应该很简单，但我没有明白这一点..有什么想法吗？
【问题讨论】：
【参考方案1】：从这 4 个向量支持到密集矩阵怎么样？例如如何打印出整个矩阵？
回到稀疏矩阵：
template <typename T, std::size_t SZ>
auto constexpr BCSRmatrix<T,SZ>::recomposeMatrix() const noexcept
std::vector<std::vector<data_type>> sparseMat(denseRows, std::vector<data_type>(denseCols, 0));
auto BA_i = 0, AJ_i = 0;
//for each BCSR row
for(auto r = 0; r < denseRows/SZ; r++)
//for each Block in row
for(auto nBlock = 0; nBlock < ai_.at(r+1)-ai_.at(r); nBlock++)
//for each subMatrix (Block)
for(auto rBlock = 0; rBlock < SZ; rBlock++)
for(auto cBlock = 0; cBlock < SZ; cBlock++)
//insert value
sparseMat.at(rBlock + r*SZ).at(cBlock + (aj_.at(AJ_i)-1)*SZ) = ba_.at(BA_i);
++BA_i;
++AJ_i;
return sparseMat;
在哪里：
BA_i 和 AJ_i 是各自向量的迭代器。
nBlock 保留ai_ 给出的行中的块数。
rBlock 和cBlock 是称为“块”的子矩阵sz*sz 的迭代器。
注意：an_ 未使用，您可以尝试用它替换 BA_i。
打印矩阵：
std::vector<std::vector<int>> sparseMat = bbcsr2.recomposeMatrix();
for(auto i = 0; i < sparseMat.size(); i++)
for(auto j = 0; j < sparseMat.at(i).size(); j++)
std::cout<<sparseMat.at(i).at(j) << '\t';
std::cout << std::endl;
我不确定我是否正确编写了模板，无论如何算法应该可以工作；如果有问题请告诉我。
编辑
在一个为节省时间和内存存储稀疏矩阵而创建的类中是否有意义，而不是使用向量来重建整个矩阵？
你是对的，我的错；我认为问题是重组矩阵。 我使用 findBlockIndex 作为参考重写了方法。
template <typename T, std::size_t SZ>
auto constexpr BCSRmatrix<T,SZ>::printSparseMatrix() const noexcept
//for each BCSR row
for(auto i=0 ; i < denseRows / SZ ; i++)
//for each Block sub row.
for(auto rBlock = 0; rBlock < SZ; rBlock++)
//for each BCSR col.
for(auto j = 1; j <= denseCols / SZ; j++)
//for each Block sub col.
for(auto cBlock = 0; cBlock < SZ; cBlock++)
std::cout<< findValue(i, j, rBlock, cBlock) <<'\t';
std::cout << std::endl;
template <typename T, std::size_t SZ>
auto constexpr BCSRmatrix<T,SZ>::findValue(const std::size_t i, const std::size_t j, const std::size_t rBlock, const std::size_t cBlock) const noexcept
auto index = findBlockIndex(i,j);
if(index != 0)
return ba_.at(an_.at(index-1)-1 + cBlock + /* rBlock*2 */ rBlock*SZ);
希望对你有所帮助 最好的问候马可。
【讨论】：
谢谢@Marco D.G. ，你的方法的扩展性如何？我的意思是 .. 它可以工作，但是在一个为节省时间和内存存储稀疏矩阵而创建的类中是否有意义，而不是使用向量来重建整个矩阵？我知道这就是我所要求的！你的代码很棒！但是您如何看待名为findIndex(i,j) 的方法，如果它不为零则返回原始矩阵的值，否则返回 0？你不认为应该更好吗？如何打印出整个矩阵而不是将其存储在 vector以上是关于如何将块压缩行转换为密集矩阵？的主要内容，如果未能解决你的问题，请参考以下文章