istlmatrix.hh

#ifndef DUNE_FEM_ISTLMATRIXWRAPPER_HH
#define DUNE_FEM_ISTLMATRIXWRAPPER_HH
 
#if HAVE_DUNE_ISTL
 
//- system includes
#include <vector>
#include <iostream>
#include <fstream>
#include <algorithm>
#include <set>
#include <map>
#include <string>
 
//- Dune common includes
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
 
//- Dune istl includes
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/preconditioners.hh>
 
//- Dune fem includes
#include <dune/fem/function/blockvectorfunction.hh>
#include <dune/fem/operator/common/localmatrix.hh>
#include <dune/fem/operator/common/localmatrixwrapper.hh>
#include <dune/fem/operator/common/operator.hh>
#include <dune/fem/operator/common/stencil.hh>
#include <dune/fem/function/common/scalarproducts.hh>
#include <dune/fem/io/parameter.hh>
#include <dune/fem/operator/matrix/columnobject.hh>
#include <dune/fem/storage/objectstack.hh>
 
#include <dune/fem/operator/matrix/istlmatrixadapter.hh>
#include <dune/fem/operator/matrix/istlpreconditioner.hh>
#include <dune/fem/operator/matrix/functor.hh>
 
namespace Dune
{
 
  namespace Fem
  {
 
 
    template< class MatrixObject >
    class ISTLLocalMatrix;
 
    template <class DomainSpaceImp, class RangeSpaceImp,
              class DomainBlock = Dune::FieldVector< typename DomainSpaceImp :: RangeFieldType, DomainSpaceImp:: localBlockSize >,
              class RangeBlock  = Dune::FieldVector< typename RangeSpaceImp  :: RangeFieldType, RangeSpaceImp :: localBlockSize > >
    class ISTLMatrixObject;
 
    // --ISTLMatrixHandle
    template <class LittleBlockType, class RowDiscreteFunctionImp, class ColDiscreteFunctionImp = RowDiscreteFunctionImp>
    class ImprovedBCRSMatrix : public BCRSMatrix<LittleBlockType>
    {
        friend struct MatrixDimension<ImprovedBCRSMatrix>;
      public:
        typedef RowDiscreteFunctionImp RowDiscreteFunctionType;
        typedef ColDiscreteFunctionImp ColDiscreteFunctionType;
 
        typedef BCRSMatrix<LittleBlockType> BaseType;
        typedef BaseType MatrixBaseType;
        typedef typename BaseType :: RowIterator RowIteratorType ;
        typedef typename BaseType :: ColIterator ColIteratorType ;
 
        typedef ImprovedBCRSMatrix< LittleBlockType, RowDiscreteFunctionImp, ColDiscreteFunctionImp > ThisType;
 
        typedef typename BaseType :: size_type size_type;
 
        //===== type definitions and constants
 
        typedef typename BaseType::field_type field_type;
 
        typedef typename BaseType::block_type block_type;
 
        typedef typename BaseType:: allocator_type allocator_type;
 
        typedef typename BaseType :: row_type row_type;
 
        typedef typename BaseType :: ColIterator ColIterator;
 
        typedef typename BaseType :: ConstColIterator ConstColIterator;
 
        typedef typename BaseType :: RowIterator RowIterator;
 
        typedef typename BaseType :: ConstRowIterator ConstRowIterator;
 
        typedef typename ColDiscreteFunctionType :: DiscreteFunctionSpaceType RangeSpaceType;
 
        typedef typename RowDiscreteFunctionType :: DofStorageType  RowBlockVectorType;
 
        typedef typename ColDiscreteFunctionType :: DofStorageType  ColBlockVectorType;
 
        typedef typename RangeSpaceType :: GridType :: Traits :: Communication   CommunicationType ;
 
        typedef typename BaseType :: BuildMode BuildMode ;
 
      public:
        ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nnz, double overflowFraction) :
          BaseType (rows, cols, nnz, overflowFraction, BaseType::implicit)
        {}
 
        ImprovedBCRSMatrix(size_type rows, size_type cols, size_type nz = 0 ) :
          BaseType (rows, cols, BaseType :: row_wise)
        {}
 
        ImprovedBCRSMatrix( ) :
          BaseType ()
        {}
 
        ImprovedBCRSMatrix(const ImprovedBCRSMatrix& org) :
          BaseType(org)
        {}
 
        ConstRowIterator slicedBegin( const size_type row ) const
        {
          return ConstRowIterator( this->r, row );
        }
 
        ConstRowIterator slicedEnd( const size_type row ) const
        {
          return ConstRowIterator( this->r, row );
        }
 
        std::pair< size_type, size_type > sliceBeginEnd( const size_type thread, const size_type numThreads ) const
        {
          const size_type sliceSize  = this->N() / numThreads ;
          const size_type sliceBegin = thread * sliceSize ;
          const size_type sliceEnd   = ( thread == numThreads-1 ) ? this->N() : (sliceBegin + sliceSize) ;
          return std::make_pair( sliceBegin, sliceEnd );
        }
 
      protected:
        template <class X, class Y, bool isMV>
        void mvThreadedImpl( [[maybe_unused]] const field_type& alpha,
                             const X& x, Y& y, std::integral_constant<bool, isMV> ) const
        {
          auto doMV = [this, &alpha, &x, &y] ()
          {
            // this is simply to avoid warnings since
            // there exist no maybe_unused for lambda captures
            if constexpr ( isMV )
            {
              [[maybe_unused]] field_type a = alpha;
            }
 
            const auto slice = sliceBeginEnd( MPIManager::thread(), MPIManager::numThreads() );
            const ConstRowIterator endi = slicedEnd( slice.second );
            for (ConstRowIterator i = slicedBegin(slice.first); i!=endi; ++i)
            {
              if constexpr ( isMV )
                y[i.index()] = 0;
 
              ConstColIterator endj = (*i).end();
              for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
              {
                auto&& xj = Dune::Impl::asVector(x[j.index()]);
                auto&& yi = Dune::Impl::asVector(y[i.index()]);
                if constexpr ( isMV )
                  Dune::Impl::asMatrix(*j).umv(xj, yi);
                else
                  Dune::Impl::asMatrix(*j).usmv(alpha, xj, yi);
              }
            }
          };
 
          try {
            // execute in parallel
            MPIManager :: run ( doMV );
          }
          catch ( const SingleThreadModeError& e )
          {
            // serial version
            if constexpr ( isMV )
              this->mv( x, y );
            else
            {
              DUNE_THROW(SingleThreadModeError,"ImprovedBCRSMatrix::usmvThreaded: Cannot recover from threading error. Disable threading!");
              // this->usmv( alpha, x, y );
            }
          }
        }
 
      public:
        template <class X, class Y>
        void usmvThreaded( const field_type& alpha, const X& x, Y& y ) const
        {
          mvThreadedImpl(alpha, x, y, std::false_type() );
        }
 
        template <class X, class Y>
        void mvThreaded( const X& x, Y& y ) const
        {
          mvThreadedImpl(1.0, x, y, std::true_type() );
        }
 
        template < class SparsityPattern >
        void createEntries(const SparsityPattern& sparsityPattern)
        {
          const auto spEnd = sparsityPattern.end();
 
          // insert map of indices into matrix
          auto endcreate = this->createend();
          for(auto create = this->createbegin(); create != endcreate; ++create)
          {
            const auto row = sparsityPattern.find( create.index() );
            // if a row is empty then do nothing
            if( row == spEnd ) continue;
 
            // insert all indices for this row
            for ( const auto& col : row->second )
              create.insert( col );
          }
        }
 
        void clear()
        {
          for (auto& row : *this)
            for (auto& entry : row)
              entry = 0;
        }
 
        void unitRow( const size_t row )
        {
          block_type idBlock( 0 );
          for (int i = 0; i < idBlock.rows; ++i)
              idBlock[i][i] = 1.0;
 
          auto& matRow = (*this)[ row ];
          auto colIt = matRow.begin();
          const auto& colEndIt = matRow.end();
          for (; colIt != colEndIt; ++colIt)
          {
              if( colIt.index() == row )
                  *colIt = idBlock;
              else
                  *colIt = 0.0;
          }
        }
 
        template <class HangingNodesType>
        void setup(ThisType& oldMatrix, const HangingNodesType& hangingNodes)
        {
          // necessary because element traversal not necessarily is in ascending order
          typedef std::set< std::pair<int, block_type> > LocalEntryType;
          typedef std::map< int , LocalEntryType > EntriesType;
          EntriesType entries;
 
          // map of indices
          std::map< int , std::set<int> > indices;
          // not insert map of indices into matrix
          auto rowend  = oldMatrix.end();
          for(auto it  = oldMatrix.begin(); it != rowend; ++it)
          {
            const auto row = it.index();
            auto& localIndices = indices[ row ];
 
            if( hangingNodes.isHangingNode( row ) )
            {
              // insert columns into other columns
              const auto& cols = hangingNodes.associatedDofs( row );
              const auto colSize = cols.size();
              for(auto i=0; i<colSize; ++i)
              {
                assert( ! hangingNodes.isHangingNode( cols[i].first ) );
 
                // get local indices of col
                auto& localColIndices = indices[ cols[i].first ];
                auto& localEntry = entries[  cols[i].first ];
 
                // copy from old matrix
                auto endj = (*it).end();
                for (auto j= (*it).begin(); j!=endj; ++j)
                {
                  localColIndices.insert( j.index () );
                  localEntry.insert( std::make_pair( j.index(), (cols[i].second * (*j)) ));
                }
              }
 
              // insert diagonal and hanging columns
              localIndices.insert( row );
              for(auto i=0; i<colSize; ++i)
                localIndices.insert( cols[i].first );
            }
            else
            {
              // copy from old matrix
              auto endj = (*it).end();
              for (auto j= (*it).begin(); j!=endj; ++j)
                localIndices.insert( j.index () );
            }
          }
 
          // create matrix from entry map
          createEntries( indices );
 
          // not insert map of indices into matrix
          auto rowit  = oldMatrix.begin();
 
          auto endcreate = this->end();
          for(auto create = this->begin(); create != endcreate; ++create, ++rowit )
          {
            assert( rowit != oldMatrix.end() );
 
            const auto row = create.index();
            if( hangingNodes.isHangingNode( row ) )
            {
              const auto& cols = hangingNodes.associatedDofs( row );
 
              std::map< const int , block_type > colMap;
              // only working for block size 1 ath the moment
              assert( block_type :: rows == 1 );
              // insert columns into map
              const auto colSize = cols.size();
              for( auto i=0; i<colSize; ++i)
                colMap[ cols[i].first ] = -cols[i].second;
              // insert diagonal into map
              colMap[ row ] = 1;
 
              auto endj = (*create).end();
              for (auto j= (*create).begin(); j!=endj; ++j)
              {
                assert( colMap.find( j.index() ) != colMap.end() );
                (*j) = colMap[ j.index() ];
              }
            }
            // if entries are equal, just copy
            else if ( entries.find( row ) == entries.end() )
            {
              auto colit = (*rowit).begin();
              auto endj = (*create).end();
              for (auto j= (*create).begin(); j!=endj; ++j, ++colit )
              {
                assert( colit != (*rowit).end() );
                (*j) = (*colit);
              }
            }
            else
            {
              std::map< int , block_type > oldCols;
 
              {
                auto colend = (*rowit).end();
                for(auto colit = (*rowit).begin(); colit != colend; ++colit)
                  oldCols[ colit.index() ] = 0;
              }
 
              auto entry = entries.find( row );
              assert( entry  != entries.end ());
 
              {
                auto endcol = (*entry).second.end();
                for( auto co = (*entry).second.begin(); co != endcol; ++co)
                  oldCols[ (*co).first ] = 0;
              }
 
              {
                auto colend = (*rowit).end();
                for(auto colit = (*rowit).begin(); colit != colend; ++colit)
                  oldCols[ colit.index() ] += (*colit);
              }
 
              {
                auto endcol = (*entry).second.end();
                for( auto co = (*entry).second.begin(); co != endcol; ++co)
                  oldCols[ (*co).first ] += (*co).second;
              }
 
              auto endj = (*create).end();
              for (auto j= (*create).begin(); j!=endj; ++j )
              {
                auto colEntry = oldCols.find( j.index() );
                if( colEntry != oldCols.end() )
                  (*j) = (*colEntry).second;
                else
                  abort();
              }
            }
          }
        }
 
        void extractDiagonal( ColBlockVectorType& diag ) const
        {
          const auto endi = this->end();
          for (auto i = this->begin(); i!=endi; ++i)
          {
            // get diagonal entry of matrix
            const auto row = i.index();
            auto entry = (*i).find( row );
            const LittleBlockType& block = (*entry);
            enum { blockSize = LittleBlockType :: rows };
            for( auto l=0; l<blockSize; ++l )
              diag[ row ][ l ] = block[ l ][ l ];
          }
        }
 
        field_type operator()(const std::size_t row, const std::size_t col) const
        {
          const std::size_t blockRow(row/(LittleBlockType :: rows));
          const std::size_t localRowIdx(row%(LittleBlockType :: rows));
          const std::size_t blockCol(col/(LittleBlockType :: cols));
          const std::size_t localColIdx(col%(LittleBlockType :: cols));
 
          const auto& matrixRow(this->operator[](blockRow));
          auto entry = matrixRow.find( blockCol );
          const LittleBlockType& block = (*entry);
          return block[localRowIdx][localColIdx];
        }
 
        void set(const std::size_t row, const std::size_t col, field_type value)
        {
          const std::size_t blockRow(row/(LittleBlockType :: rows));
          const std::size_t localRowIdx(row%(LittleBlockType :: rows));
          const std::size_t blockCol(col/(LittleBlockType :: cols));
          const std::size_t localColIdx(col%(LittleBlockType :: cols));
 
          auto& matrixRow(this->operator[](blockRow));
          auto entry = matrixRow.find( blockCol );
          LittleBlockType& block = (*entry);
          block[localRowIdx][localColIdx] = value;
        }
 
        void print(std::ostream& s=std::cout, unsigned int offset=0) const
        {
          s.precision( 6 );
          const auto endi=this->end();
          for (auto i=this->begin(); i!=endi; ++i)
          {
            const auto endj = (*i).end();
            for (auto j=(*i).begin(); j!=endj; ++j)
              if( (*j).infinity_norm() > 1.e-15)
                s << i.index()+offset << " " << j.index()+offset << " " << *j << std::endl;
          }
        }
    };
 
 
 
    // ISTLLocalMatrixTraits
    // ---------------------
 
    template< class MatrixObject >
    struct ISTLLocalMatrixTraits
    {
      typedef typename MatrixObject::DomainSpaceType DomainSpaceType;
      typedef typename MatrixObject::RangeSpaceType RangeSpaceType;
      typedef typename DomainSpaceType::RangeFieldType RangeFieldType;
 
      typedef ISTLLocalMatrix< MatrixObject > LocalMatrixType;
      typedef typename MatrixObject::MatrixType::block_type LittleBlockType;
    };
 
 
    // ISTLLocalMatrix
    // ---------------
 
    template< class MatrixObject >
    class ISTLLocalMatrix
    : public LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > >
    {
    public:
      typedef LocalMatrixDefault< ISTLLocalMatrixTraits< MatrixObject > > BaseType;
 
      typedef MatrixObject MatrixObjectType;
      typedef typename MatrixObjectType::MatrixType MatrixType;
      typedef typename MatrixType::block_type LittleBlockType;
      // type of row and column indices
      typedef typename MatrixType :: size_type size_type;
 
      typedef typename MatrixObjectType::DomainSpaceType DomainSpaceType;
      typedef typename MatrixObjectType::RangeSpaceType RangeSpaceType;
 
      typedef typename MatrixObjectType::DomainEntityType DomainEntityType;
      typedef typename MatrixObjectType::RangeEntityType RangeEntityType;
 
      typedef typename DomainSpaceType::RangeFieldType DofType;
      typedef typename MatrixType::row_type RowType;
 
      typedef typename DomainSpaceType::BlockMapperType ColMapperType;
      typedef typename RangeSpaceType::BlockMapperType RowMapperType;
 
      static const int littleCols = MatrixObjectType::littleCols;
      static const int littleRows = MatrixObjectType::littleRows;
 
      typedef typename MatrixType::size_type Index;
 
      [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
      ISTLLocalMatrix ( const MatrixObjectType& mObj, const DomainSpaceType& domainSpace, const RangeSpaceType& rangeSpace )
        : BaseType( domainSpace, rangeSpace ),
          rowMapper_( rangeSpace.blockMapper() ),
          colMapper_( domainSpace.blockMapper() ),
          numRows_( rowMapper_.maxNumDofs() ),
          numCols_( colMapper_.maxNumDofs() ),
          matrixObj_( mObj )
      {}
 
      [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
      ISTLLocalMatrix ( const ISTLLocalMatrix& org )
        : BaseType( org ),
          rowMapper_(org.rowMapper_),
          colMapper_(org.colMapper_),
          numRows_( org.numRows_ ),
          numCols_( org.numCols_ ),
          matrices_(org.matrices_),
          matrixObj_(org.matrixObj_)
      {}
 
      void init ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
      {
        // initialize base functions sets
        BaseType :: init ( domainEntity, rangeEntity );
 
        numRows_  = rowMapper_.numDofs( rangeEntity );
        numCols_  = colMapper_.numDofs( domainEntity );
 
        // resize matrix pointer storage for new row/col numbers
        matrices_.resize( numRows_ );
        for( auto& row : matrices_ )
        {
          row.resize( numCols_, nullptr );
        }
 
        if(  matrixObj_.implicitModeActive() )
        {
          auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
          {
            return matrixObj_.matrix().entry( index.first, index.second );
          };
          initBlocks( blockAccess, domainEntity, rangeEntity );
        }
        else
        {
          auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
          {
            return matrixObj_.matrix()[ index.first ][ index.second ];
          };
          initBlocks( blockAccess, domainEntity, rangeEntity );
        }
      }
 
      template <class BlockAccess>
      void initBlocks( BlockAccess& blockAccess, const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity )
      {
        auto functor = [ this, &blockAccess ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
        {
          matrices_[ local.first ][ local.second ] = &blockAccess( index );
        };
 
        rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
      }
 
    private:
      // check whether given (row,col) pair is valid
      void check(int localRow, int localCol) const
      {
#ifndef NDEBUG
        const std::size_t row = (int) localRow / littleRows;
        const std::size_t col = (int) localCol / littleCols;
        const int lRow = localRow%littleRows;
        const int lCol = localCol%littleCols;
        assert( row < matrices_.size() ) ;
        assert( col < matrices_[row].size() );
        assert( lRow < littleRows );
        assert( lCol < littleCols );
#endif
      }
 
      DofType& getValue(const int localRow, const int localCol)
      {
        const int row = (int) localRow / littleRows;
        const int col = (int) localCol / littleCols;
        const int lRow = localRow%littleRows;
        const int lCol = localCol%littleCols;
        return (*matrices_[row][col])[lRow][lCol];
      }
 
    public:
      const DofType get(const int localRow, const int localCol) const
      {
        const int row = (int) localRow / littleRows;
        const int col = (int) localCol / littleCols;
        const int lRow = localRow%littleRows;
        const int lCol = localCol%littleCols;
        return (*matrices_[row][col])[lRow][lCol];
      }
 
      void scale (const DofType& scalar)
      {
        const size_type nrows = matrices_.size();
        for(size_type i=0; i<nrows; ++i)
        {
          const size_type ncols = matrices_[i].size();
          for(size_type j=0; j<ncols; ++j)
          {
            (*matrices_[i][j]) *= scalar;
          }
        }
      }
 
      void add(const int localRow, const int localCol , const DofType value)
      {
#ifndef NDEBUG
        check(localRow,localCol);
#endif
        getValue(localRow,localCol) += value;
      }
 
      void set(const int localRow, const int localCol , const DofType value)
      {
#ifndef NDEBUG
        check(localRow,localCol);
#endif
        getValue(localRow,localCol) = value;
      }
 
      void unitRow(const int localRow)
      {
        const int row = (int) localRow / littleRows;
        const int lRow = localRow%littleRows;
 
        // clear row
        doClearRow( row, lRow );
 
        // set diagonal entry to 1
        (*matrices_[row][row])[lRow][lRow] = 1;
      }
 
      void clear ()
      {
        const size_type nrows = matrices_.size();
        for(size_type i=0; i<nrows; ++i)
        {
          const size_type ncols = matrices_[i].size();
          for(size_type j=0; j<ncols; ++j)
          {
            (*matrices_[i][j]) = (DofType) 0;
          }
        }
      }
 
      void clearRow ( const int localRow )
      {
        const int row = (int) localRow / littleRows;
        const int lRow = localRow%littleRows;
 
        // clear the row
        doClearRow( row, lRow );
      }
 
      void resort ()
      {}
 
    protected:
      void doClearRow ( const int row, const int lRow )
      {
        // get number of columns
        const auto col = this->columns();
        for(auto localCol=0; localCol<col; ++localCol)
        {
          const int col = (int) localCol / littleCols;
          const int lCol = localCol%littleCols;
          (*matrices_[row][col])[lRow][lCol] = 0;
        }
      }
 
    private:
      // special mapper omitting block size
      const RowMapperType& rowMapper_;
      const ColMapperType& colMapper_;
 
      // number of local matrices
      int numRows_;
      int numCols_;
 
      // dynamic matrix with pointers to block matrices
      std::vector< std::vector< LittleBlockType* > > matrices_;
 
      // matrix to build
      const MatrixObjectType& matrixObj_;
    };
 
 
    template <class DomainSpaceImp, class RangeSpaceImp, class DomainBlock, class RangeBlock >
    class ISTLMatrixObject
    {
    public:
      typedef DomainSpaceImp DomainSpaceType;
      typedef RangeSpaceImp RangeSpaceType;
 
      typedef ISTLMatrixObject< DomainSpaceImp, RangeSpaceImp, DomainBlock, RangeBlock > ThisType;
      typedef ThisType  PreconditionMatrixType;
 
      typedef typename DomainSpaceType::GridType GridType;
 
      typedef typename RangeSpaceType  :: EntityType  RangeEntityType ;
      typedef typename DomainSpaceType :: EntityType  DomainEntityType ;
 
      enum { littleCols = DomainSpaceType :: localBlockSize };
      enum { littleRows = RangeSpaceType  :: localBlockSize };
 
      typedef FieldMatrix<typename DomainSpaceType :: RangeFieldType, littleRows, littleCols> LittleBlockType;
      typedef LittleBlockType  block_type;
      typedef LittleBlockType  MatrixBlockType;
 
      typedef ISTLBlockVectorDiscreteFunction< RangeSpaceType, RangeBlock >     RowDiscreteFunctionType;
      typedef ISTLBlockVectorDiscreteFunction< DomainSpaceType, DomainBlock >   ColumnDiscreteFunctionType;
 
      typedef typename RowDiscreteFunctionType :: DofStorageType    RowBlockVectorType;
      typedef typename ColumnDiscreteFunctionType :: DofStorageType ColumnBlockVectorType;
 
      typedef typename RangeSpaceType :: BlockMapperType  RowMapperType;
      typedef typename DomainSpaceType :: BlockMapperType ColMapperType;
 
      typedef ImprovedBCRSMatrix< LittleBlockType , ColumnDiscreteFunctionType , RowDiscreteFunctionType > MatrixType;
      typedef ISTLParallelMatrixAdapterInterface< MatrixType >    MatrixAdapterType;
 
      typedef ISTLLocalMatrix<ThisType> ObjectType;
      friend class ISTLLocalMatrix<ThisType>;
 
      typedef ThisType LocalMatrixFactoryType;
      typedef ObjectStack< LocalMatrixFactoryType > LocalMatrixStackType;
      typedef LocalMatrixWrapper< LocalMatrixStackType > LocalMatrixType;
      typedef ColumnObject< ThisType > LocalColumnObjectType;
 
    protected:
      const DomainSpaceType& domainSpace_;
      const RangeSpaceType&  rangeSpace_;
 
      // special row mapper
      RowMapperType& rowMapper_;
      // special col mapper
      ColMapperType& colMapper_;
 
      int size_;
      int sequence_;
 
      mutable std::unique_ptr< MatrixType > matrix_;
 
      mutable LocalMatrixStackType localMatrixStack_;
 
      mutable std::unique_ptr< MatrixAdapterType >     matrixAdap_;
      // overflow fraction for implicit build mode
      const double overflowFraction_;
      const bool threading_ ;
 
      ISTLSolverParameter param_;
    public:
      ISTLMatrixObject(const ISTLMatrixObject&) = delete;
 
      ISTLMatrixObject ( const DomainSpaceType &domainSpace, const RangeSpaceType &rangeSpace, const ISTLSolverParameter& param = ISTLSolverParameter() ) :
        domainSpace_(domainSpace)
        , rangeSpace_(rangeSpace)
        , rowMapper_( rangeSpace.blockMapper() )
        , colMapper_( domainSpace.blockMapper() )
        , size_(-1)
        , sequence_(-1)
        , localMatrixStack_( *this )
        , overflowFraction_( param.overflowFraction() )
        , threading_( param.threading() )
        , param_( param )
      {}
 
    protected:
      MatrixType& matrix() const
      {
        assert( matrix_ );
        return *matrix_;
      }
 
    public:
      bool threading() const { return threading_; }
 
      void printTexInfo(std::ostream& out) const
      {
        out << "ISTL MatrixObj: ";
        out  << "\\\\ \n";
      }
 
      std::string preconditionName() const
      {
        return "";
      }
 
      void createMatrixAdapter ( const ISTLSolverParameter& param ) const
      {
        if( !matrixAdap_ )
        {
          matrixAdap_ = ISTLMatrixAdapterFactory< ThisType > :: matrixAdapter( *this, param );
        }
        assert( matrixAdap_ );
      }
 
      MatrixAdapterType& matrixAdapter( const ISTLSolverParameter& parameter ) const
      {
        const ISTLSolverParameter* param = dynamic_cast< const ISTLSolverParameter* > (&parameter);
        if( param )
          createMatrixAdapter( *param );
        else
        {
          ISTLSolverParameter newparam( parameter );
          createMatrixAdapter( newparam );
        }
 
        finalizeAssembly();
        return *matrixAdap_;
      }
 
      MatrixAdapterType& matrixAdapter() const
      {
        return matrixAdapter( param_ );
      }
 
    public:
      MatrixType& exportMatrix() const
      {
        finalizeAssembly();
        return matrix();
      }
 
      bool implicitModeActive() const
      {
        // implicit build mode is only active when the
        // build mode of the matrix is implicit and the
        // matrix is currently being build
 
        if( matrix().buildMode() == MatrixType::implicit && matrix().buildStage() == MatrixType::building )
          return true;
        else
          return false;
      }
 
      void finalizeAssembly() const
      {
        // finalize build mode
        const_cast< ThisType& > (*this).compress();
      }
 
      // compress matrix if not already done before and only in implicit build mode
      void compress( )
      {
        if( implicitModeActive() )
        {
          matrix().compress();
        }
      }
 
      void clear()
      {
        matrix().clear();
        // clean matrix adapter and other helper classes
        removeObj();
      }
 
      void unitRow( const size_t row )
      {
        matrix().unitRow( row );
      }
 
    protected:
      template <class Container> // could be std::set or std::vector
      void setUnitRowImpl( const Container &rows, const double diagonal )
      {
        for (const auto& r : rows)
        {
          const std::size_t blockRow( r/(LittleBlockType :: rows) );
          const std::size_t localRowIdx( r%(LittleBlockType :: rows) );
          auto& row = matrix()[blockRow];
          const auto endcol = row.end();
#ifndef NDEBUG
          bool set = false;
#endif
          for (auto col=row.begin(); col!=endcol; ++col)
          {
            for (auto& entry : (*col)[localRowIdx])
              entry = 0;
            if (col.index() == blockRow)
            {
              (*col)[localRowIdx][localRowIdx] = diagonal;
#ifndef NDEBUG
              set = true;
#endif
            }
          }
          assert(set);
        }
      }
 
    public:
      template <class Container>
      void setUnitRows( const Container& unitRows, const Container& auxRows )
      {
        setUnitRowImpl( unitRows, 1.0 );
        setUnitRowImpl( auxRows,  0.0 );
      }
 
      void reserve()
      {
        reserve( Stencil<DomainSpaceType,RangeSpaceType>(domainSpace_, rangeSpace_), true );
      }
 
      template <class Set>
      void reserve (const std::vector< Set >& sparsityPattern )
      {
        reserve( StencilWrapper< DomainSpaceType,RangeSpaceType, Set >( sparsityPattern ), false );
      }
 
      template <class Stencil>
      void reserve(const Stencil &stencil, const bool proposeImplicit = true )
      {
        // if grid sequence number changed, rebuild matrix
        if(sequence_ != domainSpace().sequence())
        {
          removeObj();
 
          // do not use implicit build mode when multi threading is enabled
          const bool implicit = proposeImplicit && MPIManager::numThreads() == 1;
          if( implicit )
          {
            auto nnz = stencil.maxNonZerosEstimate();
            if( nnz == 0 )
            {
              if (domainSpace_.begin() != domainSpace_.end() && rangeSpace_.begin() != rangeSpace_.end())
              {
                Stencil tmpStencil( stencil );
                tmpStencil.fill( *(domainSpace_.begin()), *(rangeSpace_.begin()) );
                nnz = tmpStencil.maxNonZerosEstimate();
              }
            }
            matrix_.reset( new MatrixType( rowMapper_.size(), colMapper_.size(), nnz, overflowFraction_ ) );
          }
          else
          {
            matrix_.reset( new MatrixType( rowMapper_.size(), colMapper_.size() ) );
            matrix().createEntries( stencil.globalStencil() );
          }
 
          sequence_ = domainSpace().sequence();
        }
      }
 
      template <class HangingNodesType>
      void changeHangingNodes(const HangingNodesType& hangingNodes)
      {
        // create new matrix
        MatrixType* newMatrix = new MatrixType(rowMapper_.size(), colMapper_.size());
 
        // setup with hanging rows
        newMatrix->setup( *matrix_ , hangingNodes );
 
        // remove old matrix
        removeObj();
 
        // store new matrix
        matrix_.reset( newMatrix );
      }
 
      void extractDiagonal( ColumnDiscreteFunctionType& diag ) const
      {
        // extract diagonal entries
        matrix().extractDiagonal( diag.blockVector() );
      }
 
      bool rightPrecondition() const
      {
        return false;
      }
 
      void apply(const ColumnDiscreteFunctionType& arg, RowDiscreteFunctionType& dest) const
      {
        matrixAdapter().apply( arg.blockVector(), dest.blockVector() );
      }
 
      template <class CDF, class RDF>
      void apply(const CDF& arg, RDF& dest) const
      {
        // this can happen when different storages are used for discrete
        // function and matrix/solver objects
        DUNE_THROW(NotImplemented,"ISTLMatrixObj::apply called for DiscreteFunctions not specified in the template list");
      }
 
      void resort()
      {}
 
      void createPreconditionMatrix()
      {}
 
      void print(std::ostream & s) const
      {
        matrix().print(s);
      }
 
      const DomainSpaceType& domainSpace() const
      {
        return domainSpace_;
      }
      const RangeSpaceType&  rangeSpace() const
      {
        return rangeSpace_;
      }
 
      const RowMapperType& rowMapper() const
      {
        return rowMapper_;
      }
      const ColMapperType& colMapper() const
      {
        return colMapper_;
      }
 
      ObjectType* newObject() const
      {
        return new ObjectType(*this, domainSpace(), rangeSpace());
      }
 
      [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
      LocalMatrixType localMatrix( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity ) const
      {
        return LocalMatrixType( localMatrixStack_, domainEntity, rangeEntity );
      }
 
      [[deprecated("Use TemporaryLocalMatrix,LocalContribution and {get,add,set}LocalMatrix")]]
      LocalMatrixType localMatrix() const
      {
        return LocalMatrixType( localMatrixStack_ );
      }
 
      LocalColumnObjectType localColumn( const DomainEntityType &domainEntity ) const
      {
        return LocalColumnObjectType ( *this, domainEntity );
      }
 
    protected:
      template< class LocalBlock, class Operation >
      void applyToBlock ( const size_t row, const size_t col,
                          const LocalBlock &localBlock,
                          Operation& operation )
      {
        LittleBlockType& block = ( implicitModeActive() ) ? matrix().entry( row, col ) : matrix()[ row ][ col ];
        for( int i  = 0; i < littleRows; ++i )
          for( int j = 0; j < littleCols; ++j )
            operation( block[ i ][ j ], localBlock[ i ][ j ] );
      }
 
      template< class LocalBlock >
      void setBlock ( const size_t row, const size_t col,
                      const LocalBlock &localBlock )
      {
        typedef typename DomainSpaceType :: RangeFieldType Field;
        auto copy = [] ( Field& a, const typename LocalBlock::field_type& b ) { a = b; };
        applyToBlock( row, col, localBlock, copy );
      }
 
      template< class LocalBlock >
      void addBlock ( const size_t row, const size_t col,
                      const LocalBlock &localBlock )
      {
        typedef typename DomainSpaceType :: RangeFieldType Field;
        auto add = [] ( Field& a, const typename LocalBlock::field_type& b ) { a += b; };
        applyToBlock( row, col, localBlock, add );
      }
 
    public:
      template< class LocalMatrix, class Operation >
      void applyToLocalMatrix ( const DomainEntityType &domainEntity,
                                const RangeEntityType &rangeEntity,
                                const LocalMatrix &localMat,
                                Operation& operation )
      {
        typedef typename MatrixType::size_type Index;
        if( implicitModeActive() )
        {
          auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
          {
            return matrix().entry( index.first, index.second );
          };
          auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
          {
            LittleBlockType& block = blockAccess( index );
            for( int i  = 0; i < littleRows; ++i )
              for( int j = 0; j < littleCols; ++j )
                operation( block[ i ][ j ], localMat.get( local.first * littleRows + i, local.second *littleCols + j ) );
          };
          rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
        }
        else
        {
          auto blockAccess = [ this ] ( const std::pair< Index, Index > &index ) -> LittleBlockType&
          {
            return matrix()[ index.first][ index.second ];
          };
          auto functor = [ &localMat, blockAccess, operation ] ( std::pair< int, int > local, const std::pair< Index, Index > &index )
          {
            LittleBlockType& block = blockAccess( index );
            for( int i  = 0; i < littleRows; ++i )
              for( int j = 0; j < littleCols; ++j )
                operation( block[ i ][ j ], localMat.get( local.first * littleRows + i, local.second *littleCols + j ) );
          };
          rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
        }
      }
 
      template< class LocalMatrix >
      void setLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
      {
        typedef typename DomainSpaceType :: RangeFieldType  RangeFieldType;
        auto operation = [] ( RangeFieldType& a, const RangeFieldType& b )
        {
          a = b;
        };
        applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
      }
 
      template< class LocalMatrix >
      void addLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat )
      {
        typedef typename DomainSpaceType :: RangeFieldType  RangeFieldType;
        auto operation = [] ( RangeFieldType& a, const RangeFieldType& b )
        {
          a += b;
        };
        applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
      }
 
      template< class LocalMatrix, class Scalar >
      void addScaledLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, const LocalMatrix &localMat, const Scalar &s  )
      {
        typedef typename DomainSpaceType :: RangeFieldType  RangeFieldType;
        auto operation = [ &s ] ( RangeFieldType& a, const RangeFieldType& b )
        {
          a += s * b;
        };
        applyToLocalMatrix( domainEntity, rangeEntity, localMat, operation );
      }
 
      template< class LocalMatrix >
      void getLocalMatrix ( const DomainEntityType &domainEntity, const RangeEntityType &rangeEntity, LocalMatrix &localMat ) const
      {
        // make sure that matrix is in compressed state if build mode is implicit
        finalizeAssembly();
 
        typedef typename MatrixType::size_type Index;
        auto functor = [ &localMat, this ] ( std::pair< int, int > local, const std::pair< Index, Index > &global )
        {
          for( std::size_t i  = 0; i < littleRows; ++i )
            for( std::size_t j = 0; j < littleCols; ++j )
              localMat.set( local.first * littleRows + i, local.second *littleCols + j,
                matrix()[ global.first ][ global.second ][i][j] );
        };
 
        rowMapper_.mapEach( rangeEntity, makePairFunctor( colMapper_, domainEntity, functor ) );
      }
 
    protected:
      void preConErrorMsg(int preCon) const
      {
        exit(1);
      }
 
      void removeObj ()
      {
        matrixAdap_.reset( nullptr );
      }
    };
 
  } // namespace Fem
 
} // namespace Dune
 
#endif // #if HAVE_DUNE_ISTL
 
#endif // #ifndef DUNE_FEM_ISTLMATRIXWRAPPER_HH