Insert element in Modified compressed sparse row matrix

[drudox84]

Hi everybody, i'm developing a Modified Compress Sparse Row Matrix Class (in C++) you can read here : http://www.iue.tuwien.ac.at/phd/wagner/node83.html a short explanation of the method here I wrote the constructor as follow :

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 InvalidSizeException("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;           
          }
      }     
}

and it works well ! I also wrote a findIndex(i,j) method that return the value at i,j of the whole matrix and it works fine, so i'm able to print out the whole matrix using the operator overloading () who return 0 if in this position the element is zero otherwise the right value

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

the operator is defined as :

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

Now I need to write an alghoritm that give 2 index of position is able to insert in the 2 vector the element in the right position could you help me about ?

I wrote this one but doesn't works fine, so that I have wrote a lot of correction for the bug individuated .. ! (the if clause)

template<typename T>
auto constexpr MCSRmatrix<T>::insertAt(itype row ,itype col , T elem) noexcept 
{
      if(elem == 0)
      {
          std::cerr << "zero element to insert ! Exit 1" << std::endl; 
      }
      
      int index = findIndex(row,col);
      if( index != -1  || row == col)
      {     
          aa_.at(index) = elem ;
      }
      else if(index == -1) 
      {

          for(auto i=row; i< dim+1 ; i++) 
          {     
                ja_.at(i) += 1 ;
          }
          
          if(ja_.at(row) > aa_.size() ) //std::cout << aa_.size() << ' ' <<ja_.at(row) << std::endl;
          {   
              
              aa_.push_back(elem);
              ja_.push_back(col);
          }
          else if ( ja_.at(row) <= aa_.size() )
          {    
                if( dim==4 && (ja_.at(2) - ja_.at(1) == 2) && (row==2 ) )
                {
                     //std::cout << "c1" << std::endl;   
                     ja_.insert( ja_.begin() + dim  + row +1 , col );
                     aa_.insert( aa_.begin() + dim  + row +1 , elem);
                }
                else if( dim==4 && (ja_.at(3) - ja_.at(2) == 2) && (row==3 ) )
                {
                     //std::cout <<  << std::endl;   
                     ja_.insert( ja_.begin() + dim  + row +1 , col );
                     aa_.insert( aa_.begin() + dim  + row +1 , elem);
                }
                else  
                {  //std::cout << "i  " << row << "dim+row  " << dim + row << std::endl; 
                  ja_.insert( ja_.begin() + dim  + row , col );
                  aa_.insert( aa_.begin() + dim  + row , elem);
                }

          }
      }
}

for testing the class i post the whole class :

# include <initializer_list>
# include <vector>
# include <iosfwd>
# include <string>
# include <cstdlib>
# include "MatrixException.H"
# include <cassert>
# include <iomanip>


template <typename data_type>
class MCSRmatrix;             // forward declaration

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


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 ;


   public: 
      constexpr MCSRmatrix( std::initializer_list<std::initializer_list<data_type>> rows);
      
      auto constexpr insertAt(itype,itype,data_type) noexcept ;
   
      std::size_t constexpr findIndex(const itype row ,  const itype col) const noexcept ;
      
      const data_type operator()(const itype r , const itype c) const noexcept ;
 
      data_type operator()(const itype r , const itype c) noexcept ;
   
   private:

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

     int dim ; 
};

template <typename T>
std::ostream& operator<<(std::ostream& os ,const MCSRmatrix<T>& m) noexcept 
{
      //m.print();
     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;
}

and the main program for doing testing :

# include "ModCSRmatrix.H"


using namespace std;

int main(){

  MCSRmatrix<int> m1 = {{4,0,1,2,0},{0,1,0,0,0},{0,0,5,7,0},{6,3,0,8,0},{0,0,0,0,9}};    
  cout << "---------------------------------------------------------------------------------------" << endl; 
  cout << m1 ;
  //MCSRmatrix<int> m2 = {{2 ,0 ,1 ,0},{0,7 ,2,0},{0,0,70,0 },{0,0,2,16}};
  //cout << "---------------------------------------------------------------------------------------" << endl;   
  //MCSRmatrix<int> m3 ={{4,0,0,2},{0,1,0,0},{0,0,5,7},{6,3,0,8}};
  cout << "---------------------------------------------------------------------------------------" << endl;   
  m1.insertAt(3,4,16);
  //cout << m1(1,3);
  cout << "---------------------------------------------------------------------------------------" << endl;   
  cout << m1 ;
  cout << "---------------------------------------------------------------------------------------" << endl;   
  MCSRmatrix<int> m2 = {{4,0,1,2},{0,1,0,2},{0,0,5,7},{6,3,0,8}}; 
  cout << "---------------------------------------------------------------------------------------" << endl;   
  //cout << m2 ;
  cout << "---------------------------------------------------------------------------------------" << endl;   
  m2.insertAt(2,3,44);
  cout << m2 ;
  
  return 0;
}