I'm implementing the BCRSmatrix class (block compressed-row storage), I've found here the algorithm for compute matrix vector product, so I defined the operator*(std::vector<type>,BCRSmatrix<type,size>) following exactly this alghorithm but something strange happens :
lets define this matrices :
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
and this vector:
3,4,0,1,6,8,1,19
the product of the matrix (with block size 2x2) give me :
81 0 0 0 778 595 1559 1759
the product of the matrix (with block size 4x4) give me :
81 0 221 335 778 595 1559 1759
octave give me the right result :
81 88 221 335 778 595 1559 1759
I don't understand why ! I try with debugger but I didn't find nothing strange!
here the minimal compilable code : EDIT I've fix the bug ! sorry , thank you @Marco
# include <iosfwd>
# include <vector>
# include <string>
# include <initializer_list>
# include <sstream>
# include <fstream>
# include <algorithm>
# include <iomanip>
template <typename data_type, std::size_t BS>
class BCRSmatrix {
template <typename T, std::size_t S>
friend std::ostream& operator<<(std::ostream& os , const BCRSmatrix<T,S>& m ) noexcept ;
template <typename T, std::size_t S>
friend std::vector<T> operator*(const BCRSmatrix<T,S>& m, const std::vector<T>& x );
public:
constexpr BCRSmatrix(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 ;
auto constexpr printBCRS() const noexcept ;
private:
std::size_t bn ;
std::size_t bBS ;
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 ;
auto constexpr findBlockIndex(const std::size_t r, const std::size_t c) const noexcept ;
auto constexpr findValue(
const std::size_t i, const std::size_t j,
const std::size_t rBlock, const std::size_t cBlock
) const noexcept ;
};
emplate <typename T, std::size_t BS>
constexpr BCRSmatrix<T,BS>::BCRSmatrix(std::initializer_list<std::vector<T>> dense_ )
{
this->denseRows = dense_.size();
auto it = *(dense_.begin());
this->denseCols = it.size();
if( (denseRows*denseCols) % BS != 0 )
{
throw std::runtime_error("Error block size is not multiple of dense matrix size");
}
std::vector<std::vector<T>> dense(dense_);
bBS = BS*BS ;
bn = denseRows*denseCols/(BS*BS) ;
ai_.resize(denseRows/BS +1);
ai_[0] = 1;
for(std::size_t i = 0; i < dense.size() / BS ; i++)
{
auto rowCount =0;
for(std::size_t j = 0; j < dense[i].size() / BS ; j++)
{
if(validate_block(dense,i,j))
{
aj_.push_back(j+1);
insert_block(dense, i, j);
rowCount ++ ;
}
}
ai_[i+1] = ai_[i] + rowCount ;
}
printBCRS();
}
template <typename T,std::size_t BS>
inline auto constexpr BCRSmatrix<T,BS>::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 * BS ; m < BS * (i + 1); ++m)
{
for(std::size_t n = j * BS ; n < BS * (j + 1); ++n)
{
if(dense[m][n] != 0) nonzero = true;
}
}
return nonzero ;
}
template <typename T,std::size_t BS>
inline auto constexpr BCRSmatrix<T,BS>::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 * BS ; m < BS * (i + 1); ++m)
{
for(std::size_t n = j * BS ; n < BS * (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 BS>
auto constexpr BCRSmatrix<T,BS>::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 ;
}
}
}
template <typename T, std::size_t BS>
auto constexpr BCRSmatrix<T,BS>::printBCRS() 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;
}
template <typename T, std::size_t BS>
auto constexpr BCRSmatrix<T,BS>::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*BS);
}
template <typename T, std::size_t BS>
std::ostream& operator<<(std::ostream& os , const BCRSmatrix<T,BS>& m ) noexcept
{
for(auto i=0 ; i < m.denseRows / BS ; i++)
{
//for each Block sub row.
for(auto rBlock = 0; rBlock < BS; rBlock++)
{
//for each BCSR col.
for(auto j = 1; j <= m.denseCols / BS; j++)
{
//for each Block sub col.
for(auto cBlock = 0; cBlock < BS; cBlock++)
{
os << m.findValue(i, j, rBlock, cBlock) <<'\t';
}
}
os << std::endl;
}
}
return os;
}
template <typename T, std::size_t BS>
std::vector<T> operator*(const BCRSmatrix<T,BS>& m, const std::vector<T>& x )
{
std::vector<T> y(x.size());
if(m.size1() != x.size())
{
std::string to = "x" ;
std::string mess = "Error occured in operator* attempt to perfor productor between op1: "
+ std::to_string(m.size1()) + to + std::to_string(m.size2()) +
" and op2: " + std::to_string(x.size());
throw std::runtime_error(mess.c_str());
}
else
{
auto brows = m.denseRows/BS ;
auto bnze = m.an_.size() ;
auto z=0;
for(auto b=0 ; b < brows ; b++)
{
// y.at(b) = 0; <-- bug
for(auto j= m.ai_.at(b) ; j <= m.ai_.at(b+1)-1; j++ )
{
for(auto k=0 ; k < BS ; k++ )
{
for(auto t=0 ; t < BS ; t++)
{
y.at(BS*b+k) += m.ba_.at(z) * x.at(BS*(m.aj_.at(j-1)-1)+t) ;
z++ ;
}
}
}
}
}
return y;
}
and here the main function :
using namespace std;
int main(){
BCRSmatrix<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}};
std::vector<int> v1 = {3,4,0,1,6,8,1,19};
std::vector<int> v2 = bbcsr2 *v1 ;
for(auto& x : v2)
cout << x << ' ' ;
cout << endl;
return 0;
}
