Modified compressed sparse row matrix code

Question

I've implemented a Modified Compressed Sparse Row matrix class! The class works, but I think that some part of the code is quite crap. There are two passages and the whole code became too expansive ... here I report the principal part for which I have some doubt about elegance and is resources consuming:

 template <typename element_type>
class MSRmatrix {

  public:

   using itype = std::size_t ;  

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

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


  public:

   constexpr MSRmatrix( std::initializer_list<std::vector<element_type>> row) noexcept; 

   constexpr MSRmatrix( std::string ) noexcept; 

   constexpr MSRmatrix( std::size_t ) noexcept; 

   auto constexpr printMSR()const noexcept ;

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

   const element_type operator()(const itype, const itype) const noexcept ;

   element_type operator()(const itype, const itype) noexcept ;

   auto operator()(const itype, const itype,const element_type val) noexcept ;

   auto constexpr print() const noexcept ;   

   auto constexpr toFile(std::string) noexcept ;   

   inline auto constexpr nnz() const noexcept { return aa_.size(); }

   inline auto constexpr Dim() const noexcept { return dim ; }
  private:

   itype dim ;  


   std::vector<itype>         ja_ ;   
   std::vector<element_type>  aa_ ;   

   enum value{ stars };

};


template <typename T>
constexpr MSRmatrix<T>::MSRmatrix( std::initializer_list<std::vector<T>> row ) noexcept  : dim{row.size()}
{     
      this->dim = row.size() ;
      auto il = *(row.begin()) ;

      if(this-> dim != il.size())
      {
         std::cerr << "Matrix contained into MSRmatrix format must be square!" << std::endl;    
         exit(-1);   
      }
      std::vector<std::vector<T>> temp(row);

      itype i=0 , j=0 , w=0 , k=0, z=dim;

      this->aa_.resize(dim+1);
      this->ja_.resize(dim+1);

      for(auto& x : row)
      {   
          j=0 ; k=0 ; 
          for(auto& y : x)
          {
             if(i==j) 
             {
                  aa_.at(i) = y ;
             }
             if(i != j && y != 0 ) 
             { 
                aa_.push_back(y);
                z++ ;
                ja_.push_back(j+1);
                if( k==0 )
                { 
                   ja_[w] = z+1; 
                   w++ ;
                   k++ ;
                }
             }
             j++ ;     
          }
          i++ ;  
      }
      while(w <= dim)
      {   
          ja_[w] = z+2 ;  
          w++ ;  
      }

# ifdef __DEBUG__
    printMSR();  
# endif
}

The following method is written to find the index in the original matrix and has the possibility to print out the matrix (reference here)

template <typename T>
std::size_t constexpr MSRmatrix<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;
     }
      itype i = -1;


      if( row == dim )// == ja_.at(row) ) // last row
      {
            i= (ja_.at(row-1)-1) ;


            if(ja_.at(i-1) == col)
              return i-1;
      }

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

and the last thing, which I think is the most important one, is the method written to insert an element into the matrix. I want to improve it.

template <typename T>
auto MSRmatrix<T>::operator()(const itype r, const itype c, const T val) noexcept
{    

   if(val != 0)
   {
     auto i = findIndex(r,c);
     if(i != -1 && i < nnz())
     {     

            aa_.at(i) = val;
     }
     else
     {
        if(r==c)
        {
            aa_.at(r) = val ; // the size doesn't change
        }
        else if( r != c )
        {   
            if( ja_.at(r-1) == ja_.at(r) ) // last line ;
            {     
                  ja_.at(r) += 1 ;

                  i = ja_.at(r-1)-1;
                  ja_.insert(ja_.begin() + i , c  ); 
                  aa_.insert(aa_.begin() + i , val);
            }
            else   // general line
            {
                  for(auto i = ja_.at(r-1)-1 ; i < ja_.at(r)-1 ; i++ )
                  {
                       if( ja_.at(i) == c  )
                       {    
                               aa_.at(i) = val ;
                       }
                       else
                       {     
                            for(auto ii = r ; ii < dim+1 ; ii++ )   
                            {
                              ja_.at(ii) +=1 ;
                            }  
                              ja_.insert(ja_.begin() + i , c  ); 
                              aa_.insert(aa_.begin() + i , val);
                              break;
                       }
                     //std::cout << ja_.at(r) << std::endl ;
                  }
            }
        }
     }
   }
}

in operator() i return T instead of const T& because if the element is not found it return 0 that is a rvalue !

Answer 1

Rethink design decisions

[T]he most important one [..] is the method [..] to insert an element into the matrix.

It isn't. The CSR format is not meant for easy insertion or modification in general. The CSR format is meant for \$\mathcal O(n) \$ matrix vector multiplication. The same holds for the MCSR format.

The sparse matrices are usually a result of a finite element approach: after the initial triangulation only some vertices/edges/faces/nodes will be coupled to others. This leads to a sparse matrix where only some elements will even interact with others. It leads to a single sparse matrix, that won't get changed afterwards.

A modification of a matrix in (M)CSR format will always lead to problems, since the sequential memory has to be moved, leading from an insert that should take \$\mathcal O(1)\$ to an operation that will take \$\mathcal O(n)\$. It's better to use another format for the intermediate manipulation, like DOK* or LIL, especially if you have many insertions. This can change a \$\mathcal O(kn)\$ algorithm into a \$\mathcal O(k*\log(n) + 2n)\$.

Again, the CSR format isn't meant for modification. Don't take my word, though, ask SciPy:

Advantages of the CSR format

• efficient arithmetic operations CSR + CSR, CSR * CSR, etc.
• efficient row slicing
• fast matrix vector products

Disadvantages of the CSR format

• slow column slicing operations (consider CSC)
• changes to the sparsity structure are expensive (consider LIL or DOK)

Therefore, you should change your class design at that point, unless you want to keep your matrix modifiable. It heavily depends on your use case, but if you want to modify your matrix often, then CSC or CSR are the wrong tools.

Use better names

That being said, ja_ and aa_ are possibly two of the worst names you could come up with for members. If you have a look at your code in 3 months you won't remember whether ja_ was a bunch of indices. Give them proper names, like values_ and index_ for example. You can probably come up with better ones.

Consider std::ostream& for output functions.

All your print functions may be eligible to use a std::ostream, but that's not clear from the code you've shown. That way you can re-use functionality, by the way.

Use call by value when you change the value, call by reference otherwise

Prefer const T& instead of T as an argument, unless you change the value in your function. I'm talking about your std::string constructor for example. Since you didn't include the implementation it's not clear whether call by value is needed there.

This only holds for complex values, of course, not for basic types.

---

^{* DOK is plain simple in C++, by the way, you might want to implement it as an exercise.}

Answer 2

Missing headers

#include <cstddef>
#include <iosfwd>
#include <string>
#include <vector>

---

Interface

• I don't think operator<<() can be noexcept, unless you're not using the std::ostream argument!

• Is there a good reason that the initializer_list constructor requires std::vector arguments? The code only requires arguments that are enumerable with range-based for.

• I don't think the findIndex() method ought to be public. It's an implementation detail that client code shouldn't need to be aware of; making it public may be a sign of missing functionality that clients need. I don't know if Dim() and nnz() need to be public, but either way, they could have better names.

• Consider making operator()() return a (possibly const-qualified) reference, to make the interface more like standard collections.

• The return types of public methods should be complete (auto without a trailing return type prevents the methods being called before their definitions are in scope).

• There's no benefit to const on arguments passed by value in the interface. Go ahead and use them in the definitions, of course:

  class MSRmatrix
  {
  public:
      // declaration
      const element_type& operator()(itype, itype) const noexcept;
  };
  
  // definition
  const element_type& MSRmatrix::operator()(const itype, const itype) const noexcept
  {
      // code...
  }
  

• It's probably worth declaring using value_type = element_type like the standard containers - it makes it easier to use in generic code. And consider renaming itype to size_type for the same reason.

• It should be possible to implement print() and toFile() as non-member non-friend functions. If not, consider adding suitable public methods to implement them.

• operator+() wouldn't need to be a friend if you provide operator+= as a member.

• The value enum is never used.

---

Implementation

Constructor

• Don't exit() from within the constructor. Prefer to throw an exception instead - that way the calling code can choose whether to emit output (and if so, by what means) and whether it can recover from the error.
• There's a redundant assignment this->dim = row.size() - the initializer has done its job by the time we get here.
• The temp vector is assigned but never used.
• The index variables can be declared with smaller scope.