diff --git a/examples/convection_diffusion.cc b/examples/convection_diffusion.cc
index 40871822f4fe7a9483e404fe2ab8bb8f38dc6124..14bea3a42cbf8d4e46213281d25ec3d3236cbbb2 100644
--- a/examples/convection_diffusion.cc
+++ b/examples/convection_diffusion.cc
@@ -1,6 +1,6 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
 #include <iostream>
 
 #include <amdis/AMDiS.hpp>
@@ -26,13 +26,13 @@ int main(int argc, char** argv)
 
   // -div(A*grad(u)) + div(b*u) + c*u = f
   auto opCD = convectionDiffusion(/*A=*/1.0, /*b=*/0.0, /*c=*/1.0, /*f=*/1.0);
-  prob.addMatrixOperator(opCD, _0, _0);
-  prob.addVectorOperator(opCD, _0);
+  prob.addMatrixOperator(opCD, 0, 0);
+  prob.addVectorOperator(opCD, 0);
 
   // set boundary condition
   auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
   auto dbcValues = [](auto const& x){ return 0.0; }; // set value
-  prob.addDirichletBC(predicate, _0, _0, dbcValues);
+  prob.addDirichletBC(predicate, 0, 0, dbcValues);
 
   AdaptInfo adaptInfo("adapt");
   prob.buildAfterCoarsen(adaptInfo, Flag(0));
diff --git a/examples/ellipt.cc b/examples/ellipt.cc
index 396737c3643a17662486482fe66a5fc45f540577..fdfc81c6a6bd1f85e41aa29816f8e4bfe2e5d65d 100644
--- a/examples/ellipt.cc
+++ b/examples/ellipt.cc
@@ -1,6 +1,6 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
 #include <iostream>
 
 #include <fmt/core.h>
diff --git a/examples/heat.cc b/examples/heat.cc
index d232264159f3ca35290fbaf364d9b3c06b032400..c2fa5e08d1c7d8afbd2fa611393bed74afa48822 100644
--- a/examples/heat.cc
+++ b/examples/heat.cc
@@ -1,6 +1,6 @@
-// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
-// vi: set et ts=4 sw=2 sts=2:
-
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
 #include <iostream>
 
 #include <amdis/AMDiS.hpp>
diff --git a/examples/navier_stokes.cc b/examples/navier_stokes.cc
index 4fcf941ba26aa0d16f10dc0ef26832dc01a881d7..762b369deba45fa66ffe002b519cb6bebf03a7ca 100644
--- a/examples/navier_stokes.cc
+++ b/examples/navier_stokes.cc
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
 #include <iostream>
 #include <ctime>
 #include <cmath>
@@ -41,8 +45,8 @@ int main(int argc, char** argv)
   AdaptInfo adaptInfo("adapt");
 
   // tree-paths for components
-  auto _v = 0_c;
-  auto _p = 1_c;
+  auto _v = Dune::Indices::_0;
+  auto _p = Dune::Indices::_1;
 
   // <1/tau * u, v>
   auto opTime = makeOperator(tag::testvec_trialvec{}, density);
diff --git a/examples/stokes0.cc b/examples/stokes0.cc
index c74c88b5dea8427ec2082d99b68e405a959d276c..e3785be1aa1dab3ed9b9cad3622b90a73b5fdedd 100644
--- a/examples/stokes0.cc
+++ b/examples/stokes0.cc
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
 #include <iostream>
 #include <ctime>
 #include <cmath>
@@ -30,8 +34,8 @@ int main(int argc, char** argv)
     Parameters::get("stokes->viscosity", viscosity);
 
     // tree-paths for components
-    auto _v = 0_c;
-    auto _p = 1_c;
+    auto _v = Dune::Indices::_0;
+    auto _p = Dune::Indices::_1;
 
     // <viscosity*grad(u_i), grad(v_i)>
     for (std::size_t i = 0; i < DOW; ++i) {
diff --git a/examples/stokes1.cc b/examples/stokes1.cc
index 8ae36a30e14f6f4a8ef60f215cdb9af81553ffe1..9126c3d406d904c9b9d65cbd59ff3a8213195af0 100644
--- a/examples/stokes1.cc
+++ b/examples/stokes1.cc
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
 #include <iostream>
 #include <ctime>
 #include <cmath>
@@ -30,8 +34,8 @@ int main(int argc, char** argv)
     Parameters::get("stokes->viscosity", viscosity);
 
     // tree-paths for components
-    auto _v = 0_c;
-    auto _p = 1_c;
+    auto _v = Dune::Indices::_0;
+    auto _p = Dune::Indices::_1;
 
     // <viscosity*grad(u_i), grad(v_i)>
     for (std::size_t i = 0; i < DOW; ++i) {
diff --git a/examples/stokes3.cc b/examples/stokes3.cc
index 54c75d0da556897759920a782f3cf0b029875902..44ccdd6a4db128ab582a7907dd4296ad7d7172c8 100644
--- a/examples/stokes3.cc
+++ b/examples/stokes3.cc
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
 #include <iostream>
 #include <ctime>
 #include <cmath>
@@ -24,8 +28,8 @@ int main(int argc, char** argv)
   Parameters::get("stokes->viscosity", viscosity);
 
   // tree-paths for components
-  auto _v = 0_c;
-  auto _p = 1_c;
+  auto _v = Dune::Indices::_0;
+  auto _p = Dune::Indices::_1;
 
   auto opStokes = makeOperator(tag::stokes{}, viscosity);
   prob.addMatrixOperator(opStokes, treepath(), treepath());
diff --git a/examples/vecellipt.cc b/examples/vecellipt.cc
index ea4e9f0541ffbf562448ac42c93e4aaf8adc946b..f078cdd724cc2f88102bcc7ed5c3a005f19e38e9 100644
--- a/examples/vecellipt.cc
+++ b/examples/vecellipt.cc
@@ -1,3 +1,7 @@
+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
 #include <iostream>
 
 #include <amdis/AMDiS.hpp>
diff --git a/src/amdis/Assembler.hpp b/src/amdis/Assembler.hpp
index ef99fb52f1cd519f54b83db8a0e34dfd3c7c219c..8b43adaffeff0a9d9723d661d76e42973ffb62dd 100644
--- a/src/amdis/Assembler.hpp
+++ b/src/amdis/Assembler.hpp
@@ -10,7 +10,6 @@
 #include <amdis/LinearAlgebra.hpp>
 #include <amdis/LocalAssemblerList.hpp>
 #include <amdis/common/Mpl.hpp>
-#include <amdis/common/TypeDefs.hpp>
 
 namespace AMDiS
 {
diff --git a/src/amdis/Assembler.inc.hpp b/src/amdis/Assembler.inc.hpp
index 60fb97c207fee0c00b0cb28833ef9e986d8d43d2..2d4610a827bdaf5b31b120481ae8cb5add4dc52c 100644
--- a/src/amdis/Assembler.inc.hpp
+++ b/src/amdis/Assembler.inc.hpp
@@ -24,8 +24,8 @@ void Assembler<Traits>::assemble(
   // 2. create a local matrix and vector
   std::size_t localSize = localView.maxSize();
 
-  Impl::ElementMatrix elementMatrix(localSize, localSize);
-  Impl::ElementVector elementVector(localSize);
+  mtl::mat::dense2D<typename SystemMatrixType::value_type> elementMatrix(localSize, localSize);
+  mtl::vec::dense_vector<typename SystemVectorType::value_type> elementVector(localSize);
 
   // 3. traverse grid and assemble operators on the elements
   for (auto const& element : elements(globalBasis_.gridView()))
diff --git a/src/amdis/DirichletBC.hpp b/src/amdis/DirichletBC.hpp
index 14267f6c5c1288113d060aa1b700862e7eb394e0..30960d343be66ef19a143304054ef5847f032d0c 100644
--- a/src/amdis/DirichletBC.hpp
+++ b/src/amdis/DirichletBC.hpp
@@ -11,12 +11,14 @@
 #include <amdis/Output.hpp>
 #include <amdis/common/Concepts.hpp>
 #include <amdis/common/ValueCategory.hpp>
-#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
+#include <amdis/linear_algebra/HierarchicWrapper.hpp>
 #include <amdis/utility/RangeType.hpp>
 #include <amdis/utility/TreeData.hpp>
 
 namespace AMDiS
 {
+  struct BoundaryType { int b; };
+
   /// Implements a boundary condition of Dirichlet-type.
   /**
    * By calling the methods \ref init() and \ref finish before and after
@@ -82,13 +84,13 @@ namespace AMDiS
       using Dune::Functions::interpolate;
       Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename RowBasis::LocalView::Tree>, Range>{},
         [&](auto id) {
-          auto rhsWrapper = wrapper(rhs.vector());
+          auto rhsWrapper = hierarchicVectorWrapper(rhs.vector());
           interpolate(id(rowBasis), rhsWrapper, values_, dirichletNodes_);
         });
 
       Dune::Hybrid::ifElse(std::is_same<RangeType_t<typename ColBasis::LocalView::Tree>, Range>{},
         [&](auto id) {
-          auto solutionWrapper = wrapper(solution.vector());
+          auto solutionWrapper = hierarchicVectorWrapper(solution.vector());
           interpolate(id(colBasis), solutionWrapper, values_, dirichletNodes_);
         });
 
diff --git a/src/amdis/DirichletBC.inc.hpp b/src/amdis/DirichletBC.inc.hpp
index 8f4fe3750c62984c61803f26431f50ebc0104730..c2536f44d9a2c2b181694a47395f5d1151e6aa88 100644
--- a/src/amdis/DirichletBC.inc.hpp
+++ b/src/amdis/DirichletBC.inc.hpp
@@ -3,9 +3,6 @@
 #include <dune/functions/functionspacebases/boundarydofs.hh>
 #include <dune/functions/functionspacebases/interpolate.hh>
 #include <dune/functions/functionspacebases/subspacebasis.hh>
-#include <amdis/LinearAlgebra.hpp>
-#include <amdis/linear_algebra/HierarchicWrapper.hpp>
-#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
 
 namespace AMDiS
 {
diff --git a/src/amdis/LinearAlgebra.hpp b/src/amdis/LinearAlgebra.hpp
index 569c7457b973b15f23645ab52a57761317b6c9d5..898404a001e91ab8cfd2df5f37f9f1f781f45b0f 100644
--- a/src/amdis/LinearAlgebra.hpp
+++ b/src/amdis/LinearAlgebra.hpp
@@ -2,30 +2,7 @@
 
 #include <amdis/linear_algebra/LinearSolverInterface.hpp>
 #include <amdis/linear_algebra/SolverInfo.hpp>
-
-#if defined(AMDIS_BACKEND_ISTL)
-
-#include <amdis/linear_algebra/istl/SystemVector.hpp>
-#include <amdis/linear_algebra/istl/SystemMatrix.hpp>
-#include <amdis/linear_algebra/istl/LinearSolver.hpp>
-
-#elif defined(AMDIS_BACKEND_MTL)
-
-#include <amdis/linear_algebra/mtl/SystemVector.hpp>
-#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
+#include <amdis/linear_algebra/mtl/DOFVector.hpp>
+#include <amdis/linear_algebra/mtl/DOFMatrix.hpp>
 #include <amdis/linear_algebra/mtl/LinearSolver.hpp>
 #include <amdis/linear_algebra/mtl/ITL_Solver.hpp>
-#include <amdis/linear_algebra/mtl/BITL_Solver.hpp>
-
-#elif defined(AMDIS_BACKEND_PETSC)
-
-#include <amdis/linear_algebra/petsc/SystemVector.hpp>
-#include <amdis/linear_algebra/petsc/SystemMatrix.hpp>
-#include <amdis/linear_algebra/petsc/LinearSolver.hpp>
-
-#else
-
-#error "Unknown linear algebra backend!. Set corresponding variable \
-        AMDIS_BACKEND_ISTL, AMDIS_BACKEND_MTL or AMDIS_BACKEND_PETSC."
-
-#endif
\ No newline at end of file
diff --git a/src/amdis/LocalAssemblerBase.hpp b/src/amdis/LocalAssemblerBase.hpp
index f6320729f7ddfc14d37557fa8635e03908948379..44e756a240510764a02121d99b93cbb63a1aa5e1 100644
--- a/src/amdis/LocalAssemblerBase.hpp
+++ b/src/amdis/LocalAssemblerBase.hpp
@@ -2,9 +2,9 @@
 
 #include <type_traits>
 
+#include <boost/numeric/mtl/mtl.hpp>
+
 #include <amdis/ContextGeometry.hpp>
-// #include <amdis/common/ConceptsBase.hpp>
-#include <amdis/common/TypeDefs.hpp>
 
 namespace AMDiS
 {
@@ -22,10 +22,13 @@ namespace AMDiS
     static_assert( numNodes == 1 || numNodes == 2,
       "VectorAssembler gets 1 Node, MatrixAssembler gets 2 Nodes!");
 
+    using ElementMatrix = mtl::mat::dense2D<double>; // TODO: choose correct value_type
+    using ElementVector = mtl::vec::dense_vector<double>;
+
     /// Either an ElementVector or an ElementMatrix (depending on the number of nodes)
     using ElementMatrixVector = std::conditional_t<
-      (sizeof...(Nodes)==1), Impl::ElementVector, std::conditional_t<
-      (sizeof...(Nodes)==2), Impl::ElementMatrix, void>>;
+      (sizeof...(Nodes)==1), ElementVector, std::conditional_t<
+      (sizeof...(Nodes)==2), ElementMatrix, void>>;
 
   public:
     /// Virtual destructor
diff --git a/src/amdis/ProblemStat.hpp b/src/amdis/ProblemStat.hpp
index e0e08a4360546c78f77b0a55e01a1c2bb119b533..0a53c64db5594422c3182544185357dde60ab29a 100644
--- a/src/amdis/ProblemStat.hpp
+++ b/src/amdis/ProblemStat.hpp
@@ -27,7 +27,6 @@
 #include <amdis/StandardProblemIteration.hpp>
 
 #include <amdis/common/TupleUtility.hpp>
-#include <amdis/common/TypeDefs.hpp>
 #include <amdis/common/Utility.hpp>
 
 #include <amdis/GridFunctions.hpp>
diff --git a/src/amdis/common/CMakeLists.txt b/src/amdis/common/CMakeLists.txt
index 03492346bd8b0e507998ca66b8be429a9b1f9a94..92e85985864553ed516303d2dad36a3edec12d01 100644
--- a/src/amdis/common/CMakeLists.txt
+++ b/src/amdis/common/CMakeLists.txt
@@ -19,7 +19,6 @@ install(FILES
     Size.hpp
     Tags.hpp
     TupleUtility.hpp
-    TypeDefs.hpp
     Utility.hpp
     ValueCategory.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/common)
diff --git a/src/amdis/common/TypeDefs.hpp b/src/amdis/common/TypeDefs.hpp
deleted file mode 100644
index 735d9cb6a50b7022c0ef0b414b8ed850a5f194ba..0000000000000000000000000000000000000000
--- a/src/amdis/common/TypeDefs.hpp
+++ /dev/null
@@ -1,21 +0,0 @@
-#pragma once
-
-#include <type_traits>
-
-#include <amdis/linear_algebra/Mtl.hpp>
-
-namespace AMDiS
-{
-  namespace Impl
-  {
-    using ElementVector = mtl::dense_vector<double>;
-    using ElementMatrix = mtl::dense2D<double>;
-
-    using SystemVector = mtl::dense_vector<double>;
-    using SystemMatrix = mtl::compressed2D<double>;
-
-  } // end namespace Impl
-
-  struct BoundaryType { int b; };
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/CMakeLists.txt b/src/amdis/linear_algebra/CMakeLists.txt
index 200338d51155f07071e20d69efbbafe1ab409853..ed57a3387c17abd2a34f6750aecf0e353d02806f 100644
--- a/src/amdis/linear_algebra/CMakeLists.txt
+++ b/src/amdis/linear_algebra/CMakeLists.txt
@@ -2,9 +2,7 @@
 
 install(FILES
     HierarchicWrapper.hpp
-    LinearAlgebraBase.hpp
     LinearSolverInterface.hpp
-    Mtl.hpp
     PreconditionerInterface.hpp
     RunnerInterface.hpp
     SolverInfo.hpp
diff --git a/src/amdis/linear_algebra/DOFVectorBase.hpp b/src/amdis/linear_algebra/DOFVectorInterface.hpp
similarity index 75%
rename from src/amdis/linear_algebra/DOFVectorBase.hpp
rename to src/amdis/linear_algebra/DOFVectorInterface.hpp
index 1b5586d1abf3d2396b5e274e362bac62ee4a2191..f2b5aa84301a305ff5682f01582793410d633132 100644
--- a/src/amdis/linear_algebra/DOFVectorBase.hpp
+++ b/src/amdis/linear_algebra/DOFVectorInterface.hpp
@@ -2,11 +2,11 @@
 
 namespace AMDiS
 {
-  class DOFVectorBase
+  class DOFVectorInterface
   {
   public:
     /// Virtual destructor
-    ~DOFVectorBase() = default;
+    ~DOFVectorInterface() = default;
 
     /// Change dimension of DOFVector to dimension of basis
     virtual void compress() = 0;
diff --git a/src/amdis/linear_algebra/HierarchicWrapper.hpp b/src/amdis/linear_algebra/HierarchicWrapper.hpp
index 94617e6de355bc2fbd9b7bb8ea3953fad3abb998..609cf4750388f4ed95b3df11e06ebcb08ffbef73 100644
--- a/src/amdis/linear_algebra/HierarchicWrapper.hpp
+++ b/src/amdis/linear_algebra/HierarchicWrapper.hpp
@@ -1,5 +1,6 @@
 #pragma once
 
+#include <boost/numeric/mtl/mtl_fwd.hpp>
 #include <dune/functions/functionspacebases/sizeinfo.hh>
 
 #include <amdis/utility/MultiIndex.hpp>
@@ -86,7 +87,7 @@ namespace AMDiS
     }
 
     template <class... Params>
-    std::size_t sizeImpl(mtl::dense_vector<Params...> const& v) const
+    std::size_t sizeImpl(mtl::vec::dense_vector<Params...> const& v) const
     {
       return mtl::size(v);
     }
diff --git a/src/amdis/linear_algebra/LinearAlgebraBase.hpp b/src/amdis/linear_algebra/LinearAlgebraBase.hpp
deleted file mode 100644
index 41c3d1601381e12200481f63d4ea7c23ac205a3c..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/LinearAlgebraBase.hpp
+++ /dev/null
@@ -1,65 +0,0 @@
-#pragma once
-
-namespace AMDiS
-{
-  namespace Impl
-  {
-    // Type-Trait needs to be specialized by linear algebra libraries
-    template <class Vector>
-    struct BaseVector { using type = Vector; };
-
-    template <class Matrix>
-    struct BaseMatrix { using type = Matrix; };
-  }
-
-  template <class Vector>
-  using BaseVector_t = typename Impl::BaseVector<Vector>::type;
-
-  template <class Matrix>
-  using BaseMatrix_t = typename Impl::BaseMatrix<Matrix>::type;
-
-  template <class T>
-  struct Triplet
-  {
-    std::size_t row;
-    std::size_t col;
-    T value;
-  };
-
-  namespace Impl
-  {
-    template <class Vector>
-    class BaseWrapper
-    {
-    public:
-      using value_type = typename Vector::value_type;
-      using size_type = typename Vector::size_type;
-
-      explicit BaseWrapper(Vector& vec) : vec(vec) {}
-
-      void resize(size_type s)
-      {
-        vec.resize(s);
-      }
-
-      size_type size() const
-      {
-        return vec.size();
-      }
-
-      value_type&       operator[](size_type i)       { return vec[i]; }
-      value_type const& operator[](size_type i) const { return vec[i]; }
-
-    protected:
-      Vector& vec;
-    };
-
-  } // end namespace Impl
-
-  template <class Vector>
-  Impl::BaseWrapper<std::remove_reference_t<Vector>> wrapper(Vector&& vec)
-  {
-    return {std::forward<Vector>(vec)};
-  }
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/LinearSolverInterface.hpp b/src/amdis/linear_algebra/LinearSolverInterface.hpp
index 5d1534f9f522db80170bdd3ee36bea85d25d56cd..3edf1d2dedee74c92922d5179f5a43b7d35d4290 100644
--- a/src/amdis/linear_algebra/LinearSolverInterface.hpp
+++ b/src/amdis/linear_algebra/LinearSolverInterface.hpp
@@ -12,7 +12,6 @@
  */
 
 #include <amdis/common/Utility.hpp>
-#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
 
 #include <amdis/linear_algebra/PreconditionerInterface.hpp>
 #include <amdis/linear_algebra/RunnerInterface.hpp>
@@ -25,12 +24,11 @@ namespace AMDiS
   class LinearSolverInterface
   {
   protected:
-    using RunnerBase = RunnerInterface<Matrix, BaseVector_t<VectorX>,
-                                               BaseVector_t<VectorB>>;
+    using RunnerBase = RunnerInterface<Matrix, VectorX, VectorB>;
 
   public:
     /// Destructor.
-    virtual ~LinearSolverInterface() {}
+    virtual ~LinearSolverInterface() = default;
 
     /// Public method to call in order to solve a linear system Ax = b.
     /**
diff --git a/src/amdis/linear_algebra/Mtl.hpp b/src/amdis/linear_algebra/Mtl.hpp
deleted file mode 100644
index 5d89ea427fd21f3cd57d2926100c19586478b532..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/Mtl.hpp
+++ /dev/null
@@ -1,32 +0,0 @@
-#pragma once
-
-#include <dune/common/ftraits.hh>
-
-#include <boost/numeric/mtl/matrix/compressed2D.hpp>
-#include <boost/numeric/mtl/matrix/dense2D.hpp>
-#include <boost/numeric/mtl/vector/dense_vector.hpp>
-
-namespace Dune
-{
-  template <class Elt, class Parameters>
-  struct FieldTraits<mtl::compressed2D<Elt,Parameters>>
-  {
-    using field_type = typename FieldTraits<Elt>::field_type;
-    using real_type = typename FieldTraits<Elt>::real_type;
-  };
-
-  template <class Value, class Parameters>
-  struct FieldTraits<mtl::dense2D<Value,Parameters>>
-  {
-    using field_type = typename FieldTraits<Value>::field_type;
-    using real_type = typename FieldTraits<Value>::real_type;
-  };
-
-  template <class Value, class Parameters>
-  struct FieldTraits<mtl::dense_vector<Value,Parameters>>
-  {
-    using field_type = typename FieldTraits<Value>::field_type;
-    using real_type = typename FieldTraits<Value>::real_type;
-  };
-
-} // end namespace Dune
diff --git a/src/amdis/linear_algebra/PreconditionerInterface.hpp b/src/amdis/linear_algebra/PreconditionerInterface.hpp
index 692a35250d904b2b5e63bf71f261ffdceaee705f..055ee01b886b6ee3e600583b811fc00c869c6265 100644
--- a/src/amdis/linear_algebra/PreconditionerInterface.hpp
+++ b/src/amdis/linear_algebra/PreconditionerInterface.hpp
@@ -7,7 +7,7 @@ namespace AMDiS
   struct PreconditionerInterface
   {
     /// Virtual destructor.
-    virtual ~PreconditionerInterface() {}
+    virtual ~PreconditionerInterface() = default;
 
     /// Is called a the beginning of a solution procedure
     virtual void init(Matrix const& fullMatrix) = 0;
diff --git a/src/amdis/linear_algebra/RunnerInterface.hpp b/src/amdis/linear_algebra/RunnerInterface.hpp
index 0c81227ce650ab7b3a18929dc7c0e9480bc5aa89..082c023452f2b7a58ebf4a69018045b67ab65fda 100644
--- a/src/amdis/linear_algebra/RunnerInterface.hpp
+++ b/src/amdis/linear_algebra/RunnerInterface.hpp
@@ -14,7 +14,7 @@ namespace AMDiS
 
   public:
     /// virtual destructor
-    virtual ~RunnerInterface() {}
+    virtual ~RunnerInterface() = default;
 
     /// Is called at the beginning of a solution procedure
     virtual void init(Matrix const& A) = 0;
diff --git a/src/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp b/src/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp
deleted file mode 100644
index 8a5739689d56b9406fc696096879fb6a6ff1dc84..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/BITL_Preconditioner.hpp
+++ /dev/null
@@ -1,38 +0,0 @@
-#pragma once
-
-#include "itl/block_diagonal.hpp"
-
-#include <amdis/CreatorMap.hpp>
-#include <amdis/linear_algebra/mtl/ITL_Preconditioner.hpp>
-
-namespace AMDiS
-{
-
-  template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
-  class DefaultCreators<PreconditionerInterface<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>>
-  {
-    using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
-    using PreconBase = PreconditionerInterface<Matrix, Vector>;
-
-    template <template<class> class ITLPrecon>
-    using PreconCreator
-      = typename Preconditioner<Matrix, Vector, ITLPrecon<Matrix>>::Creator;
-
-    using Map = CreatorMap<PreconBase>;
-
-  public:
-    static void init()
-    {
-      auto pc_diag = new PreconCreator<DiagonalPreconditioner>;
-      Map::addCreator("diag", pc_diag);
-      Map::addCreator("jacobi", pc_diag);
-
-      auto pc_id = new PreconCreator<IdentityPreconditioner>;
-      Map::addCreator("identity", pc_id);
-      Map::addCreator("no", pc_id);
-
-      Map::addCreator("default", pc_id);
-    }
-  };
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/BITL_Solver.hpp b/src/amdis/linear_algebra/mtl/BITL_Solver.hpp
deleted file mode 100644
index c011bdc769633a0d49c109d654763dc56efd6b41..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/BITL_Solver.hpp
+++ /dev/null
@@ -1,83 +0,0 @@
-#pragma once
-
-#include <amdis/CreatorMap.hpp>
-#include <amdis/linear_algebra/mtl/ITL_Solver.hpp>
-
-namespace AMDiS
-{
-  /// Adds default creators for linear solvers based on `BlockMTLMatrix`.
-  /**
-   * Adds creators for block matrix aware solvers.
-   * - *cg*: conjugate gradient method, \see cg_solver_type
-   * - *cgs*: stabilized conjugate gradient mtheod, \see cgs_solver_type
-   * - *bcgs*: stabilized bi-conjugate gradient method, \see bicgstab_type
-   * - *tfqmr*: Transposed-Free Quasi-Minimal Residual method, \see tfqmr_solver_type
-   * - *bcgsl*: stabilized BiCG(ell) method, \see bicgstab_ell_type
-   * - *gmres*: Generalized minimal residula method, \see gmres_type
-   * - *fgmres*: Flexible GMRES, \see fgmres_type
-   * - *minres*: Minimal residul method, \see minres_solver_type
-   * - *gcr*: Generalized conjugate residual method, \see gcr_type
-   * - *umfpack*: external UMFPACK solver, \see UmfpackRunner
-   **/
-  template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
-  class DefaultCreators<LinearSolverInterface<BlockMTLMatrix<SubMatrix, _N, _M>, Vector, Vector>>
-  {
-    using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
-    using SolverBase = LinearSolverInterface<Matrix, Vector, Vector>;
-
-    template <class ITLSolver>
-    using SolverCreator
-      = typename LinearSolver<Matrix, Vector,
-          KrylovRunner<ITLSolver, Matrix, BaseVector_t<Vector>>>::Creator;
-
-#ifdef HAVE_UMFPACK
-    using UmfpackSolverCreator
-      = typename LinearSolver<Matrix, Vector,
-          UmfpackRunner<Matrix, BaseVector_t<Vector>>>::Creator;
-#endif
-
-    using Map = CreatorMap<SolverBase>;
-
-  public:
-    static void init()
-    {
-      auto cg = new SolverCreator<cg_solver_type>;
-      Map::addCreator("cg", cg);
-
-      auto cgs = new SolverCreator<cgs_solver_type>;
-      Map::addCreator("cgs", cgs);
-
-      auto bcgs = new SolverCreator<bicgstab_type>;
-      Map::addCreator("bicgstab", bcgs);
-      Map::addCreator("bcgs", bcgs);
-
-      auto tfqmr = new SolverCreator<tfqmr_solver_type>;
-      Map::addCreator("tfqmr", tfqmr);
-
-      auto bcgsl = new SolverCreator<bicgstab_ell_type>;
-      Map::addCreator("bicgstab_ell", bcgsl);
-      Map::addCreator("bcgsl", bcgsl);
-
-      auto gmres = new SolverCreator<gmres_type>;
-      Map::addCreator("gmres", gmres);
-
-      auto fgmres = new SolverCreator<fgmres_type>;
-      Map::addCreator("fgmres", fgmres);
-
-      auto minres = new SolverCreator<minres_solver_type>;
-      Map::addCreator("minres", minres);
-
-      auto gcr = new SolverCreator<gcr_type>;
-      Map::addCreator("gcr", gcr);
-
-#ifdef HAVE_UMFPACK
-      auto umfpack = new UmfpackSolverCreator;
-      Map::addCreator("umfpack", umfpack);
-      Map::addCreator("direct", umfpack);
-#endif
-
-      Map::addCreator("default", gmres);
-    }
-  };
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp b/src/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp
deleted file mode 100644
index ef661b363706530a4b5672b1a9b2f808f91295ca..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/BlockMTLMatrix.hpp
+++ /dev/null
@@ -1,208 +0,0 @@
-/** \file BlockMTLMatrix.hpp */
-
-#pragma once
-
-#include <array>
-
-#include <boost/numeric/mtl/matrices.hpp>
-
-#include <amdis/common/Utility.hpp>
-#include <amdis/common/Literals.hpp>
-#include <amdis/common/Loops.hpp>
-#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
-
-namespace AMDiS
-{
-  /// A wrapper for AMDiS::SolverMatrix to be used in MTL/ITL solvers
-  template <class MTLMatrix, std::size_t _N, std::size_t _M>
-  class BlockMTLMatrix
-      : public std::array<std::array<MTLMatrix, _M>, _N>
-  {
-    using Self = BlockMTLMatrix;
-
-  public:
-    /// The index/size - type
-    using size_type  = typename MTLMatrix::size_type;
-
-    /// The type of the elements of the MTLMatrix
-    using value_type = typename MTLMatrix::value_type;
-
-    /// The underlying mtl matrix type
-    using BaseMatrix = MTLMatrix;
-
-  public:
-    /// Return the (R,C)'th matrix block
-    template <std::size_t R, std::size_t C>
-    auto& operator()(const index_t<R>, const index_t<C>)
-    {
-      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
-      return std::get<C>(std::get<R>(*this));
-    }
-
-    /// Return the (R,C)'th matrix block
-    template <std::size_t R, std::size_t C>
-    auto const& operator()(const index_t<R>, const index_t<C>) const
-    {
-      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
-      return std::get<C>(std::get<R>(*this));
-    }
-
-    /// Return the number of row blocks
-    static constexpr std::size_t N() { return _N; }
-
-    /// Return the number of column blocks
-    static constexpr std::size_t M() { return _M; }
-
-    /// perform blockwise multiplication A*b -> x
-    template <class VectorIn, class VectorOut, class Assign>
-    void mult(VectorIn const& b, VectorOut& x, Assign) const
-    {
-      // create iranges to access array blocks
-      std::array<mtl::irange, _N> r_rows;
-      std::array<mtl::irange, _M> r_cols;
-      getRanges(r_rows, r_cols);
-
-      forEach(range_<0, _N>, [&](const auto _i)
-      {
-        bool first = true;
-
-        // a reference to the i'th block of x
-        VectorOut x_i(x[r_rows[_i]]);
-        forEach(range_<0, _M>, [&](const auto _j)
-        {
-          auto const& A_ij = this->operator()(_i, _j);
-          if (num_rows(A_ij) > 0 && A_ij.nnz() > 0) {
-            // a reference to the j'th block of b
-            const VectorIn b_j(b[r_cols[_j]]);
-
-            if (first) {
-              Assign::first_update(x_i, A_ij * b_j);
-              first = false;
-            }
-            else {
-              Assign::update(x_i, A_ij * b_j);
-            }
-          }
-        });
-      });
-    }
-
-    /// A Multiplication operator returns a multiplication-expresssion.
-    /// Calls \ref mult internally.
-    template <class VectorIn>
-    mtl::vec::mat_cvec_multiplier<Self, VectorIn>
-    operator*(VectorIn const& v) const
-    {
-      return {*this, v};
-    }
-
-    /// Fill an array of irange corresponding to the row-sizes, used
-    /// to access sub-vectors
-    void getRowRanges(std::array<mtl::irange, _N>& r_rows) const
-    {
-      std::size_t start = 0;
-      forEach(range_<0, _N>, [&](const auto _r) {
-        std::size_t finish = start + num_rows((*this)(_r, 0_c));
-        r_rows[_r].set(start, finish);
-        start = finish;
-      });
-    }
-
-    /// Fill an array of irange corresponding to the column-sizes, used
-    /// to access sub-vectors
-    void getColRanges(std::array<mtl::irange, _M>& r_cols) const
-    {
-      std::size_t start = 0;
-      forEach(range_<0, _M>, [&](const auto _c) {
-        std::size_t finish = start + num_cols((*this)(0_c, _c));
-        r_cols[_c].set(start, finish);
-        start = finish;
-      });
-    }
-
-    /// Fill two arrays of irange corresponding to row and column sizes.
-    /// \see getRowRanges() and \see getColRanges()
-    void getRanges(std::array<mtl::irange, _N>& r_rows,
-                   std::array<mtl::irange, _M>& r_cols) const
-    {
-      getRowRanges(r_rows);
-      getColRanges(r_cols);
-    }
-  };
-
-
-  namespace Impl
-  {
-    /// Specialization of Impl::MTLMatrix from \file LinearAlgebraBase.hpp
-    template <class MTLMatrix, std::size_t _N, std::size_t _M>
-    struct BaseMatrix<BlockMTLMatrix<MTLMatrix, _N, _M>>
-    {
-      using type = MTLMatrix;
-    };
-  }
-
-  /// Return the number of overall rows of a BlockMTLMatrix
-  template <class MTLMatrix, std::size_t _N, std::size_t _M>
-  inline std::size_t num_rows(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
-  {
-    std::size_t nRows = 0;
-    forEach(range_<0, _N>, [&](const auto _r) {
-      nRows += num_rows(A(_r, 0_c));
-    });
-    return nRows;
-  }
-
-  /// Return the number of overall columns of a BlockMTLMatrix
-  template <class MTLMatrix, std::size_t _N, std::size_t _M>
-  inline std::size_t num_cols(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
-  {
-    std::size_t nCols = 0;
-    forEach(range_<0, _M>, [&](const auto _c) {
-      nCols += num_cols(A(0_c, _c));
-    });
-    return nCols;
-  }
-
-  /// Return the size, i.e. rows*columns of a BlockMTLMatrix
-  template <class MTLMatrix, std::size_t _N, std::size_t _M>
-  inline std::size_t size(BlockMTLMatrix<MTLMatrix, _N, _M> const& A)
-  {
-    return num_rows(A) * num_cols(A);
-  }
-
-  /// Nullify a BlockMTLMatrix, i.e. nullify each block.
-  template <class MTLMatrix, std::size_t _N, std::size_t _M>
-  inline void set_to_zero(BlockMTLMatrix<MTLMatrix, _N, _M>& A)
-  {
-    forEach(range_<0, _N>, [&](const auto _r) {
-      forEach(range_<0, _M>, [&](const auto _c) {
-        set_to_zero(A(_r,_c));
-      });
-    });
-  }
-
-} // end namespace AMDiS
-
-
-/// \cond HIDDEN_SYMBOLS
-namespace mtl
-{
-  template <class MTLMatrix, std::size_t _N, std::size_t _M>
-  struct Collection<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>>
-  {
-    using value_type = typename MTLMatrix::value_type;
-    using size_type  = typename MTLMatrix::size_type;
-  };
-
-  namespace ashape
-  {
-    template <class MTLMatrix, std::size_t _N, std::size_t _M>
-    struct ashape_aux<AMDiS::BlockMTLMatrix<MTLMatrix, _N, _M>>
-    {
-      using type = nonscal;
-    };
-
-  } // end namespace ashape
-
-} // end namespace mtl
-/// \endcond
diff --git a/src/amdis/linear_algebra/mtl/BlockMTLVector.hpp b/src/amdis/linear_algebra/mtl/BlockMTLVector.hpp
deleted file mode 100644
index 8bcc2466e34fd78c0db172f78d24a262b4bcdb89..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/BlockMTLVector.hpp
+++ /dev/null
@@ -1,186 +0,0 @@
-#pragma once
-
-#include <array>
-#include <type_traits>
-
-#include <boost/numeric/mtl/utility/irange.hpp>
-
-#include <amdis/common/Mpl.hpp>
-#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
-#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
-
-namespace AMDiS
-{
-  /// A simple block-vector class
-  template <class MTLVector, std::size_t N>
-  class BlockMTLVector
-      : public std::array<MTLVector, N>
-  {
-  public:
-    /// The type of the blocks
-    using BaseVector = MTLVector;
-
-    /// The index/size - type
-    using size_type  = typename MTLVector::size_type;
-
-    /// The type of the elements of the MTLVector
-    using value_type = typename MTLVector::value_type;
-  };
-
-
-  namespace Impl
-  {
-    /// Specialization of Impl::BaseVector from \file LinearAlgebraBase.hpp
-    template <class MTLVector, std::size_t N>
-    struct BaseVector<BlockMTLVector<MTLVector, N>>
-    {
-      using type = MTLVector;
-    };
-  }
-
-
-  /// Return the number of overall rows of a BlockMTLVector
-  template <class MTLVector, std::size_t N>
-  inline std::size_t num_rows(BlockMTLVector<MTLVector, N> const& vec)
-  {
-    std::size_t nRows = 0;
-    forEach(range_<0, N>, [&](const auto _i) {
-      nRows += num_rows(std::get<_i>(vec));
-    });
-    return nRows;
-  }
-
-  /// Return the number of overall columns of a BlockMTLVector
-  template <class MTLVector, std::size_t N>
-  inline std::size_t num_cols(BlockMTLVector<MTLVector, N> const& vec)
-  {
-    return 1;
-  }
-
-  /// Return the size, i.e. rows*columns of a BlockMTLVector
-  template <class MTLVector, std::size_t N>
-  inline std::size_t size(BlockMTLVector<MTLVector, N> const& vec)
-  {
-    return num_rows(vec);
-  }
-
-  /// Nullify a BlockMTLVector, i.e. nullify each block.
-  template <class MTLVector, std::size_t N>
-  inline void set_to_zero(BlockMTLVector<MTLVector, N>& vec)
-  {
-    forEach(range_<0, N>, [&](const auto _i) {
-      set_to_zero(std::get<_i>(vec));
-    });
-  }
-
-
-  /// A wrapper, that creates a contiguos vector corresponding to a block-vector
-  /// and copy the value on construction and eventually back on destruction, if
-  /// required.
-  template <class BlockVector, class Vector, class NonConstBlockVector>
-  class BlockVectorWrapper;
-
-  // specialization for BlockMTLVector
-  template <class BlockVector, class Vector, class MTLVector, std::size_t N>
-  class BlockVectorWrapper<BlockVector, Vector, BlockMTLVector<MTLVector, N>>
-  {
-    static_assert( std::is_same< std::remove_const_t<BlockVector>, BlockMTLVector<MTLVector, N> >::value,
-                   "This specialization is for BlockMTLVectors only.");
-
-  public:
-    explicit BlockVectorWrapper(BlockVector& blockVector,
-                                bool copyBack = !std::is_const<BlockVector>::value)
-      : blockVector(blockVector)
-      , vector_(num_rows(blockVector))
-      , copyBack(copyBack)
-    {
-      assignTo();
-    }
-
-    ~BlockVectorWrapper()
-    {
-      if (copyBack)
-        assignFrom(bool_<!std::is_const<BlockVector>::value>);
-    }
-
-    /// Return a reference to the block-vector
-    BlockVector const&  getBlockVector() const  { return blockVector; }
-
-    /// Return a reference to the contiguose-vector
-    Vector&       vector()       { return vector_; }
-    Vector const& vector() const { return vector_; }
-
-  private:
-    /// copy from block-vector to vector
-    void assignTo()
-    {
-      std::size_t start = 0;
-      for (std::size_t r = 0; r < N; ++r) {
-        std::size_t finish = start + num_rows(blockVector[r]);
-        mtl::irange range(start, finish);
-
-        vector_[range] = blockVector[r];
-        start = finish;
-      }
-    }
-
-    /// do not copy vector to block-vector since this is const
-    template <bool Copy>
-    std::enable_if_t<!Copy> assignFrom(bool_t<Copy>) { /* Do nothing */ }
-
-    /// copy from vector to block-vector
-    template <bool Copy>
-    std::enable_if_t<Copy> assignFrom(bool_t<Copy>)
-    {
-      std::size_t start = 0;
-      for (std::size_t r = 0; r < N; ++r) {
-        std::size_t finish = start + num_rows(blockVector[r]);
-        mtl::irange range(start, finish);
-
-        blockVector[r] = vector_[range];
-        start = finish;
-      }
-    }
-
-  private:
-    BlockVector& blockVector;
-    Vector vector_;
-
-    /// Copy data back to block-vector on destruction
-    bool copyBack;
-  };
-
-
-  // Specialization for non-block-vectors
-  template <class BlockVector, class Vector>
-  class BlockVectorWrapper<BlockVector, Vector, Vector>
-  {
-    static_assert( std::is_same< std::remove_const_t<BlockVector>, Vector >::value,
-                   "This specialization is for contiguose vectors only.");
-  public:
-    explicit BlockVectorWrapper(BlockVector& vector, bool = false)
-      : vector_(vector)
-    {}
-
-    /// Return a reference to the vector
-    BlockVector const& getBlockVector() const { return vector_; }
-    BlockVector&       vector()            { return vector_; }
-    BlockVector const& vector()      const { return vector_; }
-
-  private:
-    BlockVector& vector_;
-  };
-
-  template <class BlockVector, class Vector>
-  using BlockVectorWrapper_t
-    = BlockVectorWrapper<BlockVector, Vector, std::remove_const_t<BlockVector>>;
-
-
-  template <class BlockVector, class Vector = MTLDenseVector<typename BlockVector::value_type>>
-  auto blockWrapper(BlockVector& bvec)
-  {
-    return BlockVectorWrapper_t<BlockVector, Vector>{bvec};
-  }
-
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/CMakeLists.txt b/src/amdis/linear_algebra/mtl/CMakeLists.txt
index 0e6d6a6706c12d6bc4d862d2dfb09b1ff4649065..e5d7c4d7641ac63672f3b2068cebe024e5901cb6 100644
--- a/src/amdis/linear_algebra/mtl/CMakeLists.txt
+++ b/src/amdis/linear_algebra/mtl/CMakeLists.txt
@@ -1,25 +1,11 @@
-#dune_library_add_sources(amdis SOURCES
-#    SystemMatrix.cpp
-#    SystemVector.cpp
-#)
-
 install(FILES
-    BITL_Preconditioner.hpp
-    BITL_Solver.hpp
-    BlockMTLMatrix.hpp
-    BlockMTLVector.hpp
-    Copy.hpp
     DOFMatrix.hpp
     DOFVector.hpp
     ITL_Preconditioner.hpp
     ITL_Solver.hpp
     KrylovRunner.hpp
     LinearSolver.hpp
-    Mapper.hpp
-    MTLDenseVector.hpp
     Preconditioner.hpp
-    SystemMatrix.hpp
-    SystemVector.hpp
     UmfpackRunner.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linear_algebra/mtl)
 
diff --git a/src/amdis/linear_algebra/mtl/Copy.hpp b/src/amdis/linear_algebra/mtl/Copy.hpp
deleted file mode 100644
index 1724b311ac95756a312d1c5bf8d5f3b348e8e842..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/Copy.hpp
+++ /dev/null
@@ -1,52 +0,0 @@
-#pragma once
-
-#include <algorithm>
-
-#include <amdis/linear_algebra/mtl/BlockMTLVector.hpp>
-#include <amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
-
-#include <amdis/linear_algebra/mtl/SystemVector.hpp>
-#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
-
-#include <amdis/linear_algebra/mtl/Mapper.hpp>
-
-namespace AMDiS
-{
-  // requires MatrixType to have an inserter
-  template <class Value, std::size_t _N, std::size_t _M, class TargetType>
-  void copy(BlockMTLMatrix<Value, _N, _M> const& source, TargetType& target)
-  {
-    using Mapper = std::conditional_t<(_N == _M), BlockMapper, RectangularMapper>;
-    using value_type = typename TargetType::value_type;
-    using BaseInserter = mtl::mat::inserter<TargetType, mtl::operations::update_plus<value_type>>;
-    using MappedInserter = typename Mapper::template Inserter<BaseInserter>::type;
-
-    std::size_t nnz = 0;
-    for (std::size_t rb = 0; rb < _N; ++rb)
-      for (std::size_t cb = 0; cb < _M; ++cb)
-        nnz += source[rb][cb].nnz();
-
-    {
-      target.change_dim(num_rows(source), num_cols(source));
-      set_to_zero(target);
-
-      Mapper mapper(source);
-      MappedInserter ins(target, mapper, int(1.2 * nnz / target.dim1()));
-
-      for (std::size_t rb = 0; rb < _N; ++rb) {
-        mapper.setRow(rb);
-        for (std::size_t cb = 0; cb < _M; ++cb) {
-	        mapper.setCol(cb);
-          ins << source[rb][cb];
-        }
-      }
-    }
-  }
-
-  template <class FeSpaces, class Value, class TargetType>
-  void copy(SystemMatrix<FeSpaces, Value> const& source, TargetType& target)
-  {
-    copy(source.getMatrix(), target);
-  }
-
-} // end namespace AMDiS
\ No newline at end of file
diff --git a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
index 0b748f5da353c82c0eb8ad87a1e6dfdc6f375564..7adbf1b879864e0688f4a48912797915a3e7c76e 100644
--- a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
@@ -1,7 +1,6 @@
 #pragma once
 
 #include <list>
-#include <map>
 #include <memory>
 #include <string>
 #include <vector>
@@ -11,15 +10,19 @@
 #include <boost/numeric/mtl/utility/property_map.hpp>
 #include <boost/numeric/mtl/utility/range_wrapper.hpp>
 
-#include <dune/common/std/optional.hh>
-#include <dune/common/reservedvector.hh>
-#include <dune/functions/functionspacebases/flatmultiindex.hh>
-
 #include <amdis/Output.hpp>
 #include <amdis/common/ClonablePtr.hpp>
+#include <amdis/utility/MultiIndex.hpp>
 
 namespace AMDiS
 {
+  template <class T>
+  struct Triplet
+  {
+    std::size_t row, cols;
+    T value;
+  };
+
   /// \brief The basic container that stores a base matrix and a corresponding
   /// row/column feSpace.
   template <class RowBasisType,
@@ -43,9 +46,6 @@ namespace AMDiS
     /// The type of the elements of the DOFMatrix
     using value_type = ValueType;
 
-    /// The underlying field-type (typically the same as \ref value_type)
-    using field_type = typename BaseMatrix::value_type;
-
     /// The type of the inserter used when filling the matrix. NOTE: need not be public
     using Inserter = mtl::mat::inserter<BaseMatrix, mtl::operations::update_plus<value_type>>;
 
@@ -61,10 +61,10 @@ namespace AMDiS
     /// Constructor. Takes pointer of data-matrix.
     DOFMatrix(RowBasis const& rowBasis,
               ColBasis const& colBasis,
-              BaseMatrix& matrix_ref)
+              BaseMatrix& matrix)
       : rowBasis_(&rowBasis)
       , colBasis_(&colBasis)
-      , matrix_(matrix_ref)
+      , matrix_(matrix)
     {}
 
     /// Return the row-basis \ref rowBasis_ of the matrix
@@ -82,58 +82,48 @@ namespace AMDiS
     /// Return a reference to the data-matrix \ref matrix
     BaseMatrix& matrix()
     {
-      assert( !insertionMode_ );
+      assert( !inserter_ );
       return *matrix_;
     }
 
     /// Return a reference to the data-matrix \ref matrix
     BaseMatrix const& matrix() const
     {
-      assert( !insertionMode_ );
+      assert( !inserter_ );
       return *matrix_;
     }
 
     /// Return the size of the row \ref feSpace
-    size_type N() const
+    size_type rows() const
     {
       return rowBasis_->size();
     }
 
     /// Return the size of the column \ref feSpace
-    size_type M() const
+    size_type cols() const
     {
       return colBasis_->size();
     }
 
 
-    using FlatIndex = Dune::Functions::FlatMultiIndex<size_type>;
-    using BlockIndex = Dune::ReservedVector<size_type,1>;
-
     /// \brief Returns an update-proxy of the inserter, to inserte/update a value at
     /// position (\p r, \p c) in the matrix. Need an insertionMode inserter, that can
     /// be created by \ref init and must be closed by \ref finish after insertion.
-    auto operator()(FlatIndex row, FlatIndex col)
+    template <class Index>
+    auto operator()(Index row, Index col)
     {
-      size_type r = row[0], c = col[0];
-      test_exit_dbg( insertionMode_, "Inserter not initilized!");
-      test_exit_dbg( r < N() && c < M() ,
-          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
+      size_type r = flatMultiIndex(row), c = flatMultiIndex(col);
+      test_exit_dbg(inserter_, "Inserter not initilized!");
+      test_exit_dbg(r < rows() && c < cols(),
+          "Indices out of range [0,", rows(), ")x[0,", cols(), ")" );
       return (*inserter_)[r][c];
     }
 
-    auto operator()(BlockIndex row, BlockIndex col)
-    {
-      size_type r = row[0], c = col[0];
-      test_exit_dbg( insertionMode_, "Inserter not initilized!");
-      test_exit_dbg( r < N() && c < M() ,
-          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
-      return (*inserter_)[r][c];
-    }
 
     /// Create inserter. Assumes that no inserter is currently active on this matrix.
     void init(bool prepareForInsertion)
     {
-      test_exit(!insertionMode_ && !inserter_, "Matrix already in insertion mode!");
+      test_exit(!inserter_, "Matrix already in insertion mode!");
 
       calculateSlotSize();
       matrix_->change_dim(rowBasis_->dimension(), colBasis_->dimension());
@@ -142,7 +132,6 @@ namespace AMDiS
       if (prepareForInsertion) {
         set_to_zero(*matrix_);
         inserter_ = new Inserter(*matrix_, slotSize_);
-        insertionMode_ = true;
       }
     }
 
@@ -151,8 +140,7 @@ namespace AMDiS
     void finish()
     {
       delete inserter_;
-      inserter_ = NULL;
-      insertionMode_ = false;
+      inserter_ = nullptr;
     }
 
     /// Return the number of nonzeros in the matrix
@@ -287,9 +275,6 @@ namespace AMDiS
     /// A pointer to the inserter. Created in \ref init and destroyed in \ref finish
     Inserter* inserter_ = nullptr;
 
-    /// A flag that indicates whether the matrix is in insertion mode
-    bool insertionMode_ = false;
-
     /// The size of the slots used to insert values per row
     int slotSize_ = MIN_NNZ_PER_ROW;
   };
diff --git a/src/amdis/linear_algebra/mtl/DOFVector.hpp b/src/amdis/linear_algebra/mtl/DOFVector.hpp
index 1679e1ed7ab3d52d762a383ee8e425bf7156aa24..f5ff128ce2b19f3f5163b89d22837854cf5fdf65 100644
--- a/src/amdis/linear_algebra/mtl/DOFVector.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFVector.hpp
@@ -4,15 +4,12 @@
 #include <memory>
 #include <string>
 
-#include <dune/common/reservedvector.hh>
-#include <dune/functions/functionspacebases/flatmultiindex.hh>
-#include <dune/functions/functionspacebases/interpolate.hh>
+#include <boost/numeric/mtl/vector/dense_vector.hpp>
 
 #include <amdis/Output.hpp>
 #include <amdis/common/ClonablePtr.hpp>
 #include <amdis/common/ScalarTypes.hpp>
-#include <amdis/linear_algebra/DOFVectorBase.hpp>
-#include <amdis/linear_algebra/mtl/MTLDenseVector.hpp>
+#include <amdis/linear_algebra/DOFVectorInterface.hpp>
 #include <amdis/utility/MultiIndex.hpp>
 
 namespace AMDiS
@@ -20,7 +17,7 @@ namespace AMDiS
   /// The basic container that stores a base vector and a corresponding basis
   template <class GlobalBasis, class ValueType = double>
   class DOFVector
-      : public DOFVectorBase
+      : public DOFVectorInterface
   {
     using Self = DOFVector;
 
@@ -29,7 +26,7 @@ namespace AMDiS
     using Basis = GlobalBasis;
 
     /// The type of the base vector
-    using BaseVector = MTLDenseVector<ValueType>;
+    using BaseVector = mtl::vec::dense_vector<ValueType>;
 
     /// The index/size - type
     using size_type  = typename GlobalBasis::size_type;
@@ -46,24 +43,24 @@ namespace AMDiS
     }
 
     /// Constructor. Takes reference to existing BaseVector
-    DOFVector(GlobalBasis const& basis, BaseVector& vector_ref)
+    DOFVector(GlobalBasis const& basis, BaseVector& vector)
       : basis_(&basis)
-      , vector_(vector_ref)
+      , vector_(vector)
     {}
 
-    /// Return the basis \ref basis of the vector
+    /// Return the basis \ref basis_ of the vector
     Basis const& basis() const
     {
         return *basis_;
     }
 
-    /// Return the data-vector \ref vector
+    /// Return the data-vector \ref vector_
     BaseVector const& vector() const
     {
       return *vector_;
     }
 
-    /// Return the data-vector \ref vector
+    /// Return the data-vector \ref vector_
     BaseVector& vector()
     {
       return *vector_;
@@ -95,8 +92,8 @@ namespace AMDiS
     template <class Index>
     value_type const& operator[](Index idx) const
     {
-      auto const i = flatMultiIndex(idx);
-      test_exit_dbg( i < num_rows(*vector_) ,
+      size_type i = flatMultiIndex(idx);
+      test_exit_dbg(i < num_rows(*vector_),
         "Index ", i, " out of range [0,", num_rows(*vector_), ")" );
       return (*vector_)[i];
     }
@@ -105,20 +102,20 @@ namespace AMDiS
     template <class Index>
     value_type& operator[](Index idx)
     {
-      auto const i = flatMultiIndex(idx);
-      test_exit_dbg( i < num_rows(*vector_) ,
+      size_type i = flatMultiIndex(idx);
+      test_exit_dbg(i < num_rows(*vector_),
         "Index ", i, " out of range [0,", num_rows(*vector_), ")" );
       return (*vector_)[i];
     }
 
-    /// Scale each DOFVector by the factor \p s.
-    template <class Scalar>
-    std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
-    operator*=(Scalar s)
-    {
-      (*vector_) *= s;
-      return *this;
-    }
+    // /// Scale each DOFVector by the factor \p s.
+    // template <class Scalar>
+    // std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
+    // operator*=(Scalar s)
+    // {
+    //   (*vector_) *= s;
+    //   return *this;
+    // }
 
     /// Sets each DOFVector to the scalar \p s.
     template <class Scalar>
@@ -150,7 +147,7 @@ namespace AMDiS
     -> DOFVector<GlobalBasis, double>;
 
   template <class GlobalBasis, class ValueType>
-  DOFVector(GlobalBasis const& basis, MTLDenseVector<ValueType>& coefficients)
+  DOFVector(GlobalBasis const& basis, mtl::vec::dense_vector<ValueType>& coefficients)
     -> DOFVector<GlobalBasis, ValueType>;
 #endif
 
@@ -165,7 +162,7 @@ namespace AMDiS
   /// Constructor a dofvector from given basis, name, and coefficient vector
   template <class GlobalBasis, class ValueType>
   DOFVector<GlobalBasis, ValueType>
-  makeDOFVector(GlobalBasis const& basis, MTLDenseVector<ValueType>& coefficients)
+  makeDOFVector(GlobalBasis const& basis, mtl::vec::dense_vector<ValueType>& coefficients)
   {
     return {basis, coefficients};
   }
diff --git a/src/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp b/src/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp
index 2e8ac6ef19f5ed4dbd640d5628e66fb90fdcb9d7..42e25d04fb451bbd41363d9d318ed7b2aed2b4c9 100644
--- a/src/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp
+++ b/src/amdis/linear_algebra/mtl/ITL_Preconditioner.hpp
@@ -1,5 +1,3 @@
-/** \file ITL_Preconditioner.h */
-
 #pragma once
 
 // MTL4 headers
diff --git a/src/amdis/linear_algebra/mtl/ITL_Solver.hpp b/src/amdis/linear_algebra/mtl/ITL_Solver.hpp
index c60b372cd5f30d18eb71b83b38cb37adc187c174..27a76f5c460dd2cf34af7f4db01359e52b9b8d5b 100644
--- a/src/amdis/linear_algebra/mtl/ITL_Solver.hpp
+++ b/src/amdis/linear_algebra/mtl/ITL_Solver.hpp
@@ -1,21 +1,9 @@
-/** \file ITL_Solver.h */
-
 #pragma once
 
 #include <string>
 
 // MTL4 includes
 #include <boost/numeric/itl/itl.hpp>
-// #include <boost/numeric/itl/krylov/bicg.hpp>
-// #include <boost/numeric/itl/krylov/bicgstab_2.hpp>
-// #include <boost/numeric/itl/krylov/bicgstab_ell.hpp>
-// #include <boost/numeric/itl/krylov/bicgstab.hpp>
-// #include <boost/numeric/itl/krylov/cg.hpp>
-// #include <boost/numeric/itl/krylov/cgs.hpp>
-// #include <boost/numeric/itl/krylov/gmres.hpp>
-// #include <boost/numeric/itl/krylov/idr_s.hpp>
-// #include <boost/numeric/itl/krylov/qmr.hpp>
-// #include <boost/numeric/itl/krylov/tfqmr.hpp>
 
 // more solvers defined in AMDiS
 #include <amdis/linear_algebra/mtl/itl/minres.hpp>
@@ -406,12 +394,12 @@ namespace AMDiS
     template <class ITLSolver>
     using SolverCreator
       = typename LinearSolver<Matrix, Vector,
-          KrylovRunner<ITLSolver, Matrix, BaseVector_t<Vector>>>::Creator;
+          KrylovRunner<ITLSolver, Matrix, Vector>>::Creator;
 
 #ifdef HAVE_UMFPACK
     using UmfpackSolverCreator
       = typename LinearSolver<Matrix, Vector,
-          UmfpackRunner<Matrix, BaseVector_t<Vector>>>::Creator;
+          UmfpackRunner<Matrix, Vector>>::Creator;
 #endif
 
     using Map = CreatorMap<SolverBase>;
diff --git a/src/amdis/linear_algebra/mtl/KrylovRunner.hpp b/src/amdis/linear_algebra/mtl/KrylovRunner.hpp
index 45b9ea2a1475e777a9cf7a671b098cf4e34443ad..a226241f895e2eba3597577c6e53a1843884936f 100644
--- a/src/amdis/linear_algebra/mtl/KrylovRunner.hpp
+++ b/src/amdis/linear_algebra/mtl/KrylovRunner.hpp
@@ -1,5 +1,3 @@
-/** \file ITL_Runner.h */
-
 #pragma once
 
 #include <string>
@@ -14,7 +12,6 @@
 #include <amdis/linear_algebra/SolverInfo.hpp>
 
 #include <amdis/linear_algebra/mtl/ITL_Preconditioner.hpp>
-#include <amdis/linear_algebra/mtl/BITL_Preconditioner.hpp>
 
 namespace AMDiS
 {
@@ -22,10 +19,11 @@ namespace AMDiS
    *
    * \brief
    * Wrapper class for different MTL4 itl-solvers. These solvers
-   * are parametrized by Matrix- and Vector.
+   * are parametrized by Matrix and Vector.
    **/
   template <class ITLSolver, class Matrix, class Vector>
-  class KrylovRunner : public RunnerInterface<Matrix, Vector>
+  class KrylovRunner
+      : public RunnerInterface<Matrix, Vector>
   {
     using Super = RunnerInterface<Matrix, Vector>;
     using PreconBase = typename Super::PreconBase;
@@ -33,71 +31,71 @@ namespace AMDiS
   public:
     /// Constructor.
     explicit KrylovRunner(std::string prefix)
-      : itlSolver(prefix)
+      : itlSolver_(prefix)
     {
       initPrecon(prefix);
 
-      Parameters::get(prefix + "->max iteration", maxIter);
-      Parameters::get(prefix + "->print cycle", printCycle);
+      Parameters::get(prefix + "->max iteration", maxIter_);
+      Parameters::get(prefix + "->print cycle", printCycle_);
     }
 
-    /// Implements \ref RunnerInterface::getLeftPrecon(), Returns \ref L
+    /// Implements \ref RunnerInterface::getLeftPrecon(), Returns \ref L_
     virtual std::shared_ptr<PreconBase> getLeftPrecon() override
     {
-      return L;
+      return L_;
     }
 
-    /// Implements \ref RunnerInterface::getRightPrecon(), Returns \ref R
+    /// Implements \ref RunnerInterface::getRightPrecon(), Returns \ref R_
     virtual std::shared_ptr<PreconBase> getRightPrecon() override
     {
-      return R;
+      return R_;
     }
 
-    /// Set a new left preconditioner \ref L
+    /// Set a new left preconditioner \ref L_
     void setLeftPrecon(std::shared_ptr<PreconBase> const& precon)
     {
-      L = precon;
+      L_ = precon;
     }
 
-    /// Set a new right preconditioner \ref R
+    /// Set a new right preconditioner \ref R_
     void setRightPrecon(std::shared_ptr<PreconBase> const& precon)
     {
-      R = precon;
+      R_ = precon;
     }
 
     /// Implementation of \ref RunnerInterface::init()
     virtual void init(Matrix const& A) override
     {
-      L->init(A);
-      R->init(A);
+      L_->init(A);
+      R_->init(A);
     }
 
     /// Implementation of \ref RunnerInterface::exit()
     virtual void exit() override
     {
-      L->exit();
-      R->exit();
+      L_->exit();
+      R_->exit();
     }
 
     /// Implementation of \ref RunnerInterface::solve()
     virtual int solve(Matrix const& A, Vector& x, Vector const& b,
                       SolverInfo& solverInfo) override
     {
-      test_exit(bool(L), "There is no left preconditioner");
-      test_exit(bool(R), "There is no right preconditioner");
+      test_exit(bool(L_), "There is no left preconditioner");
+      test_exit(bool(R_), "There is no right preconditioner");
 
       auto r = Vector(b);
       if (two_norm(x) != 0)
         r -= A * x;
 
       // print information about the solution process
-      itl::cyclic_iteration<double> iter(r, maxIter, solverInfo.getRelTolerance(),
-                                                     solverInfo.getAbsTolerance(),
-                                                     printCycle);
+      itl::cyclic_iteration<double> iter(r, maxIter_, solverInfo.getRelTolerance(),
+                                                      solverInfo.getAbsTolerance(),
+                                                      printCycle_);
       iter.set_quite(solverInfo.getInfo() == 0);
       iter.suppress_resume(solverInfo.getInfo() == 0);
 
-      int error = itlSolver(A, x, b, *L, *R, iter);
+      int error = itlSolver_(A, x, b, *L_, *R_, iter);
       solverInfo.setAbsResidual(iter.resid());
       solverInfo.setRelResidual(iter.relresid());
 
@@ -118,23 +116,22 @@ namespace AMDiS
       auto rightCreator
         = CreatorMap<PreconBase>::getCreator(rightPrecon, prefix + "->right precon");
 
-      L = leftCreator->create();
-      R = rightCreator->create();
+      L_ = leftCreator->create();
+      R_ = rightCreator->create();
     }
 
   private:
     /// The itl solver object that performes the actual solve
-    ITLSolver itlSolver;
+    ITLSolver itlSolver_;
 
     /// Pointer to the left and right preconditioner
-    std::shared_ptr<PreconBase> L, R;
+    std::shared_ptr<PreconBase> L_, R_;
 
     /// The maximal number of iterations
-    std::size_t maxIter = 1000;
+    std::size_t maxIter_ = 1000;
 
     /// Print solver residuals every \ref printCycle 's iteration
-    std::size_t printCycle = 100;
-
+    std::size_t printCycle_ = 100;
   };
 
 } // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/LinearSolver.hpp b/src/amdis/linear_algebra/mtl/LinearSolver.hpp
index 856f44dd664320ab08b3dfcf454d7f5bb0dede01..ed44e4c2ea06d0d2d3d20715e28cd85ae1074a2c 100644
--- a/src/amdis/linear_algebra/mtl/LinearSolver.hpp
+++ b/src/amdis/linear_algebra/mtl/LinearSolver.hpp
@@ -8,7 +8,6 @@
 #include <amdis/CreatorInterface.hpp>
 #include <amdis/Output.hpp>
 #include <amdis/linear_algebra/LinearSolverInterface.hpp>
-#include <amdis/linear_algebra/mtl/BlockMTLVector.hpp>
 
 namespace AMDiS
 {
@@ -42,13 +41,13 @@ namespace AMDiS
   public:
     /// Constructor
     explicit LinearSolver(std::string prefix)
-      : runner(std::make_shared<Runner>(prefix))
+      : runner_(std::make_shared<Runner>(prefix))
     {}
 
     /// Implements \ref LinearSolverInterface::getRunner()
     virtual std::shared_ptr<RunnerBase> getRunner() override
     {
-      return runner;
+      return runner_;
     }
 
   private:
@@ -59,27 +58,23 @@ namespace AMDiS
       Dune::Timer t;
       if (solverInfo.doCreateMatrixData()) {
         // init matrix or wrap block-matrix or ...
-        runner->init(A);
+        runner_->init(A);
       }
 
-      // wrappers copy the data back on destruction
-      auto X = blockWrapper(x);
-      auto B = blockWrapper(b);
-
       if (solverInfo.getInfo() > 0) {
         msg("fill MTL4 matrix needed ", t.elapsed(), " seconds");
       }
 
-      int error = runner->solve(A, X.vector(), B.vector(), solverInfo);
+      int error = runner_->solve(A, x, b, solverInfo);
       solverInfo.setError(error);
 
       if (!solverInfo.doStoreMatrixData())
-        runner->exit();
+        runner_->exit();
     }
 
   private:
     /// redirect the implementation to a runner. Implements a \ref RunnerInterface
-    std::shared_ptr<Runner> runner;
+    std::shared_ptr<Runner> runner_;
   };
 
 } // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/MTLDenseVector.hpp b/src/amdis/linear_algebra/mtl/MTLDenseVector.hpp
deleted file mode 100644
index 4e08e4698dca62edaafda5c74647266decfc4911..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/MTLDenseVector.hpp
+++ /dev/null
@@ -1,50 +0,0 @@
-#pragma once
-
-
-#include <boost/numeric/mtl/vector/dense_vector.hpp>
-
-#include <amdis/linear_algebra/LinearAlgebraBase.hpp>
-
-namespace AMDiS
-{
-  template <class Value>
-  using MTLDenseVector = mtl::vec::dense_vector<Value>;
-
-  namespace Impl
-  {
-    template <class Vector>
-    class MTLWrapper
-      : public BaseWrapper<Vector>
-    {
-      using Super = BaseWrapper<Vector>;
-
-    public:
-      using value_type = typename Super::value_type;
-      using size_type = typename Super::size_type;
-
-      template <class Vector_>
-      MTLWrapper(Vector_&& vec)
-        : Super(std::forward<Vector_>(vec))
-      {}
-
-      void resize(size_type s)
-      {
-        this->vec.change_dim(s);
-      }
-
-      size_type size() const
-      {
-        return num_rows(this->vec);
-      }
-    };
-
-  } // end namespace Impl
-
-  /// wrap the mtl dense_vector class to provide resize and size member-functions
-  template <class Value>
-  Impl::MTLWrapper<MTLDenseVector<Value>> wrapper(MTLDenseVector<Value>& vec)
-  {
-    return {vec};
-  }
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/Mapper.hpp b/src/amdis/linear_algebra/mtl/Mapper.hpp
deleted file mode 100644
index 6249f0949554a1e13a49edcac61b74855e868df0..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/Mapper.hpp
+++ /dev/null
@@ -1,241 +0,0 @@
-#pragma once
-
-#include <vector>
-
-#include <boost/numeric/mtl/matrix/mapped_inserter.hpp>
-
-#include <amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
-#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
-
-namespace AMDiS
-{
-
-  /**
-   * \brief BaseClass for all mapper types
-   *
-   * A mapper assigned a local row/column index of a cell of
-   * a block matrix to a global matrix index.
-   * Call \ref setRow and \ref setCol first, to select a block
-   * in the block matrix. Then get with \ref row, respective \ref col,
-   * the global matrix index to the assigned local index.
-   **/
-  template <class Derived>
-  struct MapperBase
-  {
-    using size_type = std::size_t;
-
-    template <class BaseInserter>
-    struct Inserter {
-      using type = mtl::mat::mapped_inserter<BaseInserter, Derived>;
-    };
-
-    /// set the current block row
-    void setRow(size_type rb) { self().setRow(rb); }
-
-    /// set the current block columns
-    void setCol(size_type cb) { self().setCol(cb); }
-
-    /// return global matrix row, for local row in the current matrix block
-    size_type row(size_type r) const { return self().row(r); }
-
-    /// return global matrix column, for local column in the current matrix block
-    size_type col(size_type c) const { return self().col(c); }
-
-    /// return overall number of rows
-    size_type getNumRows() const { return self().getNumRows(); }
-
-    /// return overall number of columns
-    size_type getNumCols() const { return self().getNumCols(); }
-
-    size_type getNumRows(size_type comp) const { return self().getNumRows(comp); }
-    size_type getNumCols(size_type comp) const { return self().getNumCols(comp); }
-
-    /// return number of components/blocks
-    size_type getNumComponents() const { return self().getNumComponents(); }
-
-    size_type row(size_type rb, size_type cb, size_type r)
-    {
-      setRow(rb);
-      setCol(cb);
-      return row(r);
-    }
-
-    size_type col(size_type rb, size_type cb, size_type c)
-    {
-      setRow(rb);
-      setCol(cb);
-      return col(c);
-    }
-
-    Derived&       self()       { return static_cast<Derived&>(*this); }
-    Derived const& self() const { return static_cast<Derived const&>(*this); }
-  };
-
-
-  /// Mapper implementation for non-parallel block matrices
-  struct BlockMapper : public MapperBase<BlockMapper>
-  {
-    /// Default constructor
-    BlockMapper() = default;
-    BlockMapper(BlockMapper const& other) = default;
-    BlockMapper& operator=(BlockMapper const& other) = default;
-
-    /// Constructor for BlockMTLMatrix
-    template <class Value, std::size_t _N>
-    explicit BlockMapper(BlockMTLMatrix<Value, _N, _N> const& mat)
-      : nComp(_N)
-      , rowOffset(0), colOffset(0)
-      , nrow(num_rows(mat)), ncol(nrow)
-      , sizes(nComp)
-    {
-      for (std::size_t i = 0; i < _N; ++i)
-        sizes[i] = num_rows(mat[i][i]);
-    }
-
-    /// Constructor for SystemMatrix
-    template <class FeSpaces, class Value>
-    explicit BlockMapper(SystemMatrix<FeSpaces, Value> const& mat)
-      : BlockMapper(mat.getMatrix())
-    {}
-
-    /// Constructor for single DOFMatrix
-    template <class FeSpace, class Value>
-    explicit BlockMapper(DOFMatrix<FeSpace, Value> const& mat)
-      : nComp(1)
-      , rowOffset(0), colOffset(0)
-      , nrow(0), ncol(0)
-      , sizes(nComp)
-    {
-      sizes[0] = mat.N();
-      nrow += sizes[0];
-      ncol = nrow;
-    }
-
-    /// Constructor for system with equal components
-    BlockMapper(size_type nComp, size_type nDOFperComp)
-    : nComp(nComp),
-      rowOffset(0), colOffset(0),
-      nrow(nComp*nDOFperComp), ncol(nrow),
-      sizes(nComp)
-    {
-       for (size_type i = 0; i < nComp; ++i)
-         sizes[i] = nDOFperComp;
-    }
-
-    /// calculate row offset for row component \param r
-    void setRow(size_type rb)
-    {
-      test_exit_dbg(rb <= sizes.size(), "row nr out of range!");
-      rowOffset = sum(rb);
-    }
-
-    /// calculate column offset for col component \param c
-    void setCol(size_type cb)
-    {
-      test_exit_dbg(cb <= sizes.size(), "column nr out of range!");
-      colOffset = sum(cb);
-    }
-
-    size_type row(size_type r) const { return r + rowOffset; }
-    size_type col(size_type c) const { return c + colOffset; }
-
-    size_type getNumRows() const { return nrow; }
-    size_type getNumCols() const { return ncol; }
-
-    size_type getNumRows(size_type comp) const { return sizes[comp]; }
-    size_type getNumCols(size_type comp) const { return sizes[comp]; }
-
-    size_type getNumComponents() const { return nComp; }
-
-  private: // methods
-
-    ///compute the sum of sizes from [0, end)
-    size_type sum(size_type end) const
-    {
-      size_type ret = 0;
-      for (size_type i = 0; i < end; ++i)
-	ret += sizes[i];
-      return ret;
-    }
-
-  private: // variables
-
-    size_type nComp = 0;
-    size_type rowOffset = 0;
-    size_type colOffset = 0;
-    size_type nrow = 0;
-    size_type ncol = 0;
-
-    std::vector<size_type> sizes;
-  };
-
-
-  /// Mapper implementation for non-parallel rectangular block matrices
-  struct RectangularMapper : public MapperBase<RectangularMapper>
-  {
-    /// Constructor for block-matrices
-    template <class Value, std::size_t _N, std::size_t _M>
-    explicit RectangularMapper(BlockMTLMatrix<Value, _N, _M> const& mat)
-      : nRowComp(_N), nColComp(_M)
-      , rowOffset(0), colOffset(0)
-      , nrow(num_rows(mat)), ncol(num_cols(mat))
-      , sizes_rows(nRowComp), sizes_cols(nColComp)
-    {
-      for (size_type i = 0; i < _N; ++i)
-        sizes_rows[i] = num_rows(mat[i][0]);
-      for (size_type j = 0; j < _M; ++j)
-        sizes_cols[j] = num_cols(mat[0][j]);
-    }
-
-    /// calculate row offset for row component \param r
-    void setRow(size_type rb)
-    {
-      test_exit_dbg(rb <= sizes_rows.size(), "row nr out of range!");
-      rowOffset = sum(rb, sizes_rows);
-    }
-
-    /// calculate column offset for col component \param c
-    void setCol(size_type cb)
-    {
-      test_exit_dbg(cb <= sizes_cols.size(), "column nr out of range!");
-      colOffset = sum(cb, sizes_cols);
-    }
-
-    size_type row(size_type r) const { return r + rowOffset; }
-    size_type col(size_type c) const { return c + colOffset; }
-
-    size_type getNumRows() const { return nrow; }
-    size_type getNumCols() const { return ncol; }
-
-    size_type getNumRows(size_type comp) const { return sizes_rows[comp]; }
-    size_type getNumCols(size_type comp) const { return sizes_cols[comp]; }
-
-    size_type getNumComponents() const { return nRowComp; }
-    size_type getNumRowComponents() const { return nRowComp; }
-    size_type getNumColComponents() const { return nColComp; }
-
-  private: // methods
-
-    ///compute the sum of sizes from [0, end)
-    size_type sum(size_type end, std::vector<size_type> const& sizes) const
-    {
-      size_type ret = 0;
-      for (size_type i = 0; i < end; ++i)
-	ret += sizes[i];
-      return ret;
-    }
-
-  private: // variables
-
-    size_type nRowComp = 0;
-    size_type nColComp = 0;
-    size_type rowOffset = 0;
-    size_type colOffset = 0;
-    size_type nrow = 0;
-    size_type ncol = 0;
-
-    std::vector<size_type> sizes_rows;
-    std::vector<size_type> sizes_cols;
-  };
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/Preconditioner.hpp b/src/amdis/linear_algebra/mtl/Preconditioner.hpp
index 822d1221841e73676e74b282f74571fa1460fc3b..6d9922ea9d3f945b270c603a7a5caca08fa9ba04 100644
--- a/src/amdis/linear_algebra/mtl/Preconditioner.hpp
+++ b/src/amdis/linear_algebra/mtl/Preconditioner.hpp
@@ -1,5 +1,3 @@
-/** \file ITL_Preconditioner.h */
-
 #pragma once
 
 #include <amdis/linear_algebra/PreconditionerInterface.hpp>
@@ -31,31 +29,31 @@ namespace AMDiS
     /// Implementation of \ref PreconditionerBase::init()
     virtual void init(Matrix const& fullMatrix) override
     {
-      precon.reset(new PreconRunner(fullMatrix));
+      precon_.reset(new PreconRunner(fullMatrix));
     }
 
     /// Implementation of \ref PreconditionerInterface::exit()
     virtual void exit() override
     {
-      precon.reset();
+      precon_.reset();
     }
 
     /// Implementation of \ref PreconditionerBase::solve()
     virtual void solve(Vector const& vin, Vector& vout) const override
     {
-      test_exit_dbg(bool(precon), "No preconditioner initialized!");
-      precon->solve(vin, vout);
+      test_exit_dbg(bool(precon_), "No preconditioner initialized!");
+      precon_->solve(vin, vout);
     }
 
     /// Implementation of \ref PreconditionerBase::adjoint_solve()
     virtual void adjoint_solve(Vector const& vin, Vector& vout) const override
     {
-      test_exit_dbg(bool(precon), "No preconditioner initialized!");
-      precon->adjoint_solve(vin, vout);
+      test_exit_dbg(bool(precon_), "No preconditioner initialized!");
+      precon_->adjoint_solve(vin, vout);
     }
 
   private:
-    std::shared_ptr<PreconRunner> precon;
+    std::shared_ptr<PreconRunner> precon_;
   };
 
 } // namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/SystemMatrix.cpp b/src/amdis/linear_algebra/mtl/SystemMatrix.cpp
deleted file mode 100644
index 6182e50796fac5ccbea79435d9e3c444df2f91f3..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/SystemMatrix.cpp
+++ /dev/null
@@ -1,10 +0,0 @@
-#define AMDIS_NO_EXTERN_SYSTEMMATRIX
-#include <amdis/linear_algebra/mtl/SystemMatrix.hpp>
-#undef AMDIS_NO_EXTERN_SYSTEMMATRIX
-
-
-namespace AMDiS
-{
-  // explicit template instatiation
-  // template class SystemMatrix<typename TestTraits<2,1,2>::FeSpaces>;
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/SystemMatrix.hpp b/src/amdis/linear_algebra/mtl/SystemMatrix.hpp
deleted file mode 100644
index a674fd16a44fc8afa275cfd2a00faa4ebaeea3de..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/SystemMatrix.hpp
+++ /dev/null
@@ -1,207 +0,0 @@
-#pragma once
-
-#include <vector>
-#include <string>
-#include <memory>
-#include <tuple>
-
-#include <boost/numeric/mtl/matrix/dense2D.hpp>
-#include <boost/numeric/mtl/matrix/compressed2D.hpp>
-
-#include <amdis/Output.hpp>
-#include <amdis/common/TupleUtility.hpp>
-#include <amdis/common/Utility.hpp>
-#include <amdis/linear_algebra/mtl/BlockMTLMatrix.hpp>
-#include <amdis/linear_algebra/mtl/DOFMatrix.hpp>
-
-// for explicit instantiation
-#include <amdis/ProblemStatTraits.hpp>
-
-namespace AMDiS
-{
-  /// \cond HIDDEN_SYMBOLS
-  namespace Impl
-  {
-    // DOFMatrices = std::tuple< std::tuple< DOFMatrix<RowFeSpace, ColFeSpace, ValueType>... >... >
-    template <class DOFMatrices>
-    struct BuildDOFMatrices
-    {
-      template <std::size_t R>
-      using DOFMatrixRow = std::tuple_element_t<R, DOFMatrices>;
-
-      template <std::size_t R, std::size_t C>
-      using DOFMatrixRowCol = std::tuple_element_t<C, DOFMatrixRow<R>>;
-
-      // add arg to repeated constructor argument list
-      template <std::size_t... R, class RowFeSpaces, class ColFeSpaces, class MultiMatrix>
-      static DOFMatrices make(Indices<R...>,
-                              RowFeSpaces&& rowFeSpaces,
-                              ColFeSpaces&& colFeSpace,
-                              MultiMatrix&& multiMatrix)
-      {
-        return DOFMatrices{make_row(
-          index_<R>,
-          MakeSeq_t<std::tuple_size<std::decay_t<ColFeSpaces>>::value>(),
-          std::get<R>(std::forward<RowFeSpaces>(rowFeSpaces)),
-          std::forward<ColFeSpaces>(colFeSpace),
-          std::forward<MultiMatrix>(multiMatrix)
-          )...};
-      }
-
-      template <std::size_t R, std::size_t... C, class RowFeSpace, class ColFeSpaces, class MultiMatrix>
-      static DOFMatrixRow<R> make_row(index_t<R>,
-                                      Indices<C...>,
-                                      RowFeSpace&& rowFeSpace,
-                                      ColFeSpaces&& colFeSpace,
-                                      MultiMatrix&& multiMatrix)
-      {
-        return DOFMatrixRow<R>{DOFMatrixRowCol<R,C>(
-          std::forward<RowFeSpace>(rowFeSpace),
-          std::get<C>(std::forward<ColFeSpaces>(colFeSpace)),
-          "matrix[" + std::to_string(R) + "][" + std::to_string(C) + "]",
-          multiMatrix(index_<R>, index_<C>)
-          )...};
-      }
-    };
-
-    // construction method to construct a tuple of tuples of DOFMatrices
-    template <class DOFMatrices, class RowFeSpaces, class ColFeSpaces, class MultiMatrix>
-    DOFMatrices buildDOFMatrices(RowFeSpaces&& rowFeSpaces,
-                                  ColFeSpaces&& colFeSpace,
-                                  MultiMatrix&& multiMatrix)
-    {
-        return BuildDOFMatrices<DOFMatrices>::make(
-            MakeSeq_t<std::tuple_size<std::decay_t<RowFeSpaces>>::value>(), // traverse rows first
-            std::forward<RowFeSpaces>(rowFeSpaces),
-            std::forward<ColFeSpaces>(colFeSpace),
-            std::forward<MultiMatrix>(multiMatrix));
-    }
-
-  } // end namespace Impl
-  /// \endcond
-
-
-  /// \brief Container that repesents multiple data-matrices.
-  /**
-    *  Represents a fixed-size matrix of \ref mtl::compressed2D + a tuple of corresponding
-    *  feSpaces. This (R,C)'th matrix combined with the (R,C)'th FeSpace
-    *  builds a \ref DOFMatrix and can be return by the \ref operator().
-    **/
-  template <class FeSpaces, class ValueType = double>
-  class SystemMatrix
-  {
-    template <class RowFeSpace>
-    struct MultiMatrixRowGenerator
-    {
-      template <class ColFeSpace>
-      using DOFMatrixGenerator = DOFMatrix<RowFeSpace, ColFeSpace, ValueType>;
-      using DOFMatrices  = MakeTuple_t<DOFMatrixGenerator, FeSpaces>;
-    };
-
-    template <class RowFeSpace>
-    using DOFMatrixGenerator = typename MultiMatrixRowGenerator<RowFeSpace>::DOFMatrices;
-
-  public:
-    /// The number of blocks per row/column
-    static constexpr std::size_t nComponent = std::tuple_size<FeSpaces>::value;
-    AMDIS_STATIC_ASSERT( nComponent > 0 );
-
-    /// The type of the matrix in each block
-    using BaseMatrix = mtl::compressed2D<ValueType>;
-
-    /// The type of the matrix of DOFMatrices
-    using DOFMatrices = MakeTuple_t<DOFMatrixGenerator, FeSpaces>;
-
-    /// The type of the matrix of base-matrices
-    using MultiMatrix = BlockMTLMatrix<BaseMatrix, nComponent, nComponent>;
-
-
-    /// Constructor that constructs new base-matrices
-    explicit SystemMatrix(FeSpaces const& feSpaces)
-      : feSpaces(feSpaces)
-      , matrix()
-      , dofmatrices(Impl::buildDOFMatrices<DOFMatrices>(feSpaces, feSpaces, matrix))
-    {}
-
-
-    /// Return the number of blocks per row
-    static constexpr std::size_t N()
-    {
-      return std::tuple_size<FeSpaces>::value;
-    }
-
-    /// Return the number of blocks per column
-    static constexpr std::size_t M()
-    {
-      return std::tuple_size<FeSpaces>::value;
-    }
-
-    /// Return the (R,C)'th underlaying base matrix
-    template <std::size_t R, std::size_t C>
-    auto& getMatrix(const index_t<R> _r = {}, const index_t<C> _c = {})
-    {
-      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
-      return matrix(_r, _c);
-    }
-
-    /// Return the (R,C)'th underlaying base matrix
-    template <std::size_t R, std::size_t C>
-    auto const& getMatrix(const index_t<R> _r = {}, const index_t<C> _c = {}) const
-    {
-      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
-      return matrix(_r, _c);
-    }
-
-    /// Return the underlying multi matrix
-    MultiMatrix&        getMatrix()       { return matrix; }
-    MultiMatrix const&  getMatrix() const { return matrix; }
-
-    /// Return the I'th row finite element space
-    template <std::size_t I = 0>
-    auto& getRowFeSpace(const index_t<I> = {}) const
-    {
-      static_assert(I < N(), "Index out of range [0,N)");
-      return std::get<I>(feSpaces);
-    }
-
-    /// Return the I'th column finite element space
-    template <std::size_t I = 0>
-    auto const& getColFeSpace(const index_t<I> = {}) const
-    {
-      static_assert(I < M(), "Index out of range [0,M)");
-      return std::get<I>(feSpaces);
-    }
-
-    /// Return the (R,C)'th DOFMatrix
-    template <std::size_t R = 0, std::size_t C = 0>
-    auto& operator()(const index_t<R> = {}, const index_t<C> = {})
-    {
-      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
-      return std::get<C>(std::get<R>(dofmatrices));
-    }
-
-    /// Return the (R,C)'th DOFMatrix
-    template <std::size_t R = 0, std::size_t C = 0>
-    auto const& operator()(const index_t<R> = {}, const index_t<C> = {}) const
-    {
-      static_assert(R < N() && C < M(), "Indices out of range [0,N)x[0,M)");
-      return std::get<C>(std::get<R>(dofmatrices));
-    }
-
-  private:
-    /// a tuple of feSpaces for rows and columns of the matrix
-    FeSpaces const& feSpaces;
-
-    /// a tuple of data-matrices
-    MultiMatrix matrix;
-
-    /// a tuple of DOFMatrices
-    DOFMatrices dofmatrices;
-  };
-
-// #ifndef AMDIS_NO_EXTERN_SYSTEMMATRIX
-//     // explicit instantiation in SystemMatrix.cpp
-//     extern template class SystemMatrix<typename TestTraits<2,1,2>::FeSpaces>;
-// #endif
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/SystemVector.cpp b/src/amdis/linear_algebra/mtl/SystemVector.cpp
deleted file mode 100644
index 6fe23e982458095ec98ca59a2564cafb541cc638..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/SystemVector.cpp
+++ /dev/null
@@ -1,10 +0,0 @@
-#define AMDIS_NO_EXTERN_SYSTEMVECTOR
-#include <amdis/linear_algebra/mtl/SystemVector.hpp>
-#undef AMDIS_NO_EXTERN_SYSTEMVECTOR
-
-
-namespace AMDiS
-{
-  // explicit template instatiation
-  // template class SystemVector<typename TestTraits<2,1,2>::FeSpaces>;
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/SystemVector.hpp b/src/amdis/linear_algebra/mtl/SystemVector.hpp
deleted file mode 100644
index dcf7e32815db9f859e58a61d675b24b04dd62bab..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/mtl/SystemVector.hpp
+++ /dev/null
@@ -1,263 +0,0 @@
-#pragma once
-
-#include <algorithm>
-#include <memory>
-#include <string>
-#include <tuple>
-#include <vector>
-
-#include <amdis/Output.hpp>
-#include <amdis/common/Loops.hpp>
-#include <amdis/common/ScalarTypes.hpp>
-#include <amdis/common/TupleUtility.hpp>
-#include <amdis/common/Utility.hpp>
-
-#include <amdis/linear_algebra/mtl/BlockMTLVector.hpp>
-#include <amdis/linear_algebra/mtl/DOFVector.hpp>
-
-// for explicit instantiation
-#include <amdis/ProblemStatTraits.hpp>
-
-namespace AMDiS
-{
-  namespace Impl
-  {
-    // DOFVectors = std::tuple< DOFVector<FeSpace, ValueType>... >
-    template <class DOFVectors>
-    struct BuildDOFVectors
-    {
-      template <std::size_t I>
-      using DOFVectorIdx = std::tuple_element_t<I, DOFVectors>;
-
-      template <std::size_t... I, class FeSpaces, class MultiVector>
-      static DOFVectors make(Indices<I...>,
-                              FeSpaces&& feSpaces,
-                              std::vector<std::string> const& names,
-                              MultiVector&& multiVector)
-      {
-        return DOFVectors{DOFVectorIdx<I>(
-          std::get<I>(std::forward<FeSpaces>(feSpaces)),
-          names[I],
-          multiVector[I]
-          )...};
-      }
-    };
-
-    // construction method to construct a tuple of DOFVectors
-    template <class DOFVectors, class FeSpaces, class MultiVector>
-    DOFVectors buildDOFVectors(FeSpaces&& feSpaces,
-                                    std::vector<std::string> const& names,
-                                    MultiVector&& multiVector)
-    {
-      return BuildDOFVectors<DOFVectors>::make(
-          MakeSeq_t<std::tuple_size<std::decay_t<FeSpaces>>::value>(),
-          std::forward<FeSpaces>(feSpaces),
-          names,
-          std::forward<MultiVector>(multiVector));
-    }
-
-  } // end namespace Impl
-
-
-  /// \brief Container that repesents multiple data-Vectors of different value types.
-  /**
-    *  Represents a std::array of \ref mtl::dense_vector + a tuple of corresponding
-    *  feSpaces. This I'th matrix combined with the I'th FeSpace
-    *  builds a \ref DOFVector  and can be return by the \ref operator().
-    **/
-  template <class FeSpaces, class ValueType = double>
-  class SystemVector
-  {
-    using Self = SystemVector;
-
-    template <class FeSpace>
-    using DOFVectorGenerator = DOFVector<FeSpace, ValueType>;
-
-  public:
-    /// The number of blocks
-    static constexpr std::size_t nComponents = std::tuple_size<FeSpaces>::value;
-    AMDIS_STATIC_ASSERT( nComponents > 0 );
-
-    /// The type of the block-vector
-    using BaseVector = typename DOFVectorGenerator<std::tuple_element_t<0, FeSpaces>>::BaseVector;
-    using MultiVector = BlockMTLVector<BaseVector, nComponents>;
-
-    /// The type of the vector of DOFVectors
-    using DOFVectors  = MakeTuple_t<DOFVectorGenerator, FeSpaces>;
-
-    /// Constructor that constructs new base-vectors
-    SystemVector(FeSpaces const& feSpaces, std::vector<std::string> const& names)
-      : feSpaces(feSpaces)
-      , names(names)
-      , vector_{}
-      , dofvectors(Impl::buildDOFVectors<DOFVectors>(feSpaces, names, vector_))
-    {
-      compress();
-    }
-
-
-    /// Return the number of blocks
-    static constexpr std::size_t size()
-    {
-      return std::tuple_size<FeSpaces>::value;
-    }
-
-    /// Return the number of blocks
-    int count() const
-    {
-      return int(size());
-    }
-
-    /// Return a shared pointer to the I'th underlaying base vector
-    template <std::size_t I>
-    auto& getVector(const index_t<I> = {})
-    {
-      static_assert(I < size(), "Index out of range [0,SIZE)" );
-      return vector_[I];
-    }
-
-    /// Return a shared pointer to the I'th underlaying base vector
-    template <std::size_t I>
-    auto const& getVector(const index_t<I> = {}) const
-    {
-      static_assert(I < size(), "Index out of range [0,SIZE)" );
-      return vector_[I];
-    }
-
-    /// Return the underlying multi vector
-    MultiVector&        vector()       { return vector_; }
-    MultiVector const&  vector() const { return vector_; }
-
-    /// Return the I'th finite element space
-    template <std::size_t I = 0>
-    auto const& getFeSpace(const index_t<I> = {}) const
-    {
-      static_assert(I < size(), "Index out of range [0,SIZE)" );
-      return std::get<I>(feSpaces);
-    }
-
-    /// Return the name of the i'th DOFVector
-    std::string getName(std::size_t i) const
-    {
-        return names[i];
-    }
-
-    /// Return the I'th DOFVector
-    template <std::size_t I>
-    auto& operator[](const index_t<I>)
-    {
-      static_assert(I < size(), "Index out of range [0,SIZE)" );
-      return std::get<I>(dofvectors);
-    }
-
-    /// Return the I'th DOFVector
-    template <std::size_t I>
-    auto const& operator[](const index_t<I>) const
-    {
-      static_assert(I < size(), "Index out of range [0,SIZE)" );
-      return std::get<I>(dofvectors);
-    }
-
-    /// Return the I'th DOFVector
-    template <std::size_t I>
-    auto& getDOFVector(const index_t<I> = {})
-    {
-      static_assert(I < size(), "Index out of range [0,SIZE)" );
-      return std::get<I>(dofvectors);
-    }
-
-    /// Return the I'th DOFVector
-    template <std::size_t I>
-    auto const& getDOFVector(const index_t<I> = {}) const
-    {
-      static_assert(I < size(), "Index out of range [0,SIZE)" );
-      return std::get<I>(dofvectors);
-    }
-
-    /// Resize the I'th \ref vector to the sizes of the I'th \ref feSpaces.
-    template <std::size_t I = 0>
-    void compress(const index_t<I> = {})
-    {
-      std::get<I>(dofvectors).compress();
-    }
-
-    /// Resize the \ref vectors to the sizes of the \ref feSpaces.
-    void compress()
-    {
-      forEach(range_<0, size()>, [this](const auto I) {
-        std::get<I>(dofvectors).compress();
-      });
-    }
-
-    /// Copy the Systemvector \p that element-wise
-    void copy(Self const& that)
-    {
-      forEach(range_<0, size()>, [this, &that](const auto I) {
-        std::get<I>(dofvectors).copy(that[I]);
-      });
-    }
-
-    /// Scale each DOFVector by the factor \p s.
-    template <class Scalar>
-    std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
-    operator*=(Scalar s)
-    {
-      forEach(range_<0, size()>, [this, s](const auto I) {
-        vector_[I] *= s;
-      });
-      return *this;
-    }
-
-    /// Sets each DOFVector to the scalar \p s.
-    template <class Scalar>
-    std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
-    operator=(Scalar s)
-    {
-      forEach(range_<0, size()>, [this, s](const auto I) {
-        vector_[I] = s;
-      });
-      return *this;
-    }
-
-  private:
-    /// a tuple of feSpaces
-    FeSpaces const& feSpaces;
-
-    /// The names of the DOFVectors
-    std::vector<std::string> names;
-
-    /// a tuple of base vectors, i.e. mtl::dense_vectors
-    MultiVector vector_;
-
-    /// a tuple of dofvectors referencing \ref vector
-    DOFVectors dofvectors;
-  };
-
-
-  /// Construct a systemvector from given feSpace-tuple and names-vector
-  template <class FeSpacesType, class ValueType = double>
-  SystemVector<FeSpacesType, ValueType>
-  makeSystemVector(FeSpacesType const& feSpaces, std::vector<std::string> names)
-  {
-    return {feSpaces, names};
-  }
-
-  /// Construct a systemvector from given feSpace-tuple. Create vector of names
-  /// from base_name, i.e. {base_name_0, base_name_1, ...}
-  template <class FeSpaces, class ValueType = double>
-  SystemVector<FeSpaces, ValueType>
-  makeSystemVector(FeSpaces const& feSpaces, std::string base_name)
-  {
-    std::vector<std::string> names;
-    for (std::size_t i = 0; i < std::tuple_size<FeSpaces>::value; ++i)
-      names.push_back(base_name + "_" + std::to_string(i));
-    return {feSpaces, names};
-  }
-
-
-// #ifndef AMDIS_NO_EXTERN_SYSTEMVECTOR
-//     // explicit instantiation in SystemVector.cpp
-//     extern template class SystemVector<typename TestTraits<2,1,2>::FeSpaces>;
-// #endif
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/UmfpackRunner.hpp b/src/amdis/linear_algebra/mtl/UmfpackRunner.hpp
index fbca866f1b4cef556ec22cb702df3c28f21f8eb8..26b0cb0805edca1fa5c349f0268be75264b56a38 100644
--- a/src/amdis/linear_algebra/mtl/UmfpackRunner.hpp
+++ b/src/amdis/linear_algebra/mtl/UmfpackRunner.hpp
@@ -1,23 +1,17 @@
 #pragma once
 
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
 #ifdef HAVE_UMFPACK
 
 #include <algorithm>
-#include <iostream>
 #include <string>
 
 #include <boost/numeric/mtl/operation/two_norm.hpp>
 #include <boost/numeric/mtl/interface/umfpack_solve.hpp>
 
+#include <amdis/Output.hpp>
 #include <amdis/linear_algebra/RunnerInterface.hpp>
 #include <amdis/linear_algebra/SolverInfo.hpp>
 
-#include <amdis/linear_algebra/mtl/Copy.hpp>
-
 namespace AMDiS
 {
   /**
@@ -40,16 +34,16 @@ namespace AMDiS
     using Super = RunnerInterface<Matrix, Vector>;
     using PreconBase = typename Super::PreconBase;
 
-    using FullMatrix = BaseMatrix_t<Matrix>;
+    using FullMatrix = Matrix;
     using SolverType = mtl::mat::umfpack::solver<FullMatrix>;
 
   public:
     /// Constructor. Reads UMFPACK parameters from initfile
     UmfpackRunnerBase(std::string prefix)
     {
-      Parameters::get(prefix + "->store symbolic", store_symbolic); // ?
-      Parameters::get(prefix + "->symmetric strategy", symmetric_strategy);
-      Parameters::get(prefix + "->alloc init", alloc_init);
+      Parameters::get(prefix + "->store symbolic", storeSymbolic_); // ?
+      Parameters::get(prefix + "->symmetric strategy", symmetricStrategy_);
+      Parameters::get(prefix + "->alloc init", allocInit_);
     }
 
     /// Destructor
@@ -60,11 +54,11 @@ namespace AMDiS
                       SolverInfo& solverInfo) override
     {
       AMDIS_FUNCNAME("Umfpack_Runner::solve()");
-      test_exit(bool(solver), "The umfpack solver was not initialized\n");
+      test_exit(bool(solver_), "The umfpack solver was not initialized\n");
 
       int code = 0;
       try {
-	      code = (*solver)(x, b);
+	      code = (*solver_)(x, b);
       } catch (mtl::mat::umfpack::error& e) {
 	      error_exit("UMFPACK_ERROR(solve, ", e.code, ") = ", e.what());
       }
@@ -79,47 +73,17 @@ namespace AMDiS
       return code;
     }
 
-    /// Implementation of \ref RunnerInterface::adjoint_solve()
+    /// Implementation of \ref RunnerInterface::exit()
     virtual void exit() override {}
 
   protected:
-    std::shared_ptr<SolverType> solver;
+    std::shared_ptr<SolverType> solver_;
 
-    int store_symbolic = 0;
-    int symmetric_strategy = 0;
-    double alloc_init = 0.7;
+    int storeSymbolic_ = 0;
+    int symmetricStrategy_ = 0;
+    double allocInit_ = 0.7;
   };
 
-  // specialization for block-matrices
-  template <class SubMatrix, std::size_t _N, std::size_t _M, class Vector>
-  class UmfpackRunner<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>
-      : public UmfpackRunnerBase<BlockMTLMatrix<SubMatrix, _N, _M>, Vector>
-  {
-    using Matrix = BlockMTLMatrix<SubMatrix, _N, _M>;
-    using Super = UmfpackRunnerBase<Matrix, Vector>;
-
-    using SolverType = typename Super::SolverType;
-
-  public:
-    UmfpackRunner(std::string prefix)
-      : Super(prefix)
-    {}
-
-    /// Implementation of \ref RunnerInterface::init()
-    virtual void init(Matrix const& A) override
-    {
-      copy(A, fullMatrix);
-
-      try {
-        Super::solver.reset(new SolverType(fullMatrix, Super::symmetric_strategy, Super::alloc_init));
-      } catch (mtl::mat::umfpack::error const& e) {
-        error_exit("UMFPACK_ERROR(factorize, ", e.code, ") = ", e.what());
-      }
-    }
-
-  private:
-    typename Super::FullMatrix fullMatrix;
-  };
 
   // specialization for full-matrices
   template <class T, class Param, class Vector>
@@ -140,7 +104,7 @@ namespace AMDiS
     virtual void init(FullMatrix const& fullMatrix) override
     {
       try {
-        Super::solver.reset(new SolverType(fullMatrix, Super::symmetric_strategy, Super::alloc_init));
+        Super::solver_.reset(new SolverType(fullMatrix, Super::symmetricStrategy_, Super::allocInit_));
       } catch (mtl::mat::umfpack::error const& e) {
         error_exit("UMFPACK_ERROR(factorize, ", e.code, ") = ", e.what());
       }
diff --git a/src/amdis/utility/MultiIndex.hpp b/src/amdis/utility/MultiIndex.hpp
index 196f12d70829ba7f3abbbfb24f72ce625e9e02b4..424e70aedd8da13f50bd1f5b780525d1e733b37a 100644
--- a/src/amdis/utility/MultiIndex.hpp
+++ b/src/amdis/utility/MultiIndex.hpp
@@ -1,17 +1,23 @@
 #pragma once
 
 #include <dune/common/reservedvector.hh>
-#include <dune/functions/common/indexaccess.hh>
 #include <dune/functions/functionspacebases/flatmultiindex.hh>
 
-#include <amdis/common/Mpl.hpp>
-
 namespace AMDiS
 {
-  std::size_t flatMultiIndex(std::size_t idx) { return idx; }
+  std::size_t flatMultiIndex(std::size_t idx)
+  {
+    return idx;
+  }
 
-  std::size_t flatMultiIndex(Dune::Functions::FlatMultiIndex<std::size_t> const& idx) { return idx[0]; }
+  std::size_t flatMultiIndex(Dune::Functions::FlatMultiIndex<std::size_t> const& idx)
+  {
+    return idx[0];
+  }
 
-  std::size_t flatMultiIndex(Dune::ReservedVector<std::size_t,1> const& idx) { return idx[0]; }
+  std::size_t flatMultiIndex(Dune::ReservedVector<std::size_t,1> const& idx)
+  {
+    return idx[0];
+  }
 
 } // end namespace AMDiS