2017-11-21 14 views
1

안녕하세요, 저는 모두 수정 된 압축 된 드문 드문 행에 저장된 희소 행렬을 구현하고 있습니다! 생성자가 잘 작동하지만, 나는 그것을 검증 할 수 있지만 operator +는 이상한 행동을 보입니다. 만일 내가 0이 아닌 값을 가지고 있다면 합계는 올바른 결과를 계산하지 못합니다. 수정 된 압축 스파 스 열 방법은 코드를 작업하는 후속이다 here 내 미님를 설명 :연산자 + 스파 스 매트릭스 : 나는 어떤 버그도 찾지 못했습니다.

# include <initializer_list> 
# include <vector> 
# include <iosfwd> 
# include <string> 
# include <cstdlib> 
# include <cassert> 
# include <iomanip> 
# include <cmath> 
# include <set> 
# include <fstream> 
# include <algorithm> 
# include <exception> 

template <typename data_type> 
class MCSRmatrix;    

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


template <typename T> 
MCSRmatrix<T> operator+(const MCSRmatrix<T>& m1, const MCSRmatrix<T>& m2) ; 

template <typename data_type> 
class MCSRmatrix { 

    public: 
     using itype = std::size_t ; 

     template <typename T> 
     friend std::ostream& operator<<(std::ostream& os ,const MCSRmatrix<T>& m) noexcept ; 

     template <typename T> 
     friend MCSRmatrix<T> operator+(const MCSRmatrix<T>& m1, const MCSRmatrix<T>& m2) ; 

    public: 

     constexpr MCSRmatrix(std::initializer_list<std::initializer_list<data_type>> rows); 

     constexpr MCSRmatrix(const std::size_t&) noexcept; 

     const data_type& operator()(const itype r , const itype c) const noexcept ; 

     data_type operator()(const itype r , const itype c) noexcept ; 

     auto constexpr printMCSR() const noexcept ; 

    private: 

    std::vector<data_type> aa_ ; // vector of value 
    std::vector<itype>  ja_ ; // pointer vector 

    int dim ; 

    std::size_t constexpr findIndex(const itype row , const itype col) const noexcept ; 
}; 

template <typename T> 
constexpr MCSRmatrix<T>::MCSRmatrix(std::initializer_list<std::initializer_list<T>> rows) 
{ 
     this->dim = rows.size(); 
     auto _rows = *(rows.begin()); 

     aa_.resize(dim+1); 
     ja_.resize(dim+1); 

     if(dim != _rows.size()) 
     { 
      throw std::runtime_error("Error in costructor! MCSR format require square matrix!"); 
     } 

     itype w = 0 ; 
     ja_.at(w) = dim+2 ; 
     for(auto ii = rows.begin(), i=1; ii != rows.end() ; ++ii, i++) 
     { 
      for(auto ij = ii->begin(), j=1, elemCount = 0 ; ij != ii->end() ; ++ij, j++) 
      { 
       if(i==j) 
       aa_[i-1] = *ij ; 
       else if(i != j && *ij != 0) 
       { 
       ja_.push_back(j); 
       aa_.push_back(*ij); 
       elemCount++ ; 
       } 
       ja_[i] = ja_[i-1] + elemCount;   
      } 
     }  
     printMCSR(); 
} 





template <typename T> 
constexpr MCSRmatrix<T>::MCSRmatrix(const std::size_t& n) noexcept 
{ 
     this->dim = n; 
     aa_.resize(dim+1); 
     ja_.resize(dim+1); 

     for(std::size_t i = 0; i < aa_.size()-1 ; i++) 
     aa_.at(i) = 1 ;  

     for(std::size_t i = 0; i < ja_.size() ; i++) 
     ja_.at(i) = aa_.size()+1 ; 


}  


template <typename T> 
std::size_t constexpr MCSRmatrix<T>::findIndex(const itype row , const itype col) const noexcept 
{  
    assert(row > 0 && row <= dim && col > 0 && col <= dim); 
    if(row == col) 
    { 
     return row-1; 
    } 
    int i = -1; 

     for(int i = ja_.at(row-1)-1 ; i < ja_.at(row)-1 ; i++) 
     { 
      if(ja_.at(i) == col) 
      { 
       return i ; 
      } 
     } 
     return -1; 
} 


template <typename T> 
inline auto constexpr MCSRmatrix<T>::printMCSR() const noexcept 
{ 
     for(auto& x : aa_) 
      std::cout << x << ' ' ; 
     std::cout << std::endl; 

     for(auto& x : ja_) 
      std::cout << x << ' ' ; 
     std::cout << std::endl; 
} 

template <typename T> 
const T& MCSRmatrix<T>::operator()(const itype r , const itype c) const noexcept 
{  
    auto i = findIndex(r,c); 
    if(i != -1 && i < aa_.size()) 
    {  
      return aa_.at(i) ; 
    } 
    else 
    { 
      return aa_.at(dim) ; 
    }  
} 


template <typename T> 
T MCSRmatrix<T>::operator()(const itype r , const itype c) noexcept 
{  
    auto i = findIndex(r,c); 
    if(i != -1 && i < aa_.size()) 
    {  
      return aa_.at(i) ; 
    } 
    else 
    { 
      return aa_.at(dim) ; 
    }  
} 

// non member operator 

template <typename T> 
std::ostream& operator<<(std::ostream& os ,const MCSRmatrix<T>& m) noexcept 
{ 
    for(std::size_t i=1 ; i <= m.dim ; i++) 
    {  
     for(std::size_t j=1 ; j <= m.dim ; j++) 
      os << std::setw(8) << m(i,j) << ' ' ; 
     os << std::endl; 
    } 
    return os; 
} 


// perform sum by 2 Mod Comp Sparse Row matrices 
template <typename T> 
MCSRmatrix<T> operator+(const MCSRmatrix<T>& m1, const MCSRmatrix<T>& m2) 
{ 
     if(m1.dim != m2.dim) 
     { 
      throw std::runtime_error("Matrixs dimension does not match! Error in operator +"); 
     } 
     else 
     { 
     MCSRmatrix<T> res(m1.dim); 



     for(auto i=0; i < res.dim ; i++) 
      res.aa_.at(i) = m1.aa_.at(i) + m2.aa_.at(i) ; 

     res.ja_.at(0) = res.dim+2; 

     std::set<unsigned int> ctrl; 


     int n1=0, n2=0, j1=0 , j2 =0; 
     for(auto i=0 ; i < res.dim ; i++) 
     { 

      ctrl.clear(); 

      n1 = m1.ja_.at(i+1)- m1.ja_.at(i) ; 
      n2 = m2.ja_.at(i+1)- m2.ja_.at(i) ; 


      j1=0 , j2=0 ; 

      auto sum1 = 0. , sum2 = 0. , sum=0.; 


      for(auto j = 0; j < std::max(n1,n2) ; j++) 
      { 
       if(n1 && n2) 
       { 
       if(m1.ja_.at(m1.ja_.at(i)-1 + j1) == m2.ja_.at(m2.ja_.at(i)-1 + j2)) 
       { 
        ctrl.insert(m1.ja_.at(m1.ja_.at(i)-1 + j1)); 

        sum = m1.aa_.at(m1.ja_.at(i)-1 + j1) + m2.aa_.at(m2.ja_.at(i)-1 + j2) ; 


       } 
       else if(m1.ja_.at(m1.dim+1 + j1) != m2.ja_.at(m2.dim+1 + j2)) 
       { 
        ctrl.insert(m1.ja_.at(m1.ja_.at(i)-1 + j1)); 
        ctrl.insert(m2.ja_.at(m1.ja_.at(i)-1 + j2)); 

        sum1 = m1.aa_.at(m1.ja_.at(i)-1 + j1);    
        sum2 = m2.aa_.at(m2.ja_.at(i)-1 + j1);  
       } 
       }   
       else if(n1) 
       { 
        ctrl.insert(m1.ja_.at(m1.ja_.at(i)-1 + j1)); 
        sum1 = m1.aa_.at(m1.ja_.at(i)-1 + j1);    
       } 
       else if(n2) 
       { 
        ctrl.insert(m2.ja_.at(m2.ja_.at(i)-1 + j2)); 
        sum2 = m2.aa_.at(m2.ja_.at(i)-1 + j2);  
       } 
       if(sum1) 
       { 
        res.aa_.push_back(sum1); 
        res.ja_.push_back(m1.ja_.at(m1.ja_.at(i)-1 + j1)) ; 
       } 
       if(sum2) 
       { 
        res.aa_.push_back(sum2); 
        res.ja_.push_back(m2.ja_.at(m2.ja_.at(i)-1 + j2)); 
       } 
       if(sum) 
       {  

        res.aa_.push_back(sum); 
        res.ja_.push_back(m1.ja_.at(m1.ja_.at(i)-1 + j1)); 
       } 

       if(j1 < n1) j1++ ; 
       if(j2 < n2) j2++ ; 
      } 


      res.ja_.at(i+1) = res.ja_.at(i) + ctrl.size() ; 

     } 
     return res ; 
     } 
} 

이 주요 프로그램 :

# include "ModifiedCSRmatrix.H" 

using namespace std ; 


int main() { 

    MCSRmatrix<int> m01 = {{0,1,0,0},{0,3,0,0},{0,0,0,3},{0,0,0,1}}; 
    MCSRmatrix<int> m02 = {{1,1,0,14},{0,22,0,3},{0,0,33,34},{0,0,1,44}}; 

    cout << endl; 
    MCSRmatrix<int> m03 = m01+m02 ; 

    cout << m03 << endl; 
    cout << endl; 
    m03.printMCSR(); 




    return 0; 
} 

이 특정 2 매트릭스 나에게 올바른 결과를 제공!

1  2  0  17 
0  25  0  3 
0  0  33  37 
0  0  1  45 

하지만 그렇게 m01 행렬 .. 과의 소자 (3,2)의 값을 변경하는 경우 행렬되었다 : 즉

:

MCSRmatrix<int> m01 = {{0,1,0,0},{0,3,0,0},{0,2,0,3},{0,0,0,1}}; 

코드 날이 잘못된 결과를 제공

1  2  0  0 
    0  25  0  14 
    0  2  33  0 
    0  0  0  45 

그러나 벡터 aa_ (첫 번째 matrix.dim이 ele 대각선 따라 소자의 멘트가 오프 DIAG 영이 아닌 요소들이있다 행렬)이이 결과로부터 다른 것으로 보인다 같습니다

aa_ = 1 25 33 45 0 2 2 14 2 3 3 1 1 
ja_ = 6 8 9 10 11 2 2 4 2 4 4 3 3 

I 주에이 라인을 가지고 다른 한 예 경우 :

제가

예 M1에 대해 변경하는 경우
2.22  0  1.06  0 
     0  7.38  1.06  0 
     0  1.06  0.07  0 
    5.18  0  3.26  1.6 

m2의 조금 :

,617,451
MCSRmatrix<double> m3 = {{1.01, 0 , 0,0}, {0, 4.07, 0,0},{0,0,6.08,0},{1.06,0,2.2,9.9} }; 

    MCSRmatrix<double> m4 = {{ 1.21, 0 , 1.06,0 }, {0, 3.31, 1.06,0},{0,1.06,-6.01,0},{4.12,0,1.06,-8.3}}; 

    MCSRmatrix<double> m5 = m3+m4 ; 

    cout << m5 ; 

하는 날 올바른 결과를 제공

MCSRmatrix<double> m6 = {{1.01, 0 , 0,0}, {0, 4.07, 0,0},{0,0,6.08,0},{1.06,0,2.2,9.9} }; 
MCSRmatrix<double> m7 = {{ 1.21, 0 , 0,2.15 }, {0, 3.31, 1.03,0},{0,1.06,-6.01,1.01},{4.12,1.1,1.06,-8.3}}; 
cout << m6+m7 ; 

당신이 sum1sum2을 설정하는 추가 루프의 일부에서 terminate called after throwing an instance of 'std::out_of_range'

답변

0

이 프로그램 종료, 당신은 둘 다 낮은 인덱스 값만을 하나를 수행하지 않습니다. 이것은 더 높은 열 번호가 다른 행렬에있을 수 있지만 행에 더있을 수 있기 때문입니다.

또한 sum2에 대한 인덱스 계산은 j2이 아니고 j1이 아니어야합니다.

+0

답장을 보내 주셔서 감사합니다. 더 잘 설명해 주시겠습니까? 더 낮은 인덱스 값을 선택하면 sum1 또는 sum2를 사용해야하는 choese를 어떻게 선택할 수 있습니까? –

+0

sum2 계산에서 인덱스를 변경했습니다 (단, 변경 사항은 없음). 'else if (m1.ja_.at (m1.dim + 1 + j1)! = m2.ja_.at (m2.dim + 1 + j2))'안에'm1 (m1.dim + 1 + j1)! = m2.ja_.at (m2.dim + 1 + j2)''m1.aa_ [index1]'에' m2.aa_ [index2]'이므로,이 2 개의 컴포넌트는 2를 벡터'res.aa_'에 기여합니다 –