首页 文章

如何将块压缩行转换为密集矩阵?

提问于
浏览
2

我很有兴趣创建一个用于以块压缩稀疏行格式存储稀疏矩阵的类
BCSR format
这种存储方法包括将矩阵细分为大小为 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 回答

  • 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_iAJ_i 是各向量的迭代器 .

    nBlock 保持 ai_ 给出的行数 .

    rBlockcBlock 是子矩阵 sz*sz 的迭代器,称为"Block" .

    注意: an_ 仍未使用,您可以尝试更换BA_i whit .

    打印矩阵:

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

    我不确定我是否正确编写了模板,无论如何算法应该可行;如果有问题,请告诉我 .


    EDIT

    在一个为节省时间和存储器而创建的类中有意义,以某种方式存储稀疏矩阵,而不是使用向量来重建整个矩阵?

    你是对的,我的错;我认为这个问题重新组合了Matrix . 我使用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);
    }
    

    我希望对你有所帮助,最好的问候马可 .

相关问题