2
\$\begingroup\$

Here my Block Compressed sparse Row storage class, I have some doubt about the constructor that use a std::vector<std::vector> for store the initializer list from which the matrix is construct, and also in the constructor from file, is there a less expansive way ?

It works fine but I'm wondering about the sense .. . In the method BCRS the matrix is split into br×bc submatrices (in my case br=bc), (called blocks), where br and bc are fixed integers. BCRS stores the non-zero blocks (submatrices with at least one non-zero element) the matrix storage result in this way :

  • Ba contain all the element of the block (top- bottom left right order)
  • An contain the index of starting block in the BA vector
  • Aj contain the block column in the block matrix see the picture
  • Ai contain the number of block per row Ai[i+1]-Ai[i] = number of block in i-th row here the matrix

this is the class that I wrote :

template <typename T, std::size_t S >
class BCRSmatrix ;


template <typename T, std::size_t S>
std::ostream& operator<<(std::ostream& os , const BCRSmatrix<T,S>& m ) noexcept ;


template <typename T, std::size_t S>
std::vector<T> operator*(const BCRSmatrix<T,S>& m, const std::vector<T>& x );


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 );  

     constexpr BCRSmatrix(const std::string& );  

     virtual ~BCRSmatrix() = default ; 

     auto constexpr print_block(const std::vector<std::vector<data_type>>& dense,
                                  std::size_t i, std::size_t j) const noexcept ; 

     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 ; 

     auto constexpr printBlockMatrix() const noexcept  ; 

     auto constexpr size1() const noexcept { return denseRows ;}

     auto constexpr size2() const noexcept { return denseCols ;} 

     auto constexpr printBlock(std::size_t i) const noexcept ;

     auto constexpr print() const noexcept ; 

  //   auto constexpr printMatrix() 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 ;

};

template <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 InvalidSizeException("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 ;
     }

}

template <typename T, std::size_t BS>
constexpr BCRSmatrix<T,BS>::BCRSmatrix(const std::string& fname)
{
    std::ifstream f(fname , std::ios::in);
    if(!f)
    {
       throw OpeningFileException("error opening file in constructor !");
    }
    else
    {
       std::vector<std::vector<T>> dense;
       std::string line, tmp;
       T elem = 0 ;
       std::vector<T> row;
       std::size_t i=0, j=0 ;     

       while(getline(f, line))
       {
          row.clear();  
          std::istringstream ss(line);
          if(i==0) 
          {
            while(ss >> elem)
            {
              row.push_back(elem);
              j++;
            }
          }
          else
          {
            while(ss >> elem) 
              row.push_back(elem);
          }
          dense.push_back(row);  
          i++;  
       }

       this->denseRows = i;
       this->denseCols = j;

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

    }


}

template <typename T,std::size_t BS>
inline auto constexpr BCRSmatrix<T,BS>::printBlockMatrix() const noexcept  
{

      for(auto i=0 ; i < denseRows / BS ; i++)
      {
        for(auto j=1 ; j <= denseCols / BS  ; j++)
        {
            std::cout << findBlockIndex(i,j) << ' ' ;  
        }
         std::cout << std::endl;   
      }
}

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>::print() const noexcept 
{      

    for(auto i=0 ; i < denseRows / BS ; i++)
    {
       for(auto rBlock = 0; rBlock < BS; rBlock++)
        {

            for(auto j = 1; j <= denseCols / BS; j++){

                for(auto cBlock = 0; cBlock < BS; cBlock++){
                    std::cout<< findValue(i, j, rBlock, cBlock) <<'\t';
                }
            }
            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(auto rBlock = 0; rBlock < BS; rBlock++)
        {

            for(auto j = 1; j <= m.denseCols / BS; j++)
            {

                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 InvalidSizeException(mess.c_str());
      }
      else
      {
            auto brows = m.denseRows/BS ;  
            auto bnze  = m.an_.size()   ;

            auto z=0;

            for(auto b=0 ; b < brows ; b++)
            {     
               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;      

}

here the main function

# include "BCRSmatrix.H"

using namespace std;

int main(){

   BCRSmatrix<int,2> 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}};

 bbcsr3.printBlockMatrix();
cout << bbcsr3 ;
std::vector<int> v1 = {3,4,0,1,6,8,1,19};

std::vector<int> v2 = bbcsr3 *v1 ;

  for(auto& x : v2)
     cout << x << ' ' ;
  cout << endl;   

EDIT: I've forgot the MatrixException file :

  class MatrixException : public std::exception {

   public:   

    explicit MatrixException(std::string mess) : std::exception{} , message{mess}
           {} 

    virtual ~MatrixException() noexcept = default ;

    inline virtual const char* what() const noexcept { return this-> message.c_str() ;}  

   protected: 

     std::string message ;


};



class InvalidSizeException : public MatrixException { 

      public:

         InvalidSizeException(const std::string &mess) : MatrixException{mess}
                  {}

};

class InvalidCoordinateException : public MatrixException {

      public:
       InvalidCoordinateException(const std::string &mess) : MatrixException{mess}     
                  {}
};

class OpeningFileException : public MatrixException {

      public :
       OpeningFileException(const std::string &mess) : MatrixException{mess}
               {}
};
\$\endgroup\$
  • 2
    \$\begingroup\$ You have debugging code mixed with your main code (like prints in the constructor), commented out dead code, and misleading indentation. You may want to put a small polish pass to your code before gettig it reviewed. \$\endgroup\$ – Frank Nov 29 '17 at 14:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.