Saya tertarik untuk membuat kelas untuk menyimpan matriks sparse dalam format Block Compressed Sparse Row Format BCSR tautan pada dasarnya matriks disimpan menggunakan 4 vektor :

  • BA berisi elemen-elemen submatriks (blok) yang disimpan dalam urutan kiri kanan atas-bawah (blok pertama pada gambar ukuran 2x2 adalah 11,12,0,22)
  • AN berisi indeks dari setiap blok awal vektor BA (dalam kasus gambar ukuran blok adalah 2x2 sehingga mengandung 1,5 ... )
  • AJ berisi indeks kolom blok dalam matriks blok (yang lebih kecil pada gambar)

  • AI vektor penunjuk baris , ini menyimpan berapa banyak blok yang ada di baris ke-i ai[i+1]-a[i] = number of block in i-th row

Saya sedang menulis konstruktor untuk mengonversi matriks dari format padat ke format 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; 

}

Dan fungsi utama untuk menguji kelas :

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

Sekarang kembali ke pertanyaan.. Saya mendapatkan 4 vektor seperti pada gambar.. tapi bagaimana dengan membackup dari 4 vektor ini ke matriks padat? misalnya bagaimana cara mencetak seluruh matriks?

Edit : Saya telah menemukan cara untuk memplot "blok matriks" yang lebih kecil pada gambar dengan indeks relatif vektor 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 ;
         }

      }
}

Bahwa ketika di utama saya memanggil:

bbcsr3.printBlockMatrix();

Beri saya hasil yang benar:

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

Sekarang hanya seluruh matriks yang hilang Saya pikir saya melewatkan sesuatu dalam pikiran .. tetapi seharusnya sesuatu yang mudah tetapi saya tidak mengerti maksudnya .. ada ide?

2
Drudox lebowsky 25 November 2017, 16:42

1 menjawab

Jawaban Terbaik

bagaimana dengan backing dari 4 vektor ini ke matriks padat? misalnya bagaimana cara mencetak seluruh matriks?

Kembali ke matriks jarang:

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

Di mana: BA_i dan AJ_i adalah iterator dari masing-masing vektor.

nBlock menyimpan jumlah blok dalam baris yang diberikan oleh ai_.

rBlock dan cBlockadalah iterator dari sub-matriks sz*sz yang disebut "Blok".

Catatan: an_ tetap tidak digunakan, Anda dapat mencoba mengganti BA_i sedikit pun.

Cetak matriks:

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

Saya tidak yakin saya menulis template dengan benar, bagaimanapun algoritme harus bekerja; beri tahu saya jika ada masalah.


EDIT

masuk akal di kelas yang dibuat untuk menghemat waktu dan memori menyimpan matriks jarang dengan cara tertentu daripada menggunakan vektor untuk merekonstruksi seluruh matriks?

Anda benar, salahku; Saya pikir masalahnya adalah menyusun ulang Matrix. Saya menulis ulang metode menggunakan findBlockIndex sebagai referensi.

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

Saya berharap dapat membantu Anda, salam Marco.

1
Marco D.G. 27 November 2017, 08:36