PDS_TriangleCount/inc/impl.hpp

/**
 * \file    impl.hpp
 * \brief   Implementation common header file
 *
 * \author
 *    Christos Choutouridis AEM:8997
 *    <cchoutou@ece.auth.gr>
 */
#ifndef IMPL_HPP_
#define IMPL_HPP_

#include <type_traits>
#include <utility>
#include <algorithm>
#include <vector>
#include <tuple>
#include <cstddef>
#include <iostream>
#include <fstream>

using std::size_t;

/*
 * Small helper to strip types
 */
template<typename T>
struct remove_cvref {
    typedef std::remove_cv_t<std::remove_reference_t<T>> type;
};
template<typename T>
using remove_cvref_t = typename remove_cvref<T>::type;

/*!
 * Enumerator to denote the storage type of the array to use.
 */
enum class MatrixType {
   FULL,          /*!< Matrix is asymmetric */
   SYMMETRIC,     /*!< Matrix is symmetric */
};

/*
 * Forward type declarations
 */
template<typename DataType, typename IndexType, MatrixType Type = MatrixType::SYMMETRIC> struct Matrix;
template<typename DataType, typename IndexType, MatrixType Type = MatrixType::SYMMETRIC> struct SpMat;
template<typename DataType, typename IndexType> struct SpMatCol;
template<typename DataType, typename IndexType> struct SpMatRow;
template<typename DataType, typename IndexType> struct SpMatVal;

 /*!
 * 2D-array wrapper for v1 and v2 part of the exercise used as RAII
 * and copy-prevention.
 *
 * This is a very thin abstraction layer over a native array.
 * This is tested using compiler explorer and our template produce
 * almost identical assembly.
 *
 * The penalty hit we have is due to the fact that we use a one dimension array
 * and we have to calculate the actual position from an (i,j) pair.
 * The use of 1D array was our intention from the beginning, so the penalty
 * was pretty much unavoidable.
 *
 * \tparam DataType  The underling data type of the array
 * \tparam Type      The storage type of the array
 *    \arg  FULL              For full matrix
 *    \arg  SYMMETRIC         For symmetric matrix (we use only the lower part)
 */
template<typename DataType, typename IndexType, MatrixType Type>
struct Matrix {

   using dataType    = DataType;                   //!< meta:export of underling data type
   using indexType   = IndexType;                  //!< meta:export of underling index type
   static constexpr MatrixType matrixType = Type;  //!< export of array type

   /*!
    * \name Obj lifetime
    */
   //! @{

   //! Constructor using data type and size
   Matrix(DataType* t, IndexType s) noexcept : m_(t), size_(s) { }

   //! RAII deleter
   ~Matrix() { delete m_; }
   //! move ctor
   Matrix(Matrix&& a) noexcept { a.swap(*this); }
   //! move
   Matrix& operator=(Matrix&& a) noexcept { a.swap(*this); return *this; }
   Matrix(const Matrix& a)             = delete;  //!< No copy ctor
   Matrix& operator=(const Matrix& a)  = delete;  //!< No copy
   //! @}

   //! \name Data exposure
   //! @{

   //! memory capacity of the matrix with diagonal data
   template<MatrixType AT=Type> std::enable_if_t<AT==MatrixType::SYMMETRIC, IndexType>
      static constexpr capacity(IndexType N) noexcept { return (N+1)*N/2; }
   //! memory capacity for full matrix
   template<MatrixType AT=Type> std::enable_if_t<AT==MatrixType::FULL, IndexType>
      static constexpr capacity(IndexType N) noexcept { return N*N; }

   //! Get the size of each dimension
   IndexType size() noexcept { return size_; }

   /*
    * virtual 2D accessors
    */
   template<MatrixType AT=Type>
   std::enable_if_t <AT==MatrixType::SYMMETRIC, DataType>
   get (IndexType i, IndexType j) {
      if (i<j)    return m_[j*(j+1)/2 + i];  // Upper, use opposite index
      else        return m_[i*(i+1)/2 + j];  // Lower, use our notation
   }
   template<MatrixType AT=Type>
   std::enable_if_t <AT==MatrixType::FULL, DataType>
   get(IndexType i, IndexType j) {
      return m_[i*size_ + j];
   }
   template<MatrixType AT=Type>
   std::enable_if_t <AT==MatrixType::SYMMETRIC, DataType>
   set(DataType v, IndexType i, IndexType j) {
      if (i<j)    return m_[j*(j+1)/2 + i] =v;  // Upper, use opposite index
      else        return m_[i*(i+1)/2 + j] =v;  // Lower, use our notation
   }
   template<MatrixType AT=Type>
   std::enable_if_t <AT==MatrixType::FULL, DataType>
   set(DataType v, IndexType i, IndexType j) {
      return m_[i*size_ + j] =v;
   }
   DataType operator()(IndexType i, IndexType j) { return get(i, j); }

   // a basic serial iterator support
   DataType* begin() noexcept { return m_; }
   DataType* end()   noexcept { return m_ + Matrix::capacity(size_); }
   //! @}

   /*!
    * \name Safe iteration API
    *
    * This api automates the iteration over the array based on
    * MatrixType
    */
   //! @{
   template<typename F, typename... Args>
   void for_each_in (IndexType begin, IndexType end, F&& lambda, Args&&... args) {
      for (IndexType it=begin ; it<end ; ++it) {
         std::forward<F>(lambda)(std::forward<Args>(args)..., it);
      }
   }
   //! @}

   // move helper
   void swap(Matrix& src) noexcept {
      std::swap(m_, src.m_);
      std::swap(size_, src.size_);
   }
private:
   DataType*   m_;      //!< Pointer to actual data.
   IndexType   size_;   //!< the virtual size of each dimension.
};

 /*!
 * RAII allocation helper, the smart_ptr way.
 */
template<typename DataType, typename IndexType, MatrixType Type = MatrixType::SYMMETRIC>
Matrix<DataType, IndexType, Type> make_Matrix(IndexType s) {
   return Matrix<DataType, IndexType, Type>(new DataType[Matrix<DataType, IndexType, Type>::capacity(s)], s);
}


/**
 * A simple sparse matrix implementation.
 *
 * We use CSC format and provide get/set functionalities for each (i,j) item
 * on the matrix. We also provide a () overload using a proxy SpMatVal object.
 * This way the user can:
 * \code
 *    auto v = A(3,4);
 *    A(3, 4) = 7;
 * \endcode
 *
 * We also provide getCol() and getRow() functions witch return a viewer/iterator to rows and
 * columns of the matrix. In the case of a symmetric matrix instead of a row we return the
 * equivalent column. This way we gain speed due to CSC format nature.
 *
 * @tparam DataType  The type for values
 * @tparam IndexType The type for indexes
 * @tparam Type      The Matrix type (FULL or SYMMETRIC)
 */
template<typename DataType, typename IndexType, MatrixType Type>
struct SpMat {
   using dataType    = DataType;                   //!< meta:export of underling data type
   using indexType   = IndexType;                  //!< meta:export of underling index type
   static constexpr MatrixType matrixType = Type;  //!< export of array type

   friend struct SpMatCol<DataType, IndexType>;
   friend struct SpMatRow<DataType, IndexType>;
   friend struct SpMatVal<DataType, IndexType>;

   /*!
    * \name Obj lifetime
    */
   //! @{

   //! Default ctor with optional memory allocations
   SpMat(IndexType n=IndexType{}, IndexType nnz=IndexType{}) :
      values(nnz, DataType{}),
      rows(nnz, IndexType{}),
      col_ptr((n)? n+1:2, IndexType{}),
      N(n),
      NNZ(nnz) { }

   //! A ctor using csc array data
   SpMat(IndexType n, IndexType nnz, const IndexType* row, const IndexType* col) :
      values(nnz, 1),
      rows(row, row+nnz),
      col_ptr(col, col+n+1),
      N(n),
      NNZ(nnz) { }

   //! ctor using csc array data with value array
   SpMat(IndexType n, IndexType nnz, const DataType* v, const IndexType* row, const IndexType* col) :
      values(v, v+nnz),
      rows(row, row+nnz),
      col_ptr(col, col+n+1),
      N(n),
      NNZ(nnz) { }

   //! ctor vectors of row/col and default value for values array
   SpMat(IndexType n, IndexType nnz, const DataType v, const std::vector<IndexType>& row, const std::vector<IndexType>& col) :
      values(nnz, v),
      rows (row),
      col_ptr(col),
      N(n),
      NNZ(nnz) { }

   //! move ctor
   SpMat(SpMat&& a) noexcept { moves(std::move(a)); }
   //! move
   SpMat& operator=(SpMat&& a) noexcept { moves(std::move(a)); return *this; }
   SpMat(const SpMat& a)             = delete;  //!< make sure there are no copies
   SpMat& operator=(const SpMat& a)  = delete;  //!< make sure there are no copies
   //! @}

   //! \name Data exposure
   //! @{

   //! \return the dimension of the matrix
   IndexType size()     noexcept { return N; }
   //! After construction size configuration tool
   IndexType size(IndexType n) {
      col_ptr.resize(n+1);
      return N = n;
   }
   //! \return the NNZ of the matrix
   IndexType capacity() noexcept { return NNZ; }
   //! After construction NNZ size configuration tool
   IndexType capacity(IndexType nnz) {
      values.reserve(nnz);
      rows.reserve(nnz);
      return NNZ;
   }
   // getters for row arrays of the struct (unused)
   std::vector<DataType>&  getValues() noexcept { return values; }
   std::vector<IndexType>& getRows() noexcept { return rows; }
   std::vector<IndexType>& getCols() noexcept { return col_ptr; }

   /*!
    * Return a proxy SpMatVal object with read and write capabilities.
    * @param i    The row number
    * @param j    The column number
    * @return     tHE SpMatVal object
    */
   SpMatVal<DataType, IndexType> operator()(IndexType i, IndexType j) {
      return SpMatVal<DataType, IndexType>(this, get(i, j), i, j);
   }

   /*!
    * A read item functionality using binary search to find the correct row
    *
    * @param i    The row number
    * @param j    The column number
    * @return     The value of the item or DataType{} if is not present.
    */
   DataType get(IndexType i, IndexType j) {
      IndexType idx; bool found;
      std::tie(idx, found) =find_idx(rows, col_ptr[j], col_ptr[j+1], i);
      return (found) ? values[idx] : 0;
   }

   /*!
    * A read item functionality using linear search to find the correct row
    *
    * @param i    The row number
    * @param j    The column number
    * @return     The value of the item or DataType{} if is not present.
    */
   DataType get_lin(IndexType i, IndexType j) {
      IndexType idx; bool found;
      std::tie(idx, found) =find_place_idx(rows, col_ptr[j], col_ptr[j+1], i);
      return (found) ? values[idx] : 0;
   }

   /*!
    * A write item functionality.
    *
    * First we search if the matrix has already a value in (i, j) position.
    * If so we just change it to a new value. If not we add the item on the matrix.
    *
    * @note
    *    When change a value, we don't increase the NNZ value of the struct. We expect the user has already
    *    change the NNZ value to the right one using @see capacity() function. When adding a value we
    *    increase the NNZ.
    *
    * @param i    The row number
    * @param j    The column number
    * @return     The new value of the item .
    */
   DataType set(DataType v, IndexType i, IndexType j) {
      IndexType idx; bool found;
      std::tie(idx, found) = find_place_idx(rows, col_ptr[j], col_ptr[j+1], i);
      if (found)
         return values[idx] = v;    // we don't change NNZ even if we write "0"
      else {
         values.insert(values.begin()+idx, v);
         rows.insert(rows.begin()+idx, i);
         std::transform(col_ptr.begin()+j+1, col_ptr.end(), col_ptr.begin()+j+1, [](IndexType it) {
            return ++it;
         });
         ++NNZ;                     // we increase the NNZ even if we write "0"
         return v;
      }
   }

   /*!
    * Get a view of a CSC column
    * @param j    The column to get
    * @return     The SpMatCol object @see SpMatCol
    */
   SpMatCol<DataType, IndexType> getCol(IndexType j) {
      return SpMatCol<DataType, IndexType>(this, col_ptr[j], col_ptr[j+1]);
   }

   /*!
    * Get a view of a CSC row
    *
    * In case of a SYMMETRIC matrix we can return a column instead.
    *
    * @param j    The row to get
    * @return     The SpMatCol object @see SpMatCol
    */
   template<MatrixType AT= Type>
   std::enable_if_t<AT==MatrixType::SYMMETRIC, SpMatCol<DataType, IndexType>>
   getRow(IndexType i) {
      return getCol(i);
   }

   /*!
    * Get a view of a CSC row
    *
    * @param j    The row to get
    * @return     The SpMatRow object @see SpMatRow
    */
   template<MatrixType AT= Type>
   std::enable_if_t<AT==MatrixType::FULL, SpMatCol<DataType, IndexType>>
   getRow(IndexType i) {
      return SpMatRow<DataType, IndexType>(this, i);
   }

   // values only iterator support
   DataType* begin() noexcept { return values.begin(); }
   DataType* end()   noexcept { return values.end(); }
   //! @}

   //! A small iteration helper
   template<typename F, typename... Args>
   void for_each_in (IndexType begin, IndexType end, F&& lambda, Args&&... args) {
      for (IndexType it=begin ; it<end ; ++it) {
         std::forward<F>(lambda)(std::forward<Args>(args)..., it);
      }
   }

   // friend operations for printing
   template<typename D, typename I, MatrixType T> friend void print(SpMat<D, I, T>& mat);
   template<typename D, typename I, MatrixType T> friend void print_dense(SpMat<D, I, T>& mat);

private:
   /*!
    * A small binary search implementation using index for begin-end instead of iterators.
    *
    * \param   v     Reference to vector to search
    * \param   begin The vector's index to begin
    * \param   end   The vector's index to end
    * \param   match What to search
    * \return  An <index, status> pair.
    *                index    is the index of the item or end if not found
    *                status   is true if found, false otherwise
    */
   std::pair<IndexType, bool> find_idx(const std::vector<IndexType>& v, IndexType begin, IndexType end, IndexType match) {
      IndexType b = begin, e = end-1;
      while (b <= e) {
         IndexType m = (b+e)/2;
         if       (v[m] == match)   return  std::make_pair(m, true);
         else if  (b >= e)          return  std::make_pair(end, false);
         else {
            if    (v[m] <  match)   b = m +1;
            else                    e = m -1;
         }
      }
      return std::make_pair(end, false);
   }
   /*!
    * find helper for set using index for begin-end instead of iterators.
    *
    * We search for the item or a place to add the item.
    * So we return the index if we find it. Otherwise we return the place to add it.
    *
    * \param   v     Reference to vector to search
    * \param   begin The vector's index to begin
    * \param   end   The vector's index to end
    * \param   match What to search
    * \return  An <index, status> pair.
    *                index    is the index of the item or end if not found
    *                status   is true if found, false otherwise
    */
   std::pair<IndexType, bool> find_place_idx(const std::vector<IndexType>& v, IndexType begin, IndexType end, IndexType match) {
      for ( ; begin < end ; ++begin) {
         if       (match == v[begin])   return std::make_pair(begin, true);
         else if  (match <  v[begin])   return std::make_pair(begin, false);
      }
      return std::make_pair(end, false);
   }
   // move helper
   void moves(SpMat&& src) noexcept {
      values  = std::move(src.values);
      rows    = std::move(src.rows);
      col_ptr = std::move(src.col_ptr);
      N       = std::move(src.N);   // redundant for primitives
      NNZ     = std::move(src.NNZ); //
   }
   //! \name Data
   //! @{
   std::vector<DataType>   values {};     //!< vector to store the values of the matrix
   std::vector<IndexType>  rows{};        //!< vector to store the row information
   std::vector<IndexType>  col_ptr{1,0};  //!< vector to store the column pointers
   IndexType   N{0};                      //!< The dimension of the matrix (square)
   IndexType   NNZ{0};                    //!< The NNZ (capacity of the matrix)
   //! @}
};

/*!
 * A view/iterator hybrid object for SpMat columns.
 *
 * This object provides access to a column of a SpMat. The public functionalities
 * allow data access using indexes instead of iterators. We prefer indexes over iterators
 * because we can apply the same index to different inner vector of SpMat without conversion.
 *
 * @tparam DataType
 * @tparam IndexType
 */
template<typename DataType, typename IndexType>
struct SpMatCol {
   using owner_t = SpMat<DataType, IndexType>;

   /*!
    * ctor using column pointers for begin-end. own is pointer to SpMat.
    */
   SpMatCol(owner_t* own, const IndexType begin, const IndexType end) noexcept :
      owner_(own), index_(begin), begin_(begin), end_(end) {
      vindex_ = vIndexCalc(index_);
   }
   SpMatCol()                          = default;
   SpMatCol(const SpMatCol&)           = delete;   //!< make sure there are no copies
   SpMatCol& operator=(const SpMatCol&)= delete;   //!< make sure there are no copies
   SpMatCol(SpMatCol&&)                = default;
   SpMatCol& operator=(SpMatCol&&)     = default;

   //! a simple dereference operator, like an iterator
   DataType operator* () {
      return get();
   }
   //! Increment operator acts on index(), like an iterator
   SpMatCol& operator++ ()    { advance(); return *this; }
   SpMatCol& operator++ (int) { SpMatCol& p = *this; advance(); return p; }

   //! () operator acts as member access (like a view)
   DataType operator()(IndexType x) {
      return (x == index())? get() : DataType{};
   }
   //! = operator acts as member assignment (like a view)
   DataType operator= (DataType v) { return owner_->values[index_] = v; }
   // iterator like handlers
   // these return a virtual index value based on the items position on the full matrix
   // but the move of the index is just a ++ away.
   IndexType       index()       noexcept { return vindex_; }
   const IndexType index() const noexcept { return vindex_; }
   IndexType       begin()       noexcept { return vIndexCalc(begin_); }
   const IndexType begin() const noexcept { return vIndexCalc(begin_); }
   IndexType       end()         noexcept { return owner_->N; }
   const IndexType end()   const noexcept { return owner_->N; }

   /*!
    * Multiplication operator
    *
    * We follow only the non-zero values and multiply only the common indexes.
    *
    * @tparam C   Universal reference for the type right half site column
    *
    * @param c    The right hand site matrix
    * @return     The value of the inner product of two vectors
    * @note       The time complexity is \$ O(nnz1+nnz2) \$.
    *             Where the nnz is the max NNZ elements of the column of the matrix
    */
   template <typename C>
   DataType operator* (C&& c) {
      static_assert(std::is_same<remove_cvref_t<C>, SpMatCol<DataType, IndexType>>(), "");
      DataType v{};
      while (index() != end() && c.index() != c.end()) {
         if      (index() < c.index())  advance(); // advance me
         else if (index() > c.index())  ++c;       // advance other
         else { //index() == c.index()
            v += get() * *c;                       // multiply and advance both
            ++c;
            advance();
         }
      }
      return v;
   }

private:
   //! small tool to increase the index pointers to SpMat matrix
   void advance() noexcept {
      ++index_;
      vindex_ = vIndexCalc(index_);
   }
   //! tool to translate between col_ptr indexes and SpMat "virtual" full matrix indexes
   IndexType vIndexCalc(IndexType idx) {
      return (idx < end_) ? owner_->rows[idx] : end();
   }
   //! small get tool
   DataType get() { return owner_->values[index_]; }

   owner_t*    owner_   {nullptr};     //!< Pointer to owner SpMat. SpMatCol is just a view
   IndexType   vindex_  {IndexType{}}; //!< Virtual index of full matrix
   IndexType   index_   {IndexType{}}; //!< index to SpMat::rows
   IndexType   begin_   {IndexType{}}; //!< beginning index of the column in SpMat::rows
   IndexType   end_     {IndexType{}}; //!< ending index of the column in SpMat::rows
};

/*!
 * A view/iterator hybrid object for SpMat rows.
 *
 * This object provides access to a column of a SpMat. The public functionalities
 * allow data access using indexes instead of iterators. We prefer indexes over iterators
 * because we can apply the same index to different inner vector of SpMat without conversion.
 *
 * @tparam DataType
 * @tparam IndexType
 */
template<typename DataType, typename IndexType>
struct SpMatRow {
   using owner_t = SpMat<DataType, IndexType>;

   /*!
    * ctor using virtual full matrix row index. own is pointer to SpMat.
    */
   SpMatRow(owner_t* own, const IndexType row) noexcept :
      owner_(own), vindex_(IndexType{}), row_(row), index_(IndexType{}),
      begin_(IndexType{}), end_(owner_->NNZ) {
      // place begin
      while(begin_ != end_ && owner_->rows[begin_] != row_)
         ++begin_;
      // place index_ and vindex_
      if (owner_->rows[index_] != row_)
         advance();
   }
   SpMatRow()                          = default;
   SpMatRow(const SpMatRow&)           = delete;   //!< make sure there are no copies
   SpMatRow& operator=(const SpMatRow&)= delete;   //!< make sure there are no copies
   SpMatRow(SpMatRow&&)                = default;
   SpMatRow& operator=(SpMatRow&&)     = default;

   //! a simple dereference operator, like an iterator
   DataType operator* () {
      return get();
   }
   //! Increment operator acts on index(), like an iterator
   //! here the increment is a O(N) process.
   SpMatRow& operator++ ()    { advance(); return *this; }
   SpMatRow& operator++ (int) { SpMatRow& p = *this; advance(); return p; }

   //! () operator acts as member access (like a view)
   DataType operator()(IndexType x) {
      return (x == index())? get() : DataType{};
   }
   //! = operator acts as member assignment (like a view)
   DataType operator= (DataType v) { return owner_->values[index_] = v; }
   // iterator like handlers
   // these return a virtual index value based on the items position on the full matrix
   // but the move of the index is just a ++ away.
   IndexType       index()       noexcept { return vindex_; }
   const IndexType index() const noexcept { return vindex_; }
   IndexType       begin()       noexcept { return vIndexCalc(begin_); }
   const IndexType begin() const noexcept { return vIndexCalc(begin_); }
   IndexType       end()         noexcept { return owner_->N; }
   const IndexType end()   const noexcept { return owner_->N; }

   /*!
    * Multiplication operator
    *
    * We follow only the non-zero values and multiply only the common indexes.
    *
    * @tparam C   Universal reference for the type right half site column
    *
    * @param c    The right hand site matrix
    * @return     The value of the inner product of two vectors
    * @note       The time complexity is \$ O(N+nnz2) \$ and way heavier the ColxCol multiplication.
    *             Where the nnz is the max NNZ elements of the column of the matrix
    */
   template <typename C>
   DataType operator* (C&& c) {
      static_assert(std::is_same<remove_cvref_t<C>, SpMatCol<DataType, IndexType>>(), "");
      DataType v{};
      while (index() != end() && c.index() != c.end()) {
         if      (index() < c.index())  advance(); // advance me
         else if (index() > c.index())  ++c;       // advance other
         else { //index() == c.index()
            v += get() * *c;                       // multiply and advance both
            ++c;
            advance();
         }
      }
      return v;
   }
private:
   //! small tool to increase the index pointers to SpMat matrix
   //! We have to search the entire rows vector in SpMat to find the next
   //! virtual row position.
   //! time complexity O(N)
   void advance() noexcept {
      do
         ++index_;
      while(index_ != end_ && owner_->rows[index_] != row_);
      vindex_ = vIndexCalc(index_);
   }
   //! tool to translate between col_ptr indexes and SpMat "virtual" full matrix indexes
   IndexType vIndexCalc(IndexType idx) {
      for(IndexType i =0 ; i<(owner_->N+1) ; ++i)
         if (idx < owner_->col_ptr[i])
            return i-1;
      return end();
   }
   //! small get tool
   DataType get() { return owner_->values[index_]; }

   owner_t*    owner_   {nullptr};     //!< Pointer to owner SpMat. SpMatCol is just a view
   IndexType   vindex_  {IndexType{}}; //!< Virtual index of full matrix
   IndexType   row_     {IndexType{}}; //!< The virtual full matrix row of the object
   IndexType   index_   {IndexType{}}; //!< index to SpMat::rows
   IndexType   begin_   {IndexType{}}; //!< beginning index of the column in SpMat::rows
   IndexType   end_     {IndexType{}}; //!< ending index of the column in SpMat::rows
};

/*!
 * A proxy SpMat value object/view.
 *
 * This object acts as proxy to provide read/write access to an SpMat item.
 *
 * @tparam DataType     The type of the values of the SpMat matrix
 * @tparam IndexType    The type of the indexes of the SpMat matrix
 */
template<typename DataType, typename IndexType>
struct SpMatVal {
   using owner_t = SpMat<DataType, IndexType>;

   //!< ctor using all value-row-column data, plus a pointer to owner SpMat object
   SpMatVal(owner_t* own, DataType v, IndexType i, IndexType j) :
      owner_(own), v_(v), i_(i), j_(j) { }
   SpMatVal()                           = default;
   SpMatVal(const SpMatVal&)            = delete;  //!< make sure there are no copies
   SpMatVal& operator=(const SpMatVal&) = delete;  //!< make sure there are no copies
   SpMatVal(SpMatVal&&)                 = default;
   SpMatVal& operator=(SpMatVal&&)      = default;

   //! Operator to return the DataType value implicitly
   operator DataType() { return v_; }
   //! Operator to write back to owner the assigned value
   //! for ex: A(2,3) = 5;
   SpMatVal& operator=(DataType v) {
      v_ = v;
      owner_->set(v_, i_, j_);
      return *this;
   }
private:
   owner_t*    owner_{nullptr};  //!< Pointer to owner SpMat. SpMatVal is just a view.
   DataType    v_{DataType{}};   //!< The value of the row-column pair (for speed)
   IndexType   i_{IndexType{}};  //!< The row
   IndexType   j_{IndexType{}};  //!< the column
};


//! enumerator for input matrix type.
enum class InputMatrix{ UNSPECIFIED, GENERATE, MTX };
//! enumerator for output handling
enum class OutputMode{ STD, FILE };

/*!
 * Session option for each invocation of the executable
 */
struct session_t {
   InputMatrix    inputMatrix    {InputMatrix::UNSPECIFIED};   //!< Source of the matrix
   std::ifstream  mtxFile        {};                           //!< matrix file in MatrixMarket format
   std::size_t    gen_size       {};                           //!< size of the matrix if we select to generate a random one
   double         gen_prob       {};                           //!< probability of the binomial distribution for the matrix
                                                               //!< if we generate one
   OutputMode     outputMode     {OutputMode::STD};            //!< Type of the output file
   std::ofstream  outFile        {};                           //!< File to use for output
   std::size_t    max_threads    {};                           //!< Maximum threads to use
   std::size_t    repeat         {1};                          //!< How many times we execute the calculations part
   bool           timing         {false};                      //!< Enable timing prints of the program
   bool           verbose        {false};                      //!< Flag to enable verbose output to stdout
#if CODE_VERSION == 3
   bool           makeSymmetric  {false};                      //!< symmetric matrix creation flag (true by default)
#else
   bool           makeSymmetric  {true};                       //!< symmetric matrix creation flag (true by default)
#endif
   bool           validate_mtx   {false};                      //!< Flag to request mtx input data triangular validation.
   bool           print_count    {false};                      //!< Flag to request total count printing
   bool           mtx_print      {false};                      //!< matrix print flag
   std::size_t    mtx_print_size {};                           //!< matrix print size
   bool           dynamic        {false};                      //!< Selects dynamic scheduling for OpenMP and pthreads.
};

extern session_t session;


#endif /* IMPL_HPP_ */