MatrixBackend.hpp

#pragma once
 
#include <list>
#include <string>
#include <memory>
 
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/matrixindexset.hh>
 
#include <amdis/Output.hpp>
#include <amdis/linearalgebra/SparsityPattern.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
#include <amdis/linearalgebra/istl/VectorBackend.hpp>
 
namespace AMDiS
{
  template <class T, class = void>
  struct BlockMatrixType
  {
    using type = Dune::FieldMatrix<T,1,1>;
  };
 
  template <class T>
  struct BlockMatrixType<T, typename T::field_type>
  {
    using type = T;
  };
 
  template <class T, class CommType>
  class ISTLBCRSMatrix
  {
  public:
    using BaseMatrix = Dune::BCRSMatrix<typename BlockMatrixType<T>::type>;
 
    using Comm = CommType;
 
    using value_type = typename BaseMatrix::block_type;
 
    using size_type = typename BaseMatrix::size_type;
 
  public:
    ISTLBCRSMatrix(Comm const& comm, Comm const&)
      : comm_(comm)
    {}
 
    BaseMatrix const& matrix() const
    {
      return matrix_;
    }
 
    BaseMatrix& matrix()
    {
      return matrix_;
    }
 
    Comm const& comm() const { return comm_; }
 
    void init(SparsityPattern const& pattern)
    {
      pattern.applyTo(matrix_);
      initialized_ = true;
    }
 
    void init()
    {
      matrix_ = value_type(0);
      initialized_ = true;
    }
 
    void finish()
    {
      initialized_ = false;
    }
 
 
    void insert(size_type r, size_type c, value_type const& value)
    {
      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
      test_exit_dbg( r < matrix_.N() && c < matrix_.M() ,
          "Indices out of range [0,{})x[0,{})", matrix_.N(), matrix_.M() );
      matrix_[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( initialized_, "Occupation pattern not initialized!");
      for (size_type i = 0; i < size_type(rows.size()); ++i)
        for (size_type j = 0; j < size_type(cols.size()); ++j)
          matrix_[rows[i]][cols[j]] += mat[i][j];
    }
 
    template <class RowInd>
    void zeroRows(RowInd const& rowInd, bool diag)
    {
      // loop over the matrix rows
      for (size_type i : rowInd) {
        auto cIt = matrix_[i].begin();
        auto cEndIt = matrix_[i].end();
        // loop over nonzero matrix entries in current row
        for (; cIt != cEndIt; ++cIt)
          *cIt = (diag && i == cIt.index() ? T(1) : T(0));
      }
    }
 
    template <class RowInd, class ColInd>
    void zeroRows(RowInd const& rowInd, ColInd const& colInd, bool diag)
    {
      // loop over the matrix rows
      auto rowIt = rowInd.begin();
      auto rowEnd = rowInd.end();
      auto colIt = colInd.begin();
      for (; rowIt != rowEnd; ++rowIt,++colIt) {
        auto cIt = matrix_[*rowIt].begin();
        auto cEndIt = matrix_[*rowIt].end();
        // loop over nonzero matrix entries in current row
        for (; cIt != cEndIt; ++cIt)
          *cIt = (diag && *colIt == cIt.index() ? T(1) : T(0));
      }
    }
 
    template <class RowInd, class ColInd>
    void zeroRowsColumnsImpl(RowInd const& rowInd, ColInd const& colInd, bool diag, ISTLBlockVector<T> const* x = nullptr, ISTLBlockVector<T>* b = nullptr)
    {
      assert(rowInd.size() == colInd.size());
 
      std::vector<bool> colBitVec(matrix_.M(), false);
      for (std::size_t i = 0; i < colInd.size(); ++i)
        colBitVec[colInd[i]] = true;
 
      // loop over nonzero matrix entries in current row
      for (size_type i : rowInd) {
        auto cIt = matrix_[i].begin();
        auto cEndIt = matrix_[i].end();
        for (; cIt != cEndIt; ++cIt)
          *cIt = T(0);
      }
 
      // iterate over the matrix
      for (size_type i = 0; i < matrix_.N(); ++i) {
        auto cIt = matrix_[i].begin();
        auto cEndIt = matrix_[i].end();
        // loop over nonzero matrix entries in current row
        for (; cIt != cEndIt; ++cIt) {
          if (colBitVec[cIt.index()]) {
            if (diag && x && b)
              b->vector()[i] -= (*cIt) * x->vector()[cIt.index()];
            *cIt = 0;
          }
        }
      }
 
      // set diagonal entry
      if (diag) {
        for (std::size_t i = 0; i < rowInd.size(); ++i) {
          test_exit_dbg(matrix_.exists(rowInd[i], colInd[i]),
            "The entry (",rowInd[i],",",colInd[i],") does not yet exist. Must be inserted in the pattern.");
          matrix_[rowInd[i]][colInd[i]] = 1;
        }
      }
    }
    template <class RowInd>
    void zeroRowsColumns(RowInd const& rowInd, bool diag)
    {
      zeroRowsColumnsImpl(rowInd,rowInd,diag,nullptr,nullptr);
    }
 
    template <class RowInd>
    void zeroRowsColumns(RowInd const& rowInd, bool diag, ISTLBlockVector<T> const& x, ISTLBlockVector<T>& b)
    {
      zeroRowsColumnsImpl(rowInd,rowInd,diag,&x,&b);
    }
 
    template <class RowInd, class ColInd>
    void zeroRowsColumns(RowInd const& rowInd, ColInd const& colInd, bool diag)
    {
      zeroRowsColumnsImpl(rowInd,colInd,diag,nullptr,nullptr);
    }
 
    template <class RowInd, class ColInd>
    void zeroRowsColumns(RowInd const& rowInd, ColInd const& colInd, bool diag, ISTLBlockVector<T> const& x, ISTLBlockVector<T>& b)
    {
      zeroRowsColumnsImpl(rowInd,colInd,diag,&x,&b);
    }
 
    std::size_t nnz() const
    {
      return matrix_.nonzeroes();
    }
 
  private:
    BaseMatrix matrix_;
    Comm const& comm_;
 
    bool initialized_ = false;
  };
 
} // end namespace AMDiS