我很有兴趣创建一个用于以块压缩稀疏行格式存储稀疏矩阵的类
这种存储方法包括将矩阵细分为大小为 sz*sz
的方块并将此块存储在向量 BA
中,在这里您可以找到关于link基本上矩阵使用4向量存储:
-
BA
包含以自上而下左右顺序存储的子矩阵(块)的元素(图片中大小为2x2
的第一个块是11,12,0,22
) -
AN
包含向量BA
的每个起始块的索引(在图示的情况下,块大小为2x2
,因此它包含1,5 ...
) -
AJ
包含块矩阵中块的列索引(图中较小的块) -
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个向量到密集矩阵的支持呢?例如如何打印出整个矩阵?
Edit : 我已经找到了用矢量AN的相对索引绘制图片中较小的"blocks matrix"的方法:
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 回答
回到稀疏矩阵:
其中:
BA_i
和AJ_i
是各向量的迭代器 .nBlock
保持ai_
给出的行数 .rBlock
和cBlock
是子矩阵sz*sz
的迭代器,称为"Block" .注意:
an_
仍未使用,您可以尝试更换BA_i whit .打印矩阵:
我不确定我是否正确编写了模板,无论如何算法应该可行;如果有问题,请告诉我 .
EDIT
你是对的,我的错;我认为这个问题重新组合了Matrix . 我使用findBlockIndex作为参考重写了方法 .
我希望对你有所帮助,最好的问候马可 .