how to convert block compressed row to dense matrix?

[drudox84]

I'm interesting to create a class for storing sparse matrix in Block Compressed Sparse Row format : https://ibb.co/ffrKDR (picture)

this method of storage consist to subdivide the matrix into square block of size sz*sz and stored this block in a vector BA , here you can find most information about link basically the matrix is stored using 4 vector :

- BA contains the elements of the submatrices (blocks) stored in top-down left right order (the first block in the picture of size 2x2 is 11,12,0,22)
- AN contains the indices of each starting block of the vector BA (in the pictur case the block size is 2x2 so it contains 1,5 ... )

- AJ contains the column index of blocks in the matrix of blocks (the smaller one in the picture)

- AI the row pointer vector , it store how many blocks there is in the i-th row ai[i+1]-a[i] = number of block in i-th row

I'm write the constructor for convert a matrix from dense format to BCRS format :

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; 

}

And the main function for test the class :

 # 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;
}

if you want try this code you have just to add the include headers in the class file .

Now back to the question .. I obtain the 4 vector as in the picture .. but what about backing from this 4 vector to the dense matrix ? for example how to print out the whole matrix ?

[drudox84]

I've figure out the way to plot the "blocks matrix" the smaller in the picture with relative index of vector 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 ;
             }
             
          }
    }

that when in the main I call :

 bbcsr3.printBlockMatrix();

Give me the right result :

    1 0 0 0 
    2 3 0 0 
    0 0 4 5 
    0 0 0 6 

Now just the whole matrix missing I think that I missed something in may mind .. but should be something easy ..