# include <initializer_list>
# include <iomanip>
# include <iostream>
# include <cassert>
# include <cmath>
# include <vector>
# include <string>
# include <set>
# include <algorithm>
# include <map>
# include <exception>

# define __TESTING__

template<typename T>
class DIAmatrix;


template <typename U>
std::ostream& operator<<(std::ostream& os, DIAmatrix<U>& m );

//- SpMV
//
template<typename U>
std::vector<U> operator*(const DIAmatrix<U>& , const std::vector<U>& ) ;


template<typename Type >
class DIAmatrix
{
    template <typename U>
    friend std::ostream& operator<<(std::ostream& os, DIAmatrix<U>& m );

     //- SpMV
     //
     template<typename U>
     friend std::vector<U> operator*(const DIAmatrix<U>& , const std::vector<U>& ) ;


  public:

     constexpr DIAmatrix(std::initializer_list<std::vector<Type>> ) ;

     Type& operator()(const std::size_t , const std::size_t ) noexcept ;

      auto constexpr printDIA() const noexcept ;


    private:

      std::size_t denseRows ;
      std::size_t denseCols ;

      mutable Type dummy    ;

      std::size_t nnz ;

      std::size_t dim ;

      std::map<int ,  std::vector<Type>> value ;   // value
      std::set<int>                        dig ;   // off-diagonal distance

      Type constexpr findValue(const std::size_t , const std::size_t ) const noexcept ;
};

//-------       implemantation      ------- //

template<typename T>
constexpr DIAmatrix<T>::DIAmatrix(std::initializer_list<std::vector<T>> rows )
{
    denseRows = rows.size();
    denseCols = (*rows.begin()).size();

    if(denseRows != denseCols )
    {
       throw std::runtime_error("DIA-Matrix , SIZE EXCEPTION THROWN :n>>> Matrix Must be square! <<<");
    }

    dim = denseRows ;

    value[0].resize(dim);
    value[1].resize(dim-1);       // diagonal
    value[-1].resize(dim-1);      // tri-bands as default

    dig.insert(-1);
    dig.insert(0);                // tribands as default
    dig.insert(1);


    auto i =0 , j=0 ;
    for(auto row : rows )
    {
       j=0;
       for(auto val : row )
       {
         if(val != 0)
         {
                if( i==j )   value[0].at(i) = val  ;
           else if(j==i-1)   value[-1].at(j) = val  ;
           else if(j==i+1)   value[1].at(i) = val ;
           else if(j < i )
           {
              auto offset =  j-i ;
              if( !(dig.find(offset) != dig.end()) )
              {
                  dig.insert(offset);
                  value[offset].resize(dim-(fabs(j-i)));
                  value[offset].at(j) = val ;
              }
              else
                  value[offset].at(j) = val ;
           }
           else if(j > i)
           {
             auto offset = j-i ;
             if( !(dig.find(offset) != dig.end()))
             {
                 dig.insert(offset);
                 value[offset].resize(dim-(fabs(j-i)));
                 value[offset].at(i) = val ;
             }
             else
               value[offset].at(i) = val;
           }
         }
         j++;
       }
       i++;
    }
#  ifdef __TESTING__
      printDIA();
#  endif
}


template<typename T>
T constexpr DIAmatrix<T>::findValue(const std::size_t r, const std::size_t c) const noexcept
{
    T val = 0.0;
    if( value.find(c-r) != value.end() )
       val  = r < c ? value.at(c-r).at(r) : value.at(c-r).at(c) ;
    return val ;
}


template<typename T>
auto constexpr DIAmatrix<T>::printDIA() const noexcept
{

   auto it = dig.begin() ;

   for( auto& n : dig)
    {
      std::cout << std::setw(4) <<*it << " | :  " ;
      for(auto& x : value.at(n) )
      {
         std::cout << x << " " ;
      }
      ++it ;
      std::cout << std::endl;
    }
}

template <typename T>
T& DIAmatrix<T>::operator()(const std::size_t r, const std::size_t c) noexcept
{
      assert( r >0 && r <= dim && c > 0 && c <= dim ) ;
      dummy = findValue(r-1,c-1);
      return dummy ;
}


template <typename U>
std::ostream& operator<<(std::ostream& os, DIAmatrix<U>& m )
{
      for(auto i=1; i <= m.denseRows ; i++ )
      {
          for(auto j=1 ; j <= m.denseCols ; j++)
          {
               os << std::setw(6) << m(i,j) << " ";
          }
          os << std::endl;
      }
}



template<typename T>
std::vector<T> operator*(const DIAmatrix<T>& m, const std::vector<T>& x )
{
    std::vector<T> y(x.size());

    if(m.denseCols != x.size())
    {
       throw std::runtime_error("Erorr size of matrix and vector doesn't match!");
    }
    else
    {
      // for(auto i=0; i< x.size() ; i++ )
          for(auto offset : m.dig )
          {
            for(auto j = 0 ; j < m.value.at(offset).size() ; j++ )
            {
                 y.at(j) += m.value.at(offset).at(j) * x.at(j+fabs(offset))  ;
            }
         }
      // }
       //}
    }
    return y;
}

# include "DiaMatrix.H"

using namespace std;

int main() {

   DIAmatrix<int> dm = {{11,12,0,0,0,0,0,0} ,{0,22,0,0,0,0,0,0} ,{31,32,33,0,0,0,0,0},
                              {41,42,43,44,0,0,0,0}, {0,0,0,0,55,56,0,0},{0,0,0,0,0,66,67,0},{0,0,0,0,0,0,77,78},{0,0,0,0,0,0,87,88}};

     cout << dm ;

    std::vector<int> v1 = {3,4,0,1,6,8,1,19};


    vector<int> v2 = dm * v1 ;

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


  return 0;
}

122 130 43 44 778 595 3212 1672

81 88 221 335 778 595 1559 1759