istlinterface.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2018 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef ISTLINTERFACE_HH
#define ISTLINTERFACE_HH
 
#include <memory>
#include <type_traits>
 
#include <boost/fusion/algorithm.hpp>
#include <boost/fusion/sequence.hpp>
 
#include <dune/istl/operators.hh>
#include "dune/istl/preconditioners.hh"
 
#include "fem/fixdune.hh"
#include "fem/functionspace.hh"
#include "fem/linearspace.hh"
 
#include "linalg/triplet.hh"
#include "linalg/crsutil.hh"
#include "linalg/factorization.hh"
#include "linalg/symmetricOperators.hh"
#include "linalg/threadedMatrix.hh"
 
#include "utilities/duneInterface.hh"
 
//---------------------------------------------------------------------
 
namespace Kaskade
{
  namespace IstlInterfaceDetail {
 
    template <class T>
    typename T::iterator begin(T& t)
    {
      return t.begin();
    }
 
    template <class T>
    typename T::const_iterator begin(T const& t)
    {
      return t.begin();
    }
 
    template <class Scalar>
    typename std::vector<Scalar>::iterator begin(std::vector<Scalar>& vec) { return vec.begin(); }
 
    template <class Scalar>
    typename std::vector<Scalar>::const_iterator begin(std::vector<Scalar> const& vec) { return vec.begin(); }
 
    using namespace boost::fusion;
 
 
    template <class ElementType, class SparseIdx, int rid, int cid, bool sym, bool trans>
    struct Block
    {
      static int const rowId = rid;          // row index of the block
      static int const colId = cid;          // col index of the block
      static bool const symmetric = sym;     // symmetry
      static bool const transposed = trans;  // whether only its mirror image is available in the assembler
 
      typedef ElementType BlockType;
      using SparseIndex = SparseIdx;
      typedef typename BlockType::field_type Scalar;
      typedef NumaBCRSMatrix<typename Transpose<BlockType,transposed>::type,SparseIndex> Matrix; // Matrix type in assembler
      typedef Dune::BlockVector<Dune::FieldVector<Scalar,BlockType::cols>> Domain; // take care of transposition
      typedef Dune::BlockVector<Dune::FieldVector<Scalar,BlockType::rows>> Range;
      typedef NumaBCRSMatrix<BlockType,SparseIndex> TMatrix;
 
      explicit Block(std::shared_ptr<Matrix const> const& m): matrix(m) {}
 
      std::shared_ptr<Matrix const> matrix;
      std::shared_ptr<TMatrix> threadedMatrix;
    };
 
    // Provides a heterogeneous sequence of all matrix blocks appearing
    // conceptually in the Galerkin operator, even if they are not stored
    // explicitly in the assembler because of symmetry.
    template <class GOP>
    class AllBlocks
    {
      struct Copy
      {
        template <class MBlock>
        auto operator()(MBlock const& mb) const {
          typedef std::remove_reference_t<MBlock> MB;
          typedef Block<typename MB::BlockType,typename MB::SparseIndex,MB::rowId,MB::colId,MB::symmetric,false> type;
          return type(mb.globalMatrixPointer());
        }
      };
 
      // Filters out flagged blocks with row index -1
      struct Cleanup {
        template <class MB> struct apply { typedef boost::mpl::bool_<MB::rowId>=0> type; };
      };
 
      // Blocks not to be mirrored are flagged with row index -1 to be filtered out later by Cleanup.
      struct Mirror
      {        
        template <class T> struct result {};
 
        template <class MBlock>
        struct result<Mirror(MBlock)> {
          typedef typename std::remove_reference<MBlock>::type MB;
          typedef typename Transpose<typename MB::BlockType,true>::type ElementType;
          typedef Block<ElementType, typename MB::SparseIndex,
              MB::mirror?MB::colId:-1, MB::rowId, MB::symmetric, true> type;
        };
 
        template <class MB>
        typename result<Mirror(MB)>::type operator()(MB const& mb) const {
          return typename result<Mirror(MB)>::type(mb.globalMatrixPointer());
        }
      };
 
      // It would be more natural to filter the mirror blocks before
      // actually mirroring them, but this does frequently result in empty
      // filter lists. For some reasons (don't know why), the filter view
      // produces wrong type results when empty lists occur.
 
      typedef typename result_of::transform<typename GOP::MatrixBlockArray const,Copy>::type CopiedBlocks;
      typedef typename result_of::transform<typename GOP::MatrixBlockArray const,Mirror>::type MirroredBlocks;
      typedef typename result_of::join<CopiedBlocks const,MirroredBlocks const>::type JointBlocks;
 
 
    public:
      typedef typename result_of::filter_if<JointBlocks const,Cleanup>::type const A;
      typedef typename result_of::as_vector<A>::type result;
 
      static result apply(GOP const& op) {
        return filter_if<Cleanup>(join(transform(op.getMatrix(),Copy()),transform(op.getMatrix(),Mirror())));
      }
    };
 
    // A boost fusion predicate that is true for some matrix block if it
    // lies within the half-open subrange given by the template
    // parameters.
    template <int firstRow, int lastRow, int firstCol, int lastCol>
    struct RangeBlockSelector
    {
      template <class Block> struct apply
      {
        typedef boost::mpl::bool_<firstRow<=Block::rowId && Block::rowId<lastRow && firstCol<=Block::colId && Block::colId<lastCol> type;
      };
    };
 
    // Shifts row and column indices of matrix blocks downwards and initializes the threaded matrix held by this block.
    template <int firstRow, int firstCol>
    struct Translate
    {
      // Constructor takes the number of threads onto which the threaded matrix shall be distributed.
      // As a conservative default, this is 1 (sequential).
      Translate(int nThreads_=1): nThreads(nThreads_)
      {
      }
 
      template <class T> struct result {};
 
      template <class Blk>
      struct result<Translate(Blk)> {
        typedef typename std::remove_reference<Blk>::type Bl;
        typedef Block<typename Bl::BlockType,typename Bl::SparseIndex,Bl::rowId-firstRow,Bl::colId-firstCol,Bl::symmetric,Bl::transposed> type;
      };
 
      template <class Bl>
      auto operator()(Bl const& bl) const 
      { 
        Block<typename Bl::BlockType,typename Bl::SparseIndex,Bl::rowId-firstRow,Bl::colId-firstCol,Bl::symmetric,Bl::transposed>  blk(bl.matrix);
        blk.threadedMatrix.reset(new typename Bl::TMatrix(*blk.matrix,bl.symmetric,bl.transposed,false));
        return blk;
      }
 
    private:
      int nThreads;
    };
 
 
    template <class FSElement, class Vector>
    typename std::enable_if<!std::is_same<FSElement,Vector>::value,void>::type
    toVector(FSElement const& x, Vector& coefficients)
    {
      coefficients.resize(x.dim());
      vectorToSequence(x,begin(coefficients));
    }
 
    template <class FSElement, class Vector>
    typename std::enable_if<std::is_same<FSElement,Vector>::value,void>::type
    toVector(FSElement const& x, Vector& coefficients)
    {
      coefficients = x;
    }
 
    template <class FSElement, class Vector>
    typename std::enable_if<!std::is_same<FSElement,Vector>::value,void>::type
    fromVector(Vector const& coefficients, FSElement& x)
    {
      assert(coefficients.size()>=x.dim());
      vectorFromSequence(x,begin(coefficients));
    }
 
    template <class FSElement, class Vector>
    typename std::enable_if<std::is_same<FSElement,Vector>::value,void>::type
    fromVector(Vector const& coefficients, FSElement& x)
    {
      x = coefficients;
    }
 
 
    template <class T>
    void toFile(std::vector<T> const& v, std::string const& filename, size_t precision = 10)
    {
      std::ofstream os;
      os.precision(precision);
      os.open(filename);
 
      for(T t : v) os << t << "\t";
      os << std::endl;
 
      os.close();
    }
 
 
    // Container storing information on blocks
    template <int firstRow_, int lastRow_, int firstCol_,int lastCol_>
    struct BlockInfo
    {
      static int const firstRow = firstRow_;
      static int const lastRow = lastRow_;
      static int const firstCol = firstCol_;
      static int const lastCol = lastCol_;
    };
  } // End of namespace IstlInterfaceDetail
 
  //---------------------------------------------------------------------
 
/*  template <class Operator>
  class TransposedOperator
  {
  public:
    typedef typename Operator::Domain Range;
    typedef typename Operator::Range Domain;
    typedef typename Operator::Scalar Scalar;
 
    explicit TransposedOperator(Operator const& A_) : A(A_) {}
 
    void apply(const Domain &x, Range &b) const { return A.applyTransposed(x,b); }
 
    void applyscaleadd(Scalar alpha, Domain const& x, Range& b) const { return A.applyscaleaddTransposed(alpha,x,b); }
 
    Operator const& getOperator() const { return A; }
 
  private:
    Operator const& A;
  };
*/
 
 
  template <class Matrix_, class Domain_, class Range_> class [[deprecated]] MatrixRepresentedOperator;
 
  template <class Scalar_, class SparseInt, class Domain_, class Range_>
  class [[deprecated]] MatrixRepresentedOperator<MatrixAsTriplet<Scalar_,SparseInt>,Domain_,Range_>
  : public Dune::LinearOperator<Domain_,Range_>
  {
  public:
    typedef MatrixAsTriplet<Scalar_,SparseInt> Matrix;
    typedef Matrix matrix_type; // for compatibility with Dune::MatrixAdapter
    typedef Domain_ Domain;
    typedef Range_ Range;
    typedef Scalar_ Scalar;
    
 
    MatrixRepresentedOperator() = default;
 
    explicit MatrixRepresentedOperator(Matrix const& A_) : A(A_)
    {}
 
    explicit MatrixRepresentedOperator(Matrix&& A_) : A(A_)
    {}
 
    MatrixAsTriplet<Scalar,SparseInt> getTriplet() const { return get<MatrixAsTriplet<Scalar,SparseInt> >(); }
 
    template <class OtherMatrix=Matrix>
    typename std::enable_if<std::is_same<Matrix,OtherMatrix>::value, Matrix const&>::type get() const { return A; }
 
    Matrix& get_non_const() { return A; }
 
 
    template <class OtherMatrix>
    typename std::enable_if<!std::is_same<Matrix,OtherMatrix>::value, OtherMatrix>::type get() const { return OtherMatrix(A); }
 
    void apply(Domain const& x, Range& b) const
    {
            std::vector<Scalar> tmpx, tmpb(A.nrows(),0.0);
            domainToVector(x,tmpx);
 
            A.mv(tmpx,tmpb);
 
            vectorToRange(tmpb,b);
    }
 
    void applyTransposed(Range const& x, Domain& b) const
    {
            std::vector<Scalar> tmpx, tmpb(A.ncols(),0.0);
            domainToVector(x,tmpx);
 
            A.mtv(tmpx,tmpb);
 
            vectorToRange(tmpb,b);
    }
 
    template <class OtherMatrix=Matrix> 
    typename std::enable_if<std::is_same<OtherMatrix,Matrix>::value,std::unique_ptr<Matrix> >::type getPointer() const
    {
      return std::unique_ptr<Matrix>(new Matrix(A));
    }
 
    template <class OtherMatrix> 
    typename std::enable_if<std::is_same<OtherMatrix,Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1> > >::value,std::unique_ptr<OtherMatrix> >::type getPointer() const
    {
      return A.toBCRS();
    }
 
    void applyscaleadd(Scalar alpha, Domain const& x, Range& b) const
    {
//      A.usmv(alpha,x,b);
            std::vector<Scalar> tmpx, tmpb;
            domainToVector(x,tmpx);
            rangeToVector(b,tmpb);
 
            A.usmv(alpha,tmpx, tmpb);
 
            vectorToRange(tmpb,b);
    }
 
    void applyscaleaddTransposed(Scalar alpha, Range const& x, Domain& b) const
    {
//      A.usmv(alpha,x,b);
            std::vector<Scalar> tmpx, tmpb;
            domainToVector(x,tmpx);
            rangeToVector(b,tmpb);
 
            A.usmtv(alpha,tmpx, tmpb);
 
            vectorToRange(tmpb,b);
    }
 
    template <class Vector> 
    void rangeToVector(Range const& y, Vector& coefficients) const
    {
      IstlInterfaceDetail::toVector(y,coefficients);
    }
 
 
    template <class Vector> 
    void domainToVector(Domain const& x, Vector& coefficients) const
    {
      IstlInterfaceDetail::toVector(x, coefficients);
    }
 
 
    template <class Vector>
    void vectorToDomain(Vector const& coefficients, Domain& x) const
    {
      IstlInterfaceDetail::fromVector(coefficients, x);
    }
 
 
    template <class Vector>
    void vectorToRange(Vector const& coefficients, Range& x) const
    {
      IstlInterfaceDetail::fromVector(coefficients, x);
    }
 
    MatrixRepresentedOperator& operator=(MatrixRepresentedOperator const& other) = default;
 
    MatrixRepresentedOperator& operator=(Matrix const& B)
    {
      A = B;   return *this;
    }
 
    MatrixRepresentedOperator& operator*= (Scalar alpha)
    {
      A *= alpha;   return *this;
    }
    
    MatrixRepresentedOperator& operator/= (Scalar alpha)
    {
      A /= alpha;   return *this;
    }
    
    MatrixRepresentedOperator& operator+= (MatrixRepresentedOperator const& B)
    {
      A += B;   return *this;
    }
    
    MatrixRepresentedOperator& operator+= (Matrix const& B)
    {
      A += B;   return *this;
    }
    
    MatrixRepresentedOperator& operator-= (MatrixRepresentedOperator const& B)
    {
      A -= B;   return *this;
    }
    
    MatrixRepresentedOperator& operator-= (Matrix const& B)
    {
      A -= B;   return *this;
    }
    
    virtual Dune::SolverCategory::Category category() const override
    {
      return Dune::SolverCategory::sequential;
    }
 
  private:
    Matrix A;
  };
 
  template <class MatrixBlock, class Allocator, class Domain_, class Range_>
  class [[deprecated]] MatrixRepresentedOperator<Dune::BCRSMatrix<MatrixBlock,Allocator>,Domain_,Range_>
  : public Dune::LinearOperator<Domain_,Range_>
  {
  public:
    typedef Dune::BCRSMatrix<MatrixBlock,Allocator> Matrix;
    typedef Domain_ Domain;
    typedef Range_ Range;
    typedef typename GetScalar<Domain>::type Scalar;
 
 
    MatrixRepresentedOperator() = default;
 
    explicit MatrixRepresentedOperator(Matrix const& A_) : A(A_)
    {}
 
    explicit MatrixRepresentedOperator(Matrix&& A_) : A(A_)
    {}
 
    MatrixAsTriplet<Scalar> getTriplet() const { return get<MatrixAsTriplet<Scalar> >(); }
 
    template <class OtherMatrix=Matrix>
    typename std::enable_if<std::is_same<Matrix,OtherMatrix>::value, Matrix const&>::type get() const { return A; }
 
    Matrix& get_non_const() { return A; }
 
 
    template <class OtherMatrix>
    typename std::enable_if<!std::is_same<Matrix,OtherMatrix>::value, OtherMatrix>::type get() const { return OtherMatrix(A); }
 
    void apply(Domain const& x, Range& b) const
    {
      A.mv(x,b);
    }
 
    void applyTransposed(Range const& x, Domain& b) const
    {
      A.mtv(x,b);
    }
 
    template <class OtherMatrix=Matrix>
    typename std::enable_if<std::is_same<OtherMatrix,Matrix>::value,std::unique_ptr<Matrix> >::type getPointer() const
    {
      return std::unique_ptr<Matrix>(new Matrix(A));
    }
 
    template <class OtherMatrix>
    typename std::enable_if<std::is_same<OtherMatrix,Dune::BCRSMatrix<Dune::FieldMatrix<Scalar,1,1> > >::value,std::unique_ptr<OtherMatrix> >::type getPointer() const
    {
      return A.toBCRS();
    }
 
    void applyscaleadd(Scalar alpha, Domain const& x, Range& b) const
    {
      A.usmv(alpha,x,b);
    }
 
    void applyscaleaddTransposed(Scalar alpha, Range const& x, Domain& b) const
    {
      A.usmtv(alpha,x,b);
    }
 
 
    template <class Vector> void rangeToVector(Range const& y, Vector& coefficients) const
    {
      IstlInterfaceDetail::toVector(y,coefficients);
    }
 
 
    template <class Vector> void domainToVector(Domain const& x, Vector& coefficients) const
    {
      IstlInterfaceDetail::toVector(x, coefficients);
    }
 
 
    template <class Vector>
    void vectorToDomain(Vector const& coefficients, Domain& x) const
    {
      IstlInterfaceDetail::fromVector(coefficients, x);
    }
 
 
    template <class Vector>
    void vectorToRange(Vector const& coefficients, Range& x) const
    {
      IstlInterfaceDetail::fromVector(coefficients, x);
    }
 
    MatrixRepresentedOperator& operator=(MatrixRepresentedOperator const& other) = default;
 
    MatrixRepresentedOperator& operator=(Matrix const& B)
    {
      A = B;   return *this;
    }
 
    MatrixRepresentedOperator& operator*= (Scalar alpha)
    {
      A *= alpha;   return *this;
    }
    
    MatrixRepresentedOperator& operator/= (Scalar alpha)
    {
      A /= alpha;   return *this;
    }
    
    MatrixRepresentedOperator& operator+= (MatrixRepresentedOperator const& B)
    {
      A += B;   return *this;
    }
    
    MatrixRepresentedOperator& operator+= (Matrix const& B)
    {
      A += B;   return *this;
    }
    
    MatrixRepresentedOperator& operator-= (MatrixRepresentedOperator const& B)
    {
      A -= B;   return *this;
    }
    
    MatrixRepresentedOperator& operator-= (Matrix const& B)
    {
      A -= B;   return *this;
    }
    
    virtual Dune::SolverCategory::Category category() const override
    {
      return Dune::SolverCategory::sequential;
    }
    
  private:
    Matrix A;
  };
 
  //-----------------------------------------------------------------------------------------------
 
  namespace IstlInterfaceDetail
  {
 
    // A traits class for determining the correct base class for an assembled Galerkin operator.
    // For variational functionals (i.e. symmetric problems), this is SymmetricLinearOperator,
    // otherwise Dune::LinearOperator.
    template <class Assembler,
              int firstRow, int lastRow,
              int firstCol, int lastCol,
              bool symmetric=false>
    struct Base
    {
      typedef typename Assembler::AnsatzVariableSetDescription::template CoefficientVectorRepresentation<firstCol,lastCol>::type Domain;
      typedef typename Assembler::TestVariableSetDescription::template CoefficientVectorRepresentation<firstRow,lastRow>::type   Range;
 
      static bool const inferredSymmetry
        =     (Assembler::Functional::type==VariationalFunctional   // a diagonally-centered block of a symmetric block matrix
                 && firstRow==firstCol && lastRow==lastCol)
              || (firstRow+1==lastRow && firstCol+1==lastCol && false); /* true if this one block is symmetric */  // TODO: implement the symmetry check
 
      using type = std::conditional_t<symmetric,
                                      SymmetricLinearOperator<Domain,Range>,
                                      Dune::LinearOperator<Domain,Range>>;
    };
  } // end of namespace IstlInterfaceDetail
 
 
 
  template <class Assembler_,
            int firstRow=0, int lastRow=Assembler_::TestVariableSetDescription::noOfVariables,
            int firstCol=0, int lastCol=Assembler_::AnsatzVariableSetDescription::noOfVariables,
            class BlockFilter=IstlInterfaceDetail::RangeBlockSelector<firstRow,lastRow,firstCol,lastCol>,
            bool symmetric=IstlInterfaceDetail::Base<Assembler_,firstRow,lastRow,firstCol,lastCol>::inferredSymmetry>
  class AssembledGalerkinOperator
  : public IstlInterfaceDetail::Base<Assembler_,firstRow,lastRow,firstCol,lastCol,symmetric>::type
  {
  public:
    typedef typename IstlInterfaceDetail::Base<Assembler_,firstRow,lastRow,firstCol,lastCol,symmetric>::type Base;
    typedef Assembler_ Assembler;
    using Domain = typename Assembler::AnsatzVariableSetDescription::template CoefficientVector<firstCol,lastCol>;
    using Range = typename Assembler::TestVariableSetDescription::template CoefficientVector<firstRow,lastRow>;
    typedef typename Assembler::Scalar Scalar;
 
    typedef IstlInterfaceDetail::BlockInfo<firstRow,lastRow,firstCol,lastCol> BlockInfo;
 
 
    explicit AssembledGalerkinOperator(Assembler const& assembler_,
                                       bool onlyLowerTriangle_=false, int nThreads_=0):
        onlyLowerTriangle(onlyLowerTriangle_),
        assembler(assembler_),
        blocks(boost::fusion::transform(boost::fusion::filter_if<BlockFilter>(IstlInterfaceDetail::AllBlocks<Assembler>::apply(assembler)),
                                        IstlInterfaceDetail::Translate<firstRow,firstCol>(nThreads_))),
        nThreads(nThreads_)
    { }
 
    virtual ~AssembledGalerkinOperator(){}
 
    void update()
    {
      using namespace boost::fusion;
      blocks = transform(filter_if<BlockFilter>(IstlInterfaceDetail::AllBlocks<Assembler>::apply(assembler)),
                         IstlInterfaceDetail::Translate<firstRow,firstCol>(nThreads));
    }
 
    virtual MatrixAsTriplet<Scalar> getTriplet() const __attribute__((deprecated))
    { return get<MatrixAsTriplet<Scalar>>(); }
 
    template <class Matrix> 
    Matrix get(int rowFirst=0, int rowLast=lastRow-firstRow, int colFirst=0, int colLast=lastCol-firstCol) const
    {
      return assembler.template get<Matrix>(onlyLowerTriangle, firstRow, lastRow, firstCol, lastCol);
    }
 
    template <class Matrix> 
    std::unique_ptr<Matrix> getPointer() const
    {
      return assembler.template getPointer<Matrix,firstRow,lastRow,firstCol,lastCol>(onlyLowerTriangle);
    }
 
    MatrixAsTriplet<Scalar> getmat() const __attribute__((deprecated));
 
    template <class Vector> 
    void rangeToVector(Range const& y, Vector& coefficients) const
    {
      IstlInterfaceDetail::toVector(y, coefficients);
    }
 
    template <class Vector> 
    void domainToVector(Domain const& x, Vector& coefficients) const
    {
      IstlInterfaceDetail::toVector(x, coefficients);
    }
 
 
    template <class Vector>
    void vectorToDomain(Vector const& coefficients, Domain& x) const
    {
      IstlInterfaceDetail::fromVector(coefficients, x);
    }
 
    template <class Vector>
    void vectorToRange(Vector const& coefficients, Range& x) const
    {
      IstlInterfaceDetail::fromVector(coefficients, x);
    }
 
    virtual void apply(Domain const& x, Range& b) const
    {
      assert(static_cast<void const*>(&x) != static_cast<void const*>(&b));
      boost::fusion::for_each(blocks,ApplyScaleAdd(assembler,1.0,x,b,true));
    }
 
    void applyTransposed(Range const& x, Domain& b) const
    {
      assert(static_cast<void const*>(&x) != static_cast<void const*>(&b));
      boost::fusion::for_each(blocks,TransposedApplyScaleAdd(assembler,1.0,x,b,true));
    }
 
    virtual Scalar applyDp(Domain const& x, Range& b) const
    {
      assert(static_cast<void const*>(&x) != static_cast<void const*>(&b));
      Scalar dp = 0;
      boost::fusion::for_each(blocks,ApplyScaleAdd(assembler,1.0,x,b,true,&dp));
      return symmetric? dp: 0;
    }
 
    virtual Scalar dp(Domain const& x, Range const& y) const
    {
      std::cerr << "AssembledGalerkinOperator::dp not yet implemented!\n"; abort();
      return 0; // TODO: fixme
    }
 
    virtual void applyscaleadd(Scalar alpha, Domain const& x, Range& b) const
    {
      assert(static_cast<void const*>(&x) != static_cast<void const*>(&b));
      boost::fusion::for_each(blocks,ApplyScaleAdd(assembler,alpha,x,b,false));
    }
 
    virtual void applyscaleaddTransposed(Scalar alpha, Range const& x, Domain& b) const
    {
      assert(static_cast<void const*>(&x) != static_cast<void const*>(&b));
      boost::fusion::for_each(blocks,TransposedApplyScaleAdd(assembler,alpha,x,b,false));
    }
 
    Assembler const& getAssembler() const 
    { 
      return assembler; 
    }
    
    template <int i>
    auto const& ansatzSpace() const 
    {
      static_assert(0<=i && i<lastCol-firstCol,"ansatz variable index out of range");
      int const varIdx = i-firstCol;
      int const spaceIdx = Assembler::AnsatzVariableSetDescription::template spaceIndex<i>;
      return *boost::fusion::at_c<spaceIdx>(assembler.spaces());
    }
 
    bool getOnlyLowerTriangle() const { return onlyLowerTriangle; }
 
    virtual Dune::SolverCategory::Category category() const override
    {
      return Dune::SolverCategory::sequential;
    }
 
  protected:
    typedef typename boost::fusion::result_of::as_vector<
        typename IstlInterfaceDetail::AllBlocks<Assembler>::result
        >::type allblocks;
 
    typedef typename boost::fusion::result_of::as_vector<
        typename boost::fusion::result_of::filter_if<allblocks,BlockFilter>::type
        >::type filtered;
 
    bool onlyLowerTriangle;
    Assembler const& assembler;
 
    typename boost::fusion::result_of::as_vector<
    typename boost::fusion::result_of::transform<filtered,
    IstlInterfaceDetail::Translate<firstRow,firstCol>>::type
    >::type blocks;
    bool isSymmetric;
 
 
    struct ApplyScaleAdd 
    {
      Assembler const& assembler;
      // Constructs the functor for matrix-vector multiplication and addition. If init is true, the first block
      // encountered in each row is called with apply (using alpha=1 and initializing the rhs to zero), otherwise
      // applyscaleadd is used. This prevents a double access of the right hand side, first for setting to zero,
      // then for adding the MV product.
      ApplyScaleAdd(Assembler const& as, Scalar a_, Domain const& x_, Range& y_, bool init, Scalar* dp_=0)
      : assembler(as)
      , a(a_)
      , x(x_)
      , y(y_)
      , dp(dp_? dp_: &dummy)
      {
        for (int i=0; i<lastRow; ++i)
          doInit[i] = init;
        *dp = 0;
      }
 
      template <class Block>
      void operator()(Block const& b) const
      {
        using namespace boost::fusion;
 
        auto const& xi = component<Block::colId>(x);
        auto      & yi = component<Block::rowId>(y);
        if (!b.threadedMatrix.get())
          abort();
 
        if (doInit[Block::rowId])
        {
          *dp += b.threadedMatrix->mv(xi,yi);
          doInit[Block::rowId] = false;
        }
        else
          *dp += b.threadedMatrix->usmv(a,xi,yi);      
      }
 
    private:
      Scalar a;
      Domain const& x;
      Range& y;
      Scalar* dp;
      Scalar dummy;                 // something to point to in case no target pointer has been provided
      mutable bool doInit[lastRow]; // ugly, but boost::fusion::for_each only works with const functors
    };
 
    struct TransposedApplyScaleAdd
    {
      Assembler const& assembler;
      // Constructs the functor for matrix-vector multiplication and addition. If init is true, the first block
      // encountered in each row is called with apply (using alpha=1 and initializing the rhs to zero), otherwise
      // applyscaleadd is used. This prevents a double access of the right hand side, first for setting to zero,
      // then for adding the MV product.
      TransposedApplyScaleAdd(Assembler const& as, Scalar a_, Range const& x_, Domain& y_, bool init, Scalar* dp_=0): assembler(as), a(a_), x(x_), y(y_), dp(dp_? dp_: &dummy)
      {
        for (int i=0; i<lastCol; ++i)
          doInit[i] = init;
        *dp = 0;
      }
 
      template <class Block>
      void operator()(Block const& b) const {
        using namespace boost::fusion;
 
        typename result_of::at_c<typename Range::Functions const,Block::rowId>::type xi = at_c<Block::rowId>(x.data);
        typename result_of::at_c<typename Domain::Functions,Block::colId>::type  yi = at_c<Block::colId>(y.data);
        if (!b.threadedMatrix.get())
          abort();
 
        if (doInit[Block::colId])
        {
          *dp += b.threadedMatrix->mv(xi,yi);
          doInit[Block::colId] = false;
        }
        else
          *dp += b.threadedMatrix->usmv(a,xi,yi);
 
        return;
      }
 
    private:
      Scalar a;
      Range const& x;
      Domain& y;
      Scalar* dp;
      Scalar dummy;                 // something to point to in case no target pointer has been provided
      mutable bool doInit[lastCol]; // ugly, but boost::fusion::for_each only works with const functors
    };
 
    int nThreads;
  };
  
  template <class Assembler,
            int firstRow=0, int lastRow=Assembler::TestVariableSetDescription::noOfVariables,
            int firstCol=0, int lastCol=Assembler::AnsatzVariableSetDescription::noOfVariables>
  AssembledGalerkinOperator<Assembler,firstRow,lastRow,firstCol,lastCol>
  makeAssembledGalerkinOperator(Assembler const& assembler)
  {
    return AssembledGalerkinOperator<Assembler,firstRow,lastRow,firstCol,lastCol>(assembler);
  }
  
  // -------------------------------------------------------------------------------------------------------
  // -------------------------------------------------------------------------------------------------------
 
 
  template <class NormalStepAssembler, class TangentialStepAssembler, int stateId=1, int controlId=0, int adjointId=2>
  class LagrangeOperator: public Dune::LinearOperator<typename NormalStepAssembler::AnsatzVariableSetDescription::template CoefficientVectorRepresentation<>::type,typename NormalStepAssembler::TestVariableSetDescription::template CoefficientVectorRepresentation<>::type>
  {
    typedef typename NormalStepAssembler::AnsatzVariableSetDescription::template CoefficientVectorRepresentation<>::type Domain;
    typedef typename NormalStepAssembler::TestVariableSetDescription::template CoefficientVectorRepresentation<>::type Range;
    static constexpr int firstRow = 0;
    static constexpr int firstCol = 0;
    static constexpr int lastRow = NormalStepAssembler::TestVariableSetDescription::noOfVariables;
    static constexpr int lastCol = NormalStepAssembler::AnsatzVariableSetDescription::noOfVariables;
    typedef IstlInterfaceDetail::RangeBlockSelector<firstRow,lastRow,firstCol,lastCol> NormalStepBlockFilter;
    typedef IstlInterfaceDetail::RangeBlockSelector<firstRow,adjointId,firstCol,adjointId> TangentialStepBlockFilter;
 
  public:
    typedef typename NormalStepAssembler::Scalar Scalar;
    typedef MatrixAsTriplet<Scalar> Triplet;
 
    typedef IstlInterfaceDetail::BlockInfo<firstRow,lastRow,firstCol,lastCol> BlockInfo;
 
    static int const category = Dune::SolverCategory::sequential;
 
    LagrangeOperator(NormalStepAssembler const& normalStepAssembler_, TangentialStepAssembler const& tangentialStepAssembler_,  bool onlyLowerTriangle_=false, int nThreads_=0):
      onlyLowerTriangle(onlyLowerTriangle_),
      normalStepAssembler(normalStepAssembler_), tangentialStepAssembler(tangentialStepAssembler_),
      normalStepBlocks(boost::fusion::transform(boost::fusion::filter_if<NormalStepBlockFilter>(IstlInterfaceDetail::AllBlocks<NormalStepAssembler>::apply(normalStepAssembler)),IstlInterfaceDetail::Translate<firstRow,firstCol>(nThreads_))),
      tangentialStepBlocks(boost::fusion::transform(boost::fusion::filter_if<TangentialStepBlockFilter>(IstlInterfaceDetail::AllBlocks<TangentialStepAssembler>::apply(tangentialStepAssembler)),IstlInterfaceDetail::Translate<firstRow,firstCol>(nThreads_))),
      nThreads(nThreads_)
    {}
 
    virtual ~LagrangeOperator(){}
 
    void update()
    {
      assert(false);
      //      blocks = boost::fusion::transform(boost::fusion::filter_if<BlockFilter>(IstlInterfaceDetail::AllBlocks<Assembler>::apply(assembler)),
      //          IstlInterfaceDetail::Translate<firstRow,firstCol>(nThreads));
    }
 
    virtual MatrixAsTriplet<Scalar> getTriplet() const
    {
      std::vector<size_t> offset(3,0);
      Triplet result, tmp;
      for(size_t i=0; i<3; ++i)
      {
        for(size_t j=0; j<3; ++j)
        {
          if(i==adjointId || j== adjointId || (i==controlId && j==controlId)) tmp = normalStepAssembler.template get<Triplet>(onlyLowerTriangle,i,i+1,j,j+1);
          else tmp = tangentialStepAssembler.template get<Triplet>(onlyLowerTriangle,i,i+1,j,j+1);
          tmp.shiftIndices(offset[i],offset[j]);
          result += tmp;
          if(i==0)
          {
            if(j==0) offset[1] = tmp.nrows();
            if(j==1) offset[2] = tmp.nrows();
          }
        }
      }
 
      return result;
    }
 
    //    /// Access matrix. Use if Matrix does support move-semantics.
    //    template <class Matrix> Matrix get() const
    //    {
    //      return assembler.template get<Matrix>(onlyLowerTriangle, firstRow, lastRow, firstCol, lastCol);
    //    }
 
    //    /// Access matrix via unique_ptr. Use if Matrix does not support move-semantics.
    //    template <class Matrix> std::unique_ptr<Matrix> getPointer() const
    //        {
    //      return assembler.template getPointer<Matrix,firstRow,lastRow,firstCol,lastCol>(onlyLowerTriangle);
    //  }
 
    //    MatrixAsTriplet<Scalar> getmat() const __attribute__((deprecated));
 
 
    template <class Vector> void rangeToVector(Range const& y, Vector& coefficients) const
    {
      IstlInterfaceDetail::toVector(y, coefficients);
    }
 
 
    template <class Vector> void domainToVector(Domain const& x, Vector& coefficients) const
    {
      IstlInterfaceDetail::toVector(x, coefficients);
    }
 
 
    template <class Vector>
    void vectorToDomain(Vector const& coefficients, Domain& x) const
    {
      IstlInterfaceDetail::fromVector(coefficients, x);
    }
 
    template <class Vector>
    void vectorToRange(Vector const& coefficients, Range& x) const
    {
      IstlInterfaceDetail::fromVector(coefficients, x);
    }
 
    virtual void apply(Domain const& x, Range& b) const
    {
      assert(static_cast<void const*>(&x) != static_cast<void const*>(&b));
      boost::fusion::for_each(normalStepBlocks,NormalStepApplyScaleAdd(normalStepAssembler,1.0,x,b,true));
      boost::fusion::for_each(tangentialStepBlocks,TangentialStepApplyScaleAdd(tangentialStepAssembler,1.0,x,b,false));
    }
 
    virtual Scalar applyDp(Domain const& x, Range& b) const
    {
      assert(static_cast<void const*>(&x) != static_cast<void const*>(&b));
      std::cerr << "AssembledGalerkinOperator::dp not yet implemented!\n"; abort();
      return 0;
    }
 
    virtual Scalar dp(Domain const& x, Range const& y) const
    {
      std::cerr << "AssembledGalerkinOperator::dp not yet implemented!\n"; abort();
      return 0; // TODO: fixme
    }
 
    virtual void applyscaleadd(Scalar alpha, Domain const& x, Range& b) const
    {
      assert(static_cast<void const*>(&x) != static_cast<void const*>(&b));
      boost::fusion::for_each(normalStepBlocks,NormalStepApplyScaleAdd(normalStepAssembler,1.0,x,b,false));
      boost::fusion::for_each(tangentialStepBlocks,TangentialStepApplyScaleAdd(tangentialStepAssembler,1.0,x,b,false));
    }
 
    //    Assembler const& getAssembler() const
    //    {
    //      return assembler;
    //    }
    //
    //    bool getOnlyLowerTriangle() const { return onlyLowerTriangle; }
 
  protected:
    typedef typename IstlInterfaceDetail::AllBlocks<NormalStepAssembler>::result AllNormalStepBlocksSeq;
    typedef typename IstlInterfaceDetail::AllBlocks<TangentialStepAssembler>::result AllTangentialStepBlocksSeq;
    typedef typename boost::fusion::result_of::filter_if<AllNormalStepBlocksSeq,NormalStepBlockFilter>::type FilteredNormalStepBlocksSeq;
    typedef typename boost::fusion::result_of::filter_if<AllTangentialStepBlocksSeq,TangentialStepBlockFilter>::type FilteredTangentialStepBlocksSeq;
 
    typedef typename boost::fusion::result_of::as_vector<AllNormalStepBlocksSeq>::type AllNormalStepBlocks;
    typedef typename boost::fusion::result_of::as_vector<AllTangentialStepBlocksSeq>::type AllTangentialStepBlocks;
    typedef typename boost::fusion::result_of::as_vector<FilteredNormalStepBlocksSeq>::type FilteredNormalStepBlocks;
    typedef typename boost::fusion::result_of::as_vector<FilteredTangentialStepBlocksSeq>::type FilteredTangentialStepBlocks;
 
    typedef typename boost::fusion::result_of::as_vector<typename boost::fusion::result_of::transform<FilteredNormalStepBlocks,IstlInterfaceDetail::Translate<firstRow,firstCol>>::type>::type NormalBlocks;
    typedef typename boost::fusion::result_of::as_vector<typename boost::fusion::result_of::transform<FilteredTangentialStepBlocks,IstlInterfaceDetail::Translate<0,0>>::type>::type TangentialBlocks;
 
    bool onlyLowerTriangle;
    NormalStepAssembler const& normalStepAssembler;
    TangentialStepAssembler const& tangentialStepAssembler;
 
    NormalBlocks  normalStepBlocks;
    TangentialBlocks tangentialStepBlocks;
 
    bool isSymmetric;
    int nThreads;
 
    struct NormalStepApplyScaleAdd
    {
      NormalStepAssembler const& assembler;
      // Constructs the functor for matrix-vector multiplication and addition. If init is true, the first block
      // encountered in each row is called with apply (using alpha=1 and initializing the rhs to zero), otherwise
      // applyscaleadd is used. This prevents a double access of the right hand side, first for setting to zero,
      // then for adding the MV product.
      NormalStepApplyScaleAdd(NormalStepAssembler const& as, Scalar a_, Domain const& x_, Range& y_, bool init, Scalar* dp_=0): assembler(as), a(a_), x(x_), y(y_), dp(dp_? dp_: &dummy)
      {
        for (int i=0; i<lastRow; ++i)
          doInit[i] = init;
        *dp = 0;
      }
 
      template <class Block>
      void operator()(Block const& b) const {
        if( ( (Block::rowId==stateId && Block::colId==stateId) || (Block::rowId==stateId && Block::colId==controlId) || (Block::colId==stateId && Block::rowId==controlId) ) ) return;
//        if( (Block::rowId==stateId && Block::colId==stateId) ) return;
        using namespace boost::fusion;
 
        typename result_of::at_c<typename Domain::Functions const,Block::colId>::type xi = at_c<Block::colId>(x.data);
        typename result_of::at_c<typename Range::Functions,Block::rowId>::type  yi = at_c<Block::rowId>(y.data);
        if (!b.threadedMatrix.get())
          abort();
 
        if (doInit[Block::rowId])
        {
          *dp += b.threadedMatrix->mv(xi,yi);
          doInit[Block::rowId] = false;
        }
        else
          *dp += b.threadedMatrix->usmv(a,xi,yi);
 
        return;
      }
 
    private:
      Scalar a;
      Domain const& x;
      Range& y;
      Scalar* dp;
      Scalar dummy;                 // something to point to in case no target pointer has been provided
      mutable bool doInit[lastRow]; // ugly, but boost::fusion::for_each only works with const functors
    };
 
 
    struct TangentialStepApplyScaleAdd
    {
      TangentialStepAssembler const& assembler;
      // Constructs the functor for matrix-vector multiplication and addition. If init is true, the first block
      // encountered in each row is called with apply (using alpha=1 and initializing the rhs to zero), otherwise
      // applyscaleadd is used. This prevents a double access of the right hand side, first for setting to zero,
      // then for adding the MV product.
      TangentialStepApplyScaleAdd(TangentialStepAssembler const& as, Scalar a_, Domain const& x_, Range& y_, bool init, Scalar* dp_=0): assembler(as), a(a_), x(x_), y(y_), dp(dp_? dp_: &dummy)
      {
        for (int i=0; i<lastRow; ++i)
          doInit[i] = init;
        *dp = 0;
      }
 
      template <class Block>
      void operator()(Block const& b) const {
        if( !( (Block::rowId==stateId && Block::colId==stateId) || (Block::rowId==stateId && Block::colId==controlId) || (Block::colId==stateId && Block::rowId==controlId) ) ) return;
        using namespace boost::fusion;
 
        typename result_of::at_c<typename Domain::Functions const,Block::colId>::type xi = at_c<Block::colId>(x.data);
        typename result_of::at_c<typename Range::Functions,Block::rowId>::type  yi = at_c<Block::rowId>(y.data);
        if (!b.threadedMatrix.get())
          abort();
 
        if (doInit[Block::rowId])
        {
          *dp += b.threadedMatrix->mv(xi,yi);
          doInit[Block::rowId] = false;
        }
        else
          *dp += b.threadedMatrix->usmv(a,xi,yi);
 
        return;
      }
 
    private:
      Scalar a;
      Domain const& x;
      Range& y;
      Scalar* dp;
      Scalar dummy;                 // something to point to in case no target pointer has been provided
      mutable bool doInit[lastRow]; // ugly, but boost::fusion::for_each only works with const functors
    };
  };
 
  template <class Operator>
  class TransposedOperator : public Dune::LinearOperator<typename Operator::Range,typename Operator::Domain>
  {
  public:
    typedef typename Operator::Scalar Scalar;
    typedef typename Operator::Range Domain;
    typedef typename Operator::Domain Range;
    typedef typename Operator::Assembler Assembler;
 
    explicit TransposedOperator(Operator const& A_) : A(A_.getTriplet().transpose()){}
 
    virtual ~TransposedOperator(){}
 
    virtual void apply(Domain const& x, Range& b) const
    {
      A.apply(x,b);
    }
 
    virtual void applyscaleadd(Scalar alpha, Domain const& x, Range& b) const
    {
      A.applyscaleadd(alpha,x,b);
    }
 
    MatrixAsTriplet<Scalar> getTriplet() const { return A.getTriplet(); }
 
  private:
    MatrixRepresentedOperator<MatrixAsTriplet<double>,Domain,Range> A;
  };
 
} // end of namespace Kaskade
 
#endif