MatrixBackend.hpp

#pragma once

#include <list>
#include <memory>
#include <string>
#include <vector>

#include <boost/numeric/mtl/matrix/compressed2D.hpp>
#include <boost/numeric/mtl/matrix/inserter.hpp>
#include <boost/numeric/mtl/utility/property_map.hpp>
#include <boost/numeric/mtl/utility/range_wrapper.hpp>

#include <amdis/Output.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
#include <amdis/linearalgebra/mtl/SlotSize.hpp>

namespace AMDiS
{
  template <class T>
  class MTLSparseMatrix
  {
  public:
    using BaseMatrix = mtl::compressed2D<T>;

    using value_type = typename BaseMatrix::value_type;

    using size_type = typename BaseMatrix::size_type;

  private:
    using Inserter = mtl::mat::inserter<BaseMatrix, mtl::operations::update_plus<value_type>>;

    using Pattern = SlotSize;

  public:
    template <class Basis>
    explicit MTLSparseMatrix(Basis const& /*rowBasis*/, Basis const& /*colBasis*/) {}

    BaseMatrix& matrix()
    {
      assert( !inserter_ );
      return matrix_;
    }

    BaseMatrix const& matrix() const
    {
      assert( !inserter_ );
      return matrix_;
    }

    void init(Pattern const& pattern)
    {
      test_exit(!inserter_, "Matrix already in insertion mode!");

      std::size_t slotSize = nnz() > 0 ? 6*nnz() / (5*num_rows(matrix_))
                                       : pattern.rowSizeEstimate();
      matrix_.change_dim(pattern.rows(), pattern.cols());
      set_to_zero(matrix_);

      symmetry_ = pattern.symmetry();
      inserter_ = new Inserter(matrix_, slotSize);
    }

    void init()
    {
      test_exit(!inserter_, "Matrix already in insertion mode!");

      std::size_t slotSize = 6*nnz() / (5*num_rows(matrix_));
      set_to_zero(matrix_);

      inserter_ = new Inserter(matrix_, slotSize);
    }


    void finish()
    {
      delete inserter_;
      inserter_ = nullptr;
    }


    void insert(size_type r, size_type c, value_type const& value)
    {
      test_exit_dbg(inserter_, "Inserter not initilized!");
      test_exit_dbg(r < num_rows(matrix_) && c < num_cols(matrix_),
          "Indices out of range [0,{})x[0,{})", num_rows(matrix_), num_cols(matrix_));
      if (value != value_type(0) || r == c)
        (*inserter_)[r][c] += value;
    }

    template <class Ind, class LocalMat>
    void scatter(Ind const& idx, LocalMat const& mat)
    {
      scatter(idx, idx, mat);
    }

    template <class RowInd, class ColInd, class LocalMat>
    void scatter(RowInd const& rows, ColInd const& cols, LocalMat const& mat)
    {
      test_exit_dbg(inserter_, "Inserter not initilized!");
      for (size_type i = 0; i < size_type(rows.size()); ++i)
        for (size_type j = 0; j < size_type(cols.size()); ++j)
          if (mat[i][j] != value_type(0) || i == j)
            (*inserter_)[rows[i]][cols[j]] += mat[i][j];
    }

    std::size_t nnz() const
    {
      return matrix_.nnz();
    }

    SymmetryStructure symmetry() const
    {
      return symmetry_;
    }

  private:
    BaseMatrix matrix_;

    Inserter* inserter_ = nullptr;

    SymmetryStructure symmetry_ = SymmetryStructure::unknown;
  };

} // end namespace AMDiS