Matrix Mod. sparse row Operator*(matrix , vector)

Question

I'm implementing a modified compressed sparse row matrix [reference], but I have a problem with Matrix * vector multiplication, I wrote the function but I don't reach to find the bug !

the class used 2 container (std::vector) for store

Diagonal element (aa_[0] to aa_[dim])
the non zero value off-diagonal (aa_[dim+2] to aa_[size_of_non_zero])
pointer of the first element in the row (ja_[0] to ja_[dim] )
in the previous pointer this rules is used : ja_[0]=dim+1 ; ja_[i+1]-ja[i]= number of element in i-th row
column index stored in ja_[ja_[row]] for ja_[row] described above is range is ja[0] to ja[dim+1] ,so the colum index are in ja_[dim+2] to ja_[size_of_non_zero elment]

here the minimal code :

# include <initializer_list>
# include <vector>
# include <iosfwd>
# include <string>
# include <cstdlib>
# include <cassert>
# include <iomanip>
# include <cmath> for(auto i=0; i< A.dim ; i++)
 {
     //for(auto k=A.ja_.at(i) ; k <= A.ja_.at(i+1)-1 ; k++ )
     auto k=A.ja_.at(i)-1; 
     do 
     {    
          b.at(i) += A.aa_.at(k)* x.at(A.ja_.at(k)-1);
          k++ ; for(auto i=0; i< A.dim ; i++)
 {
     //for(auto k=A.ja_.at(i) ; k <= A.ja_.at(i+1)-1 ; k++ )
     auto k=A.ja_.at(i)-1; 
     do 
     {    
          b.at(i) += A.aa_.at(k)* x.at(A.ja_.at(k)-1);
          k++ ;
     }while (k < A.ja_.at(i+1)-1 ); // ;
 }
 return b;

     }while (k < A.ja_.at(i+1)-1 ); // ;
 }
 return b;

# include <set>
# include <fstream>



  template <typename data_type>
    class MCSRmatrix {
       public:
             using itype = std::size_t ;

    template <typename T>
          friend std::vector<T> operator*(const MCSRmatrix<T>& A, const std::vector<T>& x ) noexcept ;

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


    private:

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

         int dim ; 
    };

    //constructor 
    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()) for(auto i=0; i< A.dim ; i++)
 {
     //for(auto k=A.ja_.at(i) ; k <= A.ja_.at(i+1)-1 ; k++ )
     auto k=A.ja_.at(i)-1; 
     do 
     {    
          b.at(i) += A.aa_.at(k)* x.at(A.ja_.at(k)-1);
          k++ ;
     }while (k < A.ja_.at(i+1)-1 ); // ;
 }
 return b;

          {
              throw std::runtime_error("error matrix must be square");
          }

          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;           
              }
          }     
      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>
    std::vector<T> operator*(const MCSRmatrix<T>& A, const std::vector<T>& x ) noexcept 
    {     

         std::vector<T> b(A.dim); 
         for(auto i=0; i < A.dim ; i++ )
             b.at(i) = A.aa_.at(i)* x.at(i) ;   


         for(auto i=0; i< A.dim ; i++)
         {
             for(auto k=A.ja_.at(i) ; k < A.ja_.at(i+1)-1 ; k++ )
             {    
                  b.at(i) += A.aa_.at(k)* x.at(A.ja_.at(k));
             }   
         }
         return b;
    }

and finally the main

# include "ModCSRmatrix.H"


using namespace std;

int main(){
   std::vector<double> v1={0,1.3,4.2,0.8};
   MCSRmatrix<double> m1  = {{1.01, 0 , 2.34,0}, {0, 4.07, 0,0},{3.12,0,6.08,0},{1.06,0,2.2,9.9} }; 
    std::vector<double> v2 = m1*v1 ;

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

but the result is different from the result obtain in octave !

I've correct the code and now compile ! it give me the result :

0 5.291 25.536 9.68

but the correct result obtained using octave is :

9.8280 5.2910 25.5360 17.1600

the strange thing is that the same code written in Fortran works!

MODULE MSR
 IMPLICIT NONE

CONTAINS
     subroutine amuxms (n, x, y, a,ja)
      real*8  x(*), y(*), a(*)
      integer n, ja(*)
      integer i, k
      do 10 i=1, n
        y(i) = a(i)*x(i)
 10     continue
      do 100 i = 1,n

         do 99 k=ja(i), ja(i+1)-1
            y(i) = y(i) + a(k) *x(ja(k))
 99      continue
 100  continue

      return

      end

END MODULE

PROGRAM MSRtest
USE MSR
IMPLICIT NONE
INTEGER :: i
REAL(KIND(0.D0)), DIMENSION(4) :: y, x= (/0.,1.3,4.2,0.8/)

REAL(KIND(0.D0)), DIMENSION(9) :: AA = (/ 1.01, 4.07, 6.08, 9.9, 0., 2.34, 3.12, 1.06, 2.2/) 
INTEGER , DIMENSION(9)         :: JA = (/6, 7, 7, 8, 10, 3, 1, 1, 3/)



WRITE(6,FMT='(4F8.3)') (x(I), I=1,4)    

CALL amuxms(4,x,y,aa,ja)

WRITE(6,FMT='(4F8.3)') (y(I), I=1,4)    

END PROGRAM

in the above code the value of aa and ja is given by the c++ constructor putting this member

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

and call it at the end of constructor! now I have added the lines of the member at the end of constructor so if you try the constructor you get exactly the same vector written in the fortran code

thanks I followed your advice @Paul H. and rewrite the operator + as follow: (I didn't change the ja_ indexing because in my class I have a lot of already more or less un-bugged method )

template <typename T>
std::vector<T> operator*(const MCSRmatrix<T>& A, const std::vector<T>& x ) noexcept 
{     

     std::vector<T> b(A.dim); 
     for(auto i=0; i < A.dim ; i++ )
         b.at(i) = A.aa_.at(i)* x.at(i) ;   


     for(auto i=0; i< A.dim ; i++)
     {
         //for(auto k=A.ja_.at(i) ; k <= A.ja_.at(i+1)-1 ; k++ )
         auto k=A.ja_.at(i)-1; 
         do 
         {    
              b.at(i) += A.aa_.at(k)* x.at(A.ja_.at(k)-1);
              k++ ;
         }while (k < A.ja_.at(i+1)-1 ); // ;
     }
     return b;
}

as You can see I have subtracts 1 from all ja_ using as indices :

x.at(A.ja_.at(k)-1) instead of x.at(A.ja_.at(k))
different start of index K k=A.ja_.at(i)-1
and different end of cicle (I've used a do while instead of for)

Paul H. Paul H. · Accepted Answer · 2017-11-18T23:53:05

The debugger is your friend! For future reference, here is a link to a very good blog post on debugging small programs: How to debug small programs. There are a couple of off by one mistakes in your code. If you create the 4 x 4 matrix used as an example in the reference you linked to, you will see that the ja_ values you calculate are all off by one. The reason your Fortran version works is because arrays in Fortran are by default indexed starting from 1, not 0. So in class MCSRmatrix change

ja_.at(w) = dim+2;

to

ja_.at(w) = dim+1;

and

ja_.push_back(j);

to

ja_.push_back(j-1);

Then in your operator* method change

for(auto k=A.ja_.at(i) ; k < A.ja_.at(i+1)-1 ; k++ )

to

for(auto k = A.ja_.at(i); k < A.ja_.at(i+1); k++)

Matrix Mod. sparse row Operator*(matrix , vector)

1 Answers