c++ - 运算符 + 在稀疏矩阵 : I don't find any bug

Question

大家好，我实现了一个存储在修改后的压缩稀疏行中的稀疏矩阵!构造函数工作正常，我可以验证它，但 operator+ 有一个奇怪的行为:如果我有一个非零值，则总和不会计算出正确的结果。描述了改进的压缩稀疏行方法 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 

但查看 vector aa_(值 vector ，其中第一个 matrix.dim 是对角线元素，后面的元素是对角线元素非零元素看起来矩阵与此结果不同)他们是:

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 

另一个例子，如果我把这一行放到主目录中:

 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 ;

给我正确的结果:

    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 

如果我稍微改变一下 m1 和 m2 :

 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 ;

此程序退出并在抛出“std::out_of_range”实例后调用terminate

Answer 1

在设置 sum1 和 sum2 的添加循环部分，您不想两者都做，只做索引值较低的那个.这是因为较高的列号可能在另一个矩阵中，但在行中更远。

此外，sum2 的索引计算应该使用 j2 而不是 j1。