mappedmatrix.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_MAPPEDMATRIX_HH
#define DUNE_FUFEM_MAPPEDMATRIX_HH
 
#include <dune/common/iteratorfacades.hh>
 
#include <dune/istl/matrixutils.hh>
 
 
template<class Map>
class MappedMatrixColIterator :
    public Dune::ForwardIteratorFacade<
        MappedMatrixColIterator<Map>,
        typename Map::block_type,
        const typename Map::block_type>
 
{
    typedef typename Map::Matrix Matrix;
    static const int rowFactor = Map::rowFactor;
    static const int colFactor = Map::colFactor;
 
public:
    typedef typename Map::block_type Value;
    typedef typename Map::block_type Reference;
 
    MappedMatrixColIterator(const Map& map, typename Matrix::row_type::ConstIterator&& it, int row, int virtualRow, int virtualCol=0):
        map_(&map),
        it_(it),
        row_(row),
        virtualRow_(virtualRow),
        virtualCol_(virtualCol)
    {}
 
    void increment()
    {
        if (virtualCol_<colFactor-1)
        {
            ++virtualCol_;
        }
        else
        {
            virtualCol_ = 0;
            ++it_;
        }
    }
 
    bool equals(const MappedMatrixColIterator& other) const
    {
        return (it_==other.it_) and (virtualRow_==other.virtualRow_) and (virtualCol_==other.virtualCol_);
    }
 
    const Reference dereference() const
    {
        return map_->apply(*it_, row_, it_.index(), virtualRow_, virtualCol_);
    }
 
    int index() const
    {
        return it_.index()*colFactor+virtualCol_;
    }
 
 
protected:
    const Map* map_;
    typename Matrix::row_type::ConstIterator it_;
    const int row_;
    int virtualRow_;
    int virtualCol_;
};
 
 
 
template<class Map>
class MappedMatrixRow
{
    typedef typename Map::Matrix Matrix;
    using size_type = typename Matrix::size_type;
    static const int rowFactor = Map::rowFactor;
    static const int colFactor = Map::colFactor;
public:
    typedef MappedMatrixColIterator<Map> ConstIterator;
 
    MappedMatrixRow(const Map& map, const Matrix& matrix, int row, int virtualRow):
        map_(&map),
        matrix_(&matrix),
        row_(row),
        virtualRow_(virtualRow)
    {}
 
    ConstIterator begin() const
    {
        return ConstIterator(*map_, (*matrix_)[row_].begin(), row_, virtualRow_);
    }
 
    ConstIterator end() const
    {
        return ConstIterator(*map_, (*matrix_)[row_].end(), row_, virtualRow_);
    }
 
    size_type size() const
    {
        return (*matrix_)[row_].size() * colFactor;
    }
 
    ConstIterator find(typename Matrix::size_type col) const
    {
        auto originalEntry = (*matrix_)[row_].find(col/colFactor);
        return (originalEntry==(*matrix_)[row_].end())
            ? end()
            : ConstIterator(*map_, (*matrix_)[row_].find(col/colFactor), row_, virtualRow_, col % colFactor);
    }
 
protected:
    const Map* map_;
    const Matrix* matrix_;
    int row_;
    int virtualRow_;
};
 
namespace Dune
{
  // Forward declaration, so MatrixDimension can be friend of MappedMatrix
  template<typename M>
  struct MatrixDimension;
}
 
 
template<class Map>
class MappedMatrix
{
    template <typename M>
    friend class MatrixDimension;
 
    typedef typename Map::Matrix Matrix;
    static const int rowFactor = Map::rowFactor;
    static const int colFactor = Map::colFactor;
 
public:
    using size_type = typename Matrix::size_type;
    typedef MappedMatrixRow<Map> row_type;
    typedef typename Map::block_type block_type;
    typedef MappedMatrixColIterator<Map> ConstColIterator;
 
    template<class T>
    class RealRowIterator
      : public Dune::RandomAccessIteratorFacade<RealRowIterator<T>, T, const T>
    {
    public:
      typedef typename std::remove_const<T>::type ValueType;
 
      RealRowIterator (const Map& map, const Matrix& matrix, typename Matrix::const_iterator originalRow)
        : map_(&map), matrix_(&matrix), originalRow_(originalRow), virtualRow_(0)
      {}
 
      RealRowIterator () = default;
 
      size_type index () const
      {
        return originalRow_.index() * rowFactor + virtualRow_;
      }
 
      std::ptrdiff_t distanceTo(const RealRowIterator<const ValueType>& other) const
      {
        return originalRow_.distanceTo(other.originalRow_) * rowFactor + other.virtualRow_ - virtualRow_;
      }
 
      bool equals (const RealRowIterator<const ValueType>& other) const
      {
        return originalRow_==other.originalRow_ and virtualRow_==other.virtualRow_;
      }
 
      void increment()
      {
        ++virtualRow_;
        if (virtualRow_ == rowFactor)
        {
          ++ originalRow_;
          virtualRow_ = 0;
        }
      }
 
      void decrement()
      {
        if (virtualRow_ == 0)
        {
          virtualRow_ = rowFactor - 1;
          --originalRow_;
        }
        else
          --virtualRow_;
      }
 
      void advance(std::ptrdiff_t diff)
      {
        originalRow_ += diff/rowFactor;
 
        if (diff > 0)
          for (std::ptrdiff_t i=0; i<diff%rowFactor; ++i)
            increment();
        else
          for (std::ptrdiff_t i=0; i<-diff%rowFactor; ++i)
            decrement();
      }
 
      row_type elementAt(std::ptrdiff_t diff) const
      {
        auto newOriginalRow = originalRow_.index() + diff/rowFactor;
        auto newVirtualRow = (virtualRow_ + diff) % rowFactor;
        return MappedMatrixRow<Map>(*map_, *matrix_, newOriginalRow, newVirtualRow);
      }
 
      row_type dereference () const
      {
        return MappedMatrixRow<Map>(*map_, *matrix_, originalRow_.index(), virtualRow_);
      }
 
      const Map* map_;
      const Matrix* matrix_;
      typename Matrix::const_iterator originalRow_;
      size_type virtualRow_;
    };
 
    using const_iterator = RealRowIterator<const row_type>;
 
    MappedMatrix(const Map& map, const Matrix& matrix):
        map_(&map),
        matrix_(&matrix)
    {}
 
    auto begin () const
    {
      return RealRowIterator<const row_type>(*map_, *matrix_, matrix_->begin());
    }
 
    auto end () const
    {
      return const_iterator(*map_, *matrix_, matrix_->end());
    }
 
    size_type N() const
    {
        return matrix_->N()*rowFactor;
    }
 
    size_type M() const
    {
        return matrix_->M()*colFactor;
    }
 
    const row_type operator[](int row) const
    {
        return MappedMatrixRow<Map>(*map_, *matrix_, row/rowFactor, row % rowFactor);
    }
 
    template<class X , class Y >
    void mv(const X &x, Y &y) const
    {
        int rows = this->N();
        y.resize(rows);
        for(int i=0; i<rows; ++i)
        {
            y[i] = 0.0;
            const row_type& Mi = (*this)[i];
            auto it = Mi.begin();
            auto end =  Mi.end();
            for(; it!=end; ++it)
                (*it).umv(x[it.index()], y[i]);
        }
    }
 
protected:
    const Map* map_;
    const Matrix* matrix_;
};
 
namespace Dune
{
  template<typename Map>
  struct MatrixDimension<MappedMatrix<Map> >
  {
    using Matrix = MappedMatrix<Map>;
    using size_type = typename Matrix::size_type;
 
    static size_type rowdim (const Matrix& /*A*/, size_type /*i*/)
    {
      return Matrix::block_type::rows;
    }
 
    static size_type coldim (const Matrix& /*A*/, size_type /*c*/)
    {
      return Matrix::block_type::cols;
    }
 
    static size_type rowdim (const Matrix& A) {
      return A.N()*Matrix::block_type::rows;
    }
 
    static size_type coldim (const Matrix& A) {
      return A.M()*Matrix::block_type::cols;
    }
  };
 
}  // namespace Dune
 
 
#endif // DUNE_FUFEM_MAPPEDMATRIX_HH