이 alghorithm을 찾을 수없는 이유는 무엇입니까?

Question

나는 BCRSmatrix 클래스 (블록 압축 열 저장)를 구현하고있어가, 나는 컴퓨팅 행렬 벡터 제품 여기 the algorithm 발견했습니다, 그래서 정확히 alghorithm하지만 뭔가 이상한 다음 operator*(std::vector<type>,BCRSmatrix<type,size>)을 정의 발생합니다 이 행렬을 정의 할 수 있습니다 :

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 

이 벡터 :

3,4,0,1,6,8,1,19 

행렬의 곱 (블록 사이즈 2 × 2)와 줘 :

,913,210

(블록 사이즈 × 4)와 행렬의 곱하면 줘 :

81 88 221 335 778 595 1,559 1,759

I 돈 :

81 0 221 335 778 595 1559 1759 

옥타브 나에게 적절한 결과를 제공 왜 그런지 이해하지 못해! 디버거로 시도했지만 이상한 것을 찾지 못했습니다!

최소 컴파일 가능한 코드 : 수정 버그를 수정했습니다! 죄송합니다, 당신이

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

} 

여기에 주요 기능을 @Marco 감사 :

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

Answer 1

당신에 대한 첫 번째 내부 y.at(b) = 0;을 제거해야합니다!