diff --git a/examples/ellipt.cc b/examples/ellipt.cc
index fdfc81c6a6bd1f85e41aa29816f8e4bfe2e5d65d..fc7830e8fee7e0b9b46ec87198a53e30593fb627 100644
--- a/examples/ellipt.cc
+++ b/examples/ellipt.cc
@@ -15,7 +15,7 @@ using namespace AMDiS;
 using namespace Dune::Indices;
 
 // 1 component with polynomial degree 1
-using Param   = YaspGridBasis<AMDIS_DIM, 1>;
+using Param   = YaspGridBasis<AMDIS_DIM, 2>;
 using ElliptProblem = ProblemStat<Param>;
 
 int main(int argc, char** argv)
diff --git a/examples/init/ellipt.dat.2d b/examples/init/ellipt.dat.2d
index f12565729e95e06f7902cf81744b29450c3445e1..0ef4d635178cbd825b31cbd818eab17a036d5c72 100644
--- a/examples/init/ellipt.dat.2d
+++ b/examples/init/ellipt.dat.2d
@@ -3,10 +3,10 @@ elliptMesh->global refinements: 0
 ellipt->mesh:                   elliptMesh
 
 ellipt->solver->name:           bicgstab
-ellipt->solver->max iteration:  1000
-ellipt->solver->relative tolerance: 1.e-8
-ellipt->solver->info:           10
-ellipt->solver->left precon:    diag
+ellipt->solver->max iteration:  10000
+ellipt->solver->relative tolerance: 1.e-10
+ellipt->solver->info:           1
+ellipt->solver->left precon:    ilu
 ellipt->solver->right precon:   no
 
 ellipt->output[0]->filename:    ellipt.2d
diff --git a/examples/stokes0.cc b/examples/stokes0.cc
index e3785be1aa1dab3ed9b9cad3622b90a73b5fdedd..40a8c425474d611be7ca3d5845f55270f2055eab 100644
--- a/examples/stokes0.cc
+++ b/examples/stokes0.cc
@@ -80,10 +80,12 @@ int main(int argc, char** argv)
     // assemble and solve system
     prob.buildAfterCoarsen(adaptInfo, Flag(0));
 
+#ifdef DEBUG_MTL
     // write matrix to file
     mtl::io::matrix_market_ostream out("matrix_stokes0.mtx");
     out << prob.getSystemMatrix().matrix();
     std::cout << prob.getSystemMatrix().matrix() << '\n';
+#endif
 
     prob.solve(adaptInfo);
 
diff --git a/examples/stokes1.cc b/examples/stokes1.cc
index 9126c3d406d904c9b9d65cbd59ff3a8213195af0..bdb95a3b5c11de137c2ba754850a7c708ad38fc0 100644
--- a/examples/stokes1.cc
+++ b/examples/stokes1.cc
@@ -81,10 +81,12 @@ int main(int argc, char** argv)
     // assemble and solve system
     prob.buildAfterCoarsen(adaptInfo, Flag(0));
 
+#ifdef DEBUG_MTL
     // write matrix to file
     mtl::io::matrix_market_ostream out("matrix_stokes1.mtx");
     out << prob.getSystemMatrix().matrix();
     std::cout << prob.getSystemMatrix().matrix() << '\n';
+#endif
 
     prob.solve(adaptInfo);
 
diff --git a/examples/vecellipt.cc b/examples/vecellipt.cc
index f078cdd724cc2f88102bcc7ed5c3a005f19e38e9..a4dfb40735a29230dbad1e02bedb43c9e4cf8927 100644
--- a/examples/vecellipt.cc
+++ b/examples/vecellipt.cc
@@ -42,6 +42,7 @@ int main(int argc, char** argv)
 
   prob.buildAfterCoarsen(adaptInfo, Flag(0));
 
+#ifdef DEBUG_MTL
   // write matrix to file
   if (Parameters::get<int>("elliptMesh->global refinements").value_or(0) < 4) {
     mtl::io::matrix_market_ostream out("matrix.mtx");
@@ -49,6 +50,7 @@ int main(int argc, char** argv)
 
     std::cout << prob.getSystemMatrix().matrix() << '\n';
   }
+#endif
 
   prob.solve(adaptInfo);
   prob.writeFiles(adaptInfo, true);
diff --git a/src/amdis/Assembler.inc.hpp b/src/amdis/Assembler.inc.hpp
index 2d4610a827bdaf5b31b120481ae8cb5add4dc52c..ac0ccb050f6ac3ace69290cae99331a663c0b948 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();
 
-  mtl::mat::dense2D<typename SystemMatrixType::value_type> elementMatrix(localSize, localSize);
-  mtl::vec::dense_vector<typename SystemVectorType::value_type> elementVector(localSize);
+  mtl::mat::dense2D<double> elementMatrix(localSize, localSize);
+  mtl::vec::dense_vector<double> elementVector(localSize);
 
   // 3. traverse grid and assemble operators on the elements
   for (auto const& element : elements(globalBasis_.gridView()))
diff --git a/src/amdis/LinearAlgebra.hpp b/src/amdis/LinearAlgebra.hpp
index 898404a001e91ab8cfd2df5f37f9f1f781f45b0f..d0ec769f15b0ce3d8c9b72a40afc026009218e0b 100644
--- a/src/amdis/LinearAlgebra.hpp
+++ b/src/amdis/LinearAlgebra.hpp
@@ -1,8 +1,16 @@
 #pragma once
 
-#include <amdis/linear_algebra/LinearSolverInterface.hpp>
+#include <amdis/linear_algebra/LinearSolver.hpp>
 #include <amdis/linear_algebra/SolverInfo.hpp>
+
+#if HAVE_MTL
 #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/ITL_Preconditioner.hpp>
+#else
+#include <amdis/linear_algebra/istl/DOFVector.hpp>
+#include <amdis/linear_algebra/istl/DOFMatrix.hpp>
+#include <amdis/linear_algebra/istl/ISTL_Solver.hpp>
+#include <amdis/linear_algebra/istl/ISTL_Preconditioner.hpp>
+#endif
\ No newline at end of file
diff --git a/src/amdis/linear_algebra/DOFMatrixBase.hpp b/src/amdis/linear_algebra/DOFMatrixBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e7bf16d923776bf233d6114d4706986d9929d9b3
--- /dev/null
+++ b/src/amdis/linear_algebra/DOFMatrixBase.hpp
@@ -0,0 +1,14 @@
+#pragma once
+
+#include <cstddef>
+
+namespace AMDiS
+{
+  template <class T>
+  struct Triplet
+  {
+    std::size_t row, cols;
+    T value;
+  };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/DOFVectorBase.hpp b/src/amdis/linear_algebra/DOFVectorBase.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..68cb1c990c2cb8555db5b40571101ade7df8eb38
--- /dev/null
+++ b/src/amdis/linear_algebra/DOFVectorBase.hpp
@@ -0,0 +1,53 @@
+#pragma once
+
+#include <algorithm>
+#include <memory>
+#include <string>
+
+#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/DOFVectorInterface.hpp>
+#include <amdis/utility/MultiIndex.hpp>
+
+namespace AMDiS
+{
+  /// The basic container that stores a base vector and a corresponding basis
+  template <class BasisType>
+  class DOFVectorBase
+      : public DOFVectorInterface
+  {
+    using Self = DOFVectorBase;
+
+  public:
+    /// The type of the \ref basis
+    using Basis = BasisType;
+
+    /// The index/size - type
+    using size_type  = typename BasisType::size_type;
+
+    /// Constructor. Constructs new BaseVector.
+    DOFVectorBase(BasisType const& basis)
+      : basis_(&basis)
+    {}
+
+    /// Return the basis \ref basis_ of the vector
+    Basis const& basis() const
+    {
+        return *basis_;
+    }
+
+    /// Return the size of the \ref basis
+    size_type size() const
+    {
+      return basis_->dimension();
+    }
+
+  private:
+    /// The finite element space / basis associated with the data vector
+    Basis const* basis_;
+  };
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/LinearSolver.hpp b/src/amdis/linear_algebra/LinearSolver.hpp
similarity index 72%
rename from src/amdis/linear_algebra/mtl/LinearSolver.hpp
rename to src/amdis/linear_algebra/LinearSolver.hpp
index ed44e4c2ea06d0d2d3d20715e28cd85ae1074a2c..b132ef31ba780e5591000c50bf304a5848b29837 100644
--- a/src/amdis/linear_algebra/mtl/LinearSolver.hpp
+++ b/src/amdis/linear_algebra/LinearSolver.hpp
@@ -14,18 +14,18 @@ namespace AMDiS
   /** \ingroup Solver
    *
    * \brief
-   * Wrapper class for various MTL4 solvers. These solvers
+   * Wrapper class for various solvers. These solvers
    * are parametrized by MatrixType and VectorType. The different
    * algorithms, like krylov subspace methods or other external
-   * solvers where MTL4 provides an interface, can be assigned
+   * solvers where the backend provides an interface, can be assigned
    * by different Runner objects.
    **/
-  template <class Matrix, class Vector, class Runner>
+  template <class Matrix, class VectorX, class VectorB, class Runner>
   class LinearSolver
-      : public LinearSolverInterface<Matrix, Vector, Vector>
+      : public LinearSolverInterface<Matrix, VectorX, VectorB>
   {
     using Self = LinearSolver;
-    using Super = LinearSolverInterface<Matrix, Vector, Vector>;
+    using Super = LinearSolverInterface<Matrix, VectorX, VectorB>;
     using RunnerBase = typename Super::RunnerBase;
 
   public:
@@ -44,15 +44,15 @@ namespace AMDiS
       : runner_(std::make_shared<Runner>(prefix))
     {}
 
-    /// Implements \ref LinearSolverInterface::getRunner()
-    virtual std::shared_ptr<RunnerBase> getRunner() override
+    /// Implements \ref LinearSolverInterface::runner()
+    virtual std::shared_ptr<RunnerBase> runner() override
     {
       return runner_;
     }
 
   private:
     /// Implements \ref LinearSolverInterface::solveSystemImpl()
-    virtual void solveImpl(Matrix const& A, Vector& x, Vector const& b,
+    virtual void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
                            SolverInfo& solverInfo) override
     {
       Dune::Timer t;
@@ -61,9 +61,8 @@ namespace AMDiS
         runner_->init(A);
       }
 
-      if (solverInfo.getInfo() > 0) {
-        msg("fill MTL4 matrix needed ", t.elapsed(), " seconds");
-      }
+      if (solverInfo.getInfo() > 0)
+        msg("fill matrix needed ", t.elapsed(), " seconds");
 
       int error = runner_->solve(A, x, b, solverInfo);
       solverInfo.setError(error);
diff --git a/src/amdis/linear_algebra/LinearSolverInterface.hpp b/src/amdis/linear_algebra/LinearSolverInterface.hpp
index 3edf1d2dedee74c92922d5179f5a43b7d35d4290..da59c028bd41f123edc659603917d38b74efb71f 100644
--- a/src/amdis/linear_algebra/LinearSolverInterface.hpp
+++ b/src/amdis/linear_algebra/LinearSolverInterface.hpp
@@ -48,7 +48,7 @@ namespace AMDiS
     }
 
     // return the runner/worker corresponding to this solver (optional)
-    virtual std::shared_ptr<RunnerBase> getRunner() { return {}; };
+    virtual std::shared_ptr<RunnerBase> runner() { return {}; };
 
   private:
     /// main methods that all solvers must implement
diff --git a/src/amdis/linear_algebra/istl/CMakeLists.txt b/src/amdis/linear_algebra/istl/CMakeLists.txt
index cd5ee97f5f7fbaf3fa22c39f74a54b6eeb33f9ce..26238d75de2f02d584150a6e059b496557369d57 100644
--- a/src/amdis/linear_algebra/istl/CMakeLists.txt
+++ b/src/amdis/linear_algebra/istl/CMakeLists.txt
@@ -11,7 +11,4 @@ install(FILES
     ISTLRunner.hpp
     LinearSolver.hpp
     Preconditioner.hpp
-    SystemMatrix.hpp
-    SystemVector.hpp
-    UmfpackRunner.hpp
 DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/amdis/linear_algebra/istl)
diff --git a/src/amdis/linear_algebra/istl/DOFMatrix.hpp b/src/amdis/linear_algebra/istl/DOFMatrix.hpp
index fbaaa6357c4857225b4febac361c306c1d7d9185..09739652d3b2649e12946e69876bc6ce5b64428e 100644
--- a/src/amdis/linear_algebra/istl/DOFMatrix.hpp
+++ b/src/amdis/linear_algebra/istl/DOFMatrix.hpp
@@ -8,111 +8,126 @@
 
 #include <amdis/Output.hpp>
 #include <amdis/common/ClonablePtr.hpp>
+#include <amdis/linear_algebra/DOFMatrixBase.hpp>
+#include <amdis/utility/MultiIndex.hpp>
 
 namespace AMDiS
 {
-  template <class RowFeSpaceType, class ColFeSpaceType,
-            class ValueType = Dune::FieldMatrix<double,1,1>>
+  template <class T, class = void>
+  struct BlockMatrixType
+  {
+    using type = Dune::FieldMatrix<T,1,1>;
+  };
+
+  template <class T>
+  struct BlockMatrixType<T, typename T::field_type>
+  {
+    using type = T;
+  };
+
+  template <class RowBasisType, class ColBasisType,
+            class ValueType = double>
   class DOFMatrix
   {
   public:
-    using RowFeSpace = RowFeSpaceType;
-    using ColFeSpace = ColFeSpaceType;
-    using BaseMatrix = Dune::BCRSMatrix<ValueType>;
+    /// The type of the finite element space / basis of the row
+    using RowBasis = RowBasisType;
+
+    /// The type of the finite element space / basis of the column
+    using ColBasis = ColBasisType;
 
-    using size_type  = typename RowFeSpaceType::size_type;
-    using field_type = typename ValueType::field_type;
+    /// The type of the elements of the DOFMatrix
+    using value_type = typename BlockMatrixType<ValueType>::type;
 
-    using value_type = ValueType;
+    /// The matrix type of the underlying base matrix
+    using BaseMatrix = Dune::BCRSMatrix<value_type>;
 
+    /// The index/size - type
+    using size_type = typename BaseMatrix::size_type;
+
+  public:
     /// Constructor. Constructs new BaseVector.
-    DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name)
-      : rowFeSpace(rowFeSpace)
-      , colFeSpace(colFeSpace)
-      , name(name)
-      , matrix(ClonablePtr<BaseMatrix>::make())
+    DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis)
+      : rowBasis_(&rowBasis)
+      , colBasis_(&colBasis)
+      , matrix_(ClonablePtr<BaseMatrix>::make())
     {}
 
-    /// Constructor. Takes pointer of data-vector.
-    DOFMatrix(RowFeSpace const& rowFeSpace, ColFeSpace const& colFeSpace, std::string name,
-              BaseMatrix& matrix_ref)
-      : rowFeSpace(rowFeSpace)
-      , colFeSpace(colFeSpace)
-      , name(name)
-      , matrix(matrix_ref)
+    /// Constructor. Takes reference to a base matrix.
+    DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis, BaseMatrix& matrix)
+      : rowBasis_(&rowBasis)
+      , colBasis_(&colBasis)
+      , matrix_(matrix)
     {}
 
-    /// Return the row-basis \ref rowFeSpace of the matrix
-    RowFeSpace const& getRowFeSpace() const
+    /// Return the row-basis \ref rowBasis of the matrix
+    RowBasis const& rowBasis() const
     {
-      return rowFeSpace;
+      return *rowBasis_;
     }
 
-    /// Return the col-basis \ref colFeSpace of the matrix
-    ColFeSpace const& getColFeSpace() const
+    /// Return the col-basis \ref colBasis of the matrix
+    ColBasis const& colBasis() const
     {
-      return colFeSpace;
+      return *colBasis_;
     }
 
     /// Return the data-vector \ref vector
-    BaseMatrix const& getMatrix() const
+    BaseMatrix const& matrix() const
     {
-      return *matrix;
+      return *matrix_;
     }
 
     /// Return the data-vector \ref vector
-    BaseMatrix& getMatrix()
+    BaseMatrix& matrix()
     {
-      return *matrix;
+      return *matrix_;
     }
 
     /// Return the size of the \ref feSpace
-    size_type N() const
-    {
-      return rowFeSpace.size();
-    }
-
-    size_type M() const
+    size_type rows() const
     {
-      return colFeSpace.size();
+      return rowBasis_->size();
     }
 
-    /// Return the \ref name of this vector
-    std::string getName() const
+    size_type cols() const
     {
-      return name;
+      return colBasis_->size();
     }
 
 
     /// Access the entry \p i of the \ref vector with read-access.
-    value_type const& operator()(size_type r, size_type c) const
+    template <class Index>
+    value_type const& operator()(Index row, Index col) const
     {
-      test_exit_dbg( initialized, "Occupation pattern not initialized!");
-      test_exit_dbg( r < N() && c < M() ,
-          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
-      return (*matrix)[r][c];
+      size_type r = flatMultiIndex(row), c = flatMultiIndex(col);
+      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
+      test_exit_dbg( r < rows() && c < cols() ,
+          "Indices out of range [0,", rows(), ")x[0,", cols(), ")" );
+      return (*matrix_)[r][c];
     }
 
     /// Access the entry \p i of the \ref vector with write-access.
-    value_type& operator()(size_type r, size_type c)
+    template <class Index>
+    value_type& operator()(Index row, Index col)
     {
-      test_exit_dbg( initialized, "Occupation pattern not initialized!");
-      test_exit_dbg( r < N() && c < M() ,
-          "Indices out of range [0,", N(), ")x[0,", M(), ")" );
-      return (*matrix)[r][c];
+      size_type r = flatMultiIndex(row), c = flatMultiIndex(col);
+      test_exit_dbg( initialized_, "Occupation pattern not initialized!");
+      test_exit_dbg( r < rows() && c < cols() ,
+          "Indices out of range [0,", rows(), ")x[0,", cols(), ")" );
+      return (*matrix_)[r][c];
     }
 
     /// create occupation pattern and apply it to the matrix
-    void init()
+    void init(bool prepareForInsertion)
     {
-      occupationPattern.resize(rowFeSpace.size(), colFeSpace.size());
-      auto meshView = rowFeSpace.gridView();
+      occupationPattern_.resize(rowBasis_->size(), colBasis_->size());
 
       // A loop over all elements of the grid
-      auto rowLocalView = rowFeSpace.localView();
-      auto colLocalView = colFeSpace.localView();
+      auto rowLocalView = rowBasis_->localView();
+      auto colLocalView = colBasis_->localView();
 
-      for (const auto& element : elements(meshView)) {
+      for (const auto& element : elements(rowBasis_->gridView())) {
         rowLocalView.bind(element);
         colLocalView.bind(element);
 
@@ -122,54 +137,50 @@ namespace AMDiS
           for (std::size_t j = 0; j < colLocalView.size(); ++j) {
             // The global index of the j-th vertex of the element
             auto col = colLocalView.index(j);
-            occupationPattern.add(row, col);
+            occupationPattern_.add(row, col);
           }
         }
       }
-      occupationPattern.exportIdx(*matrix);
+      occupationPattern_.exportIdx(*matrix_);
 
-      initialized = true;
+      initialized_ = true;
     }
 
     void finish()
     {
-      initialized = false;
+      initialized_ = false;
     }
 
-    std::size_t nnz()
+    std::size_t nnz() const
     {
-      return matrix->nonzeroes();
+      return matrix_->nonzeroes();
     }
 
-    auto clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
+    auto applyDirichletBC(std::vector<char> const& dirichletNodes)
     {
       // loop over the matrix rows
-      for (std::size_t i = 0; i < matrix->N(); ++i) {
+      for (std::size_t i = 0; i < matrix_->N(); ++i) {
         if (dirichletNodes[i]) {
-            auto cIt = (*matrix)[i].begin();
-            auto cEndIt = (*matrix)[i].end();
+            auto cIt = (*matrix_)[i].begin();
+            auto cEndIt = (*matrix_)[i].end();
             // loop over nonzero matrix entries in current row
             for (; cIt != cEndIt; ++cIt)
-                *cIt = (apply && i == cIt.index() ? 1 : 0);
+                *cIt = (i == cIt.index() ? 1 : 0);
         }
       }
 
-      std::vector<std::map<size_t,value_type>> columns; // symmetric dbc not yet implemented
+      std::list<Triplet<value_type>> columns; // symmetric dbc not yet implemented
       return columns;
     }
 
   private:
-    RowFeSpace const&	rowFeSpace;
-    ColFeSpace const&	colFeSpace;
-    std::string	        name;
-
-    ClonablePtr<BaseMatrix> matrix;
-    Dune::MatrixIndexSet    occupationPattern;
+    RowBasis const*	rowBasis_;
+    ColBasis const*	colBasis_;
 
-    bool initialized = false;
+    ClonablePtr<BaseMatrix> matrix_;
+    Dune::MatrixIndexSet occupationPattern_;
 
-    // friend class declarations
-    template <class, class, class> friend class SystemMatrix;
+    bool initialized_ = false;
   };
 
 } // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/DOFVector.hpp b/src/amdis/linear_algebra/istl/DOFVector.hpp
index 63a5aebb576426830f9354ed8b3456044f5c93b3..ddcfd1516bb746e92feb275bf58c0b9217950933 100644
--- a/src/amdis/linear_algebra/istl/DOFVector.hpp
+++ b/src/amdis/linear_algebra/istl/DOFVector.hpp
@@ -1,124 +1,126 @@
 #pragma once
 
-#include <string>
-#include <memory>
-
 #include <dune/istl/bvector.hh>
 
-#include <dune/functions/functionspacebases/interpolate.hh>
-
 #include <amdis/Output.hpp>
 #include <amdis/common/ClonablePtr.hpp>
+#include <amdis/linear_algebra/DOFVectorBase.hpp>
+#include <amdis/utility/MultiIndex.hpp>
 
 namespace AMDiS
 {
-    template <class FeSpaceType, class ValueType = Dune::FieldVector<double,1>>
-    class DOFVector
-    {
-        using Self = DOFVector;
+  template <class T, class = void>
+  struct BlockVectorType
+  {
+    using type = Dune::FieldVector<T,1>;
+  };
+
+  template <class T>
+  struct BlockVectorType<T, typename T::field_type>
+  {
+    using type = T;
+  };
+
+	template <class BasisType, class ValueType = double>
+	class DOFVector
+      : public DOFVectorBase<BasisType>
+	{
+	  using Self = DOFVector;
+    using Super = DOFVectorBase<BasisType>;
 
-    public:
-	using FeSpace    = FeSpaceType;
-	using BaseVector = Dune::BlockVector<ValueType>;
+	public:
+    /// The type of the finite element space / basis
+	  using Basis = BasisType;
 
-	using size_type  = typename FeSpace::size_type;
-	using field_type = typename ValueType::field_type;
+    /// The type of the elements of the DOFVector
+	  using value_type = typename BlockVectorType<ValueType>::type;
 
-	using value_type = ValueType;
+    /// The vector type of the underlying base vector
+	  using BaseVector = Dune::BlockVector<value_type>;
 
-	/// Constructor. Constructs new BaseVector.
-	DOFVector(FeSpace const& feSpace, std::string name)
-	  : feSpace(feSpace)
-	  , name(name)
-	  , vector(ClonablePtr<BaseVector>::make())
-	{
-	    compress();
-	}
-
-	/// Constructor. Takes pointer of data-vector.
-	DOFVector(FeSpace const& feSpace, std::string name,
-		  BaseVector& vector_ref)
-	  : feSpace(feSpace)
-	  , name(name)
-	  , vector(vector_ref)
-	{}
-
-	/// Return the basis \ref feSpace of the vector
-	FeSpace const& getFeSpace() const
-	{
-	    return feSpace;
-	}
+    /// The index/size - type
+	  using size_type  = typename Basis::size_type;
 
-	/// Return the data-vector \ref vector
-	BaseVector const& vector() const
-	{
-	    return *vector;
-	}
-
-	/// Return the data-vector \ref vector
-	BaseVector& vector()
-	{
-	    return *vector;
-	}
-
-	/// Return the size of the \ref feSpace
-	size_type getSize() const
-	{
-	    return feSpace.size();
-	}
+    /// Constructor. Constructs new BaseVector.
+    DOFVector(BasisType const& basis)
+      : Super(basis)
+      , vector_(ClonablePtr<BaseVector>::make())
+    {
+      compress();
+    }
 
-	/// Return the \ref name of this vector
-	std::string getName() const
-	{
-	    return name;
-	}
+    /// Constructor. Takes pointer of data-vector.
+    DOFVector(BasisType const& basis, Dune::BlockVector<ValueType>& vector)
+      : Super(basis)
+      , vector_(vector)
+    {}
 
-	/// Resize the \ref vector to the size of the \ref feSpace.
-	void compress()
-	{
-	    vector->resize(getSize() / value_type::dimension);
-	}
+    /// Return the data-vector \ref vector
+    BaseVector const& vector() const
+    {
+      return *vector_;
+    }
 
+    /// Return the data-vector \ref vector
+    BaseVector& vector()
+    {
+      return *vector_;
+    }
 
-	/// Access the entry \p i of the \ref vector with read-access.
-	value_type const& operator[](size_type i) const
-	{
-	    test_exit_dbg(i < vector->size() ,
-                          "Index ", i, " out of range [0,", vector->size(), ")" );
-	    return (*vector)[i];
-	}
 
-	/// Access the entry \p i of the \ref vector with write-access.
-	value_type& operator[](size_type i)
-	{
-	    test_exit_dbg(i < vector->size() ,
-                          "Index ", i, " out of range [0,", vector->size(), ")" );
-	    return (*vector)[i];
-	}
-
-	/// interpolate a function \p f to the basis \ref feSpace and store the
-	/// coefficients in \ref vector.
-	template <class F>
-	void interpol(F const& f)
-	{
-	    Dune::Functions::interpolate(feSpace, *vector, f);
-	}
+    /// Resize the \ref vector to the size of the \ref basis.
+    virtual void compress() override
+    {
+      if (vector_->size() != this->size()) {
+        resize(this->size());
+        *vector_ = value_type(0);
+      }
+    }
+
+    template <class SizeInfo>
+    void resize(SizeInfo s)
+    {
+      vector_->resize(size_type(s));
+    }
 
-	/// Calls the copy assignment operator of the BaseVector \ref vector
-	void copy(Self const& that)
-	{
-		*vector = that.vector();
-	}
 
-    private:
-	FeSpace const&	feSpace;
-	std::string	name;
+    /// Access the entry \p i of the \ref vector with read-access.
+    template <class Index>
+    value_type const& operator[](Index idx) const
+    {
+      size_type i = flatMultiIndex(idx);
+      test_exit_dbg(i < vector_->size(),
+        "Index ", i, " out of range [0,", vector_->size(), ")" );
+      return (*vector_)[i];
+    }
+
+    /// Access the entry \p i of the \ref vector with write-access.
+    template <class Index>
+    value_type& operator[](Index idx)
+    {
+      size_type i = flatMultiIndex(idx);
+      test_exit_dbg(i < vector_->size(),
+        "Index ", i, " out of range [0,", vector_->size(), ")" );
+      return (*vector_)[i];
+    }
+
+    /// Sets each DOFVector to the scalar \p s.
+    template <class Scalar>
+    std::enable_if_t< Concepts::Arithmetic<Scalar>, Self&>
+    operator=(Scalar s)
+    {
+      (*vector_) = s;
+      return *this;
+    }
 
-	ClonablePtr<BaseVector> vector;
+    /// Calls the copy assignment operator of the BaseVector \ref vector
+    void copy(Self const& that)
+    {
+      *vector_ = that.vector();
+    }
 
-    // friend class declarations
-    template <class, class>
-    friend class SystemVector;
-    };
+	private:
+    ClonablePtr<BaseVector> vector_;
+	};
 
 } // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/DirectRunner.hpp b/src/amdis/linear_algebra/istl/DirectRunner.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..187abadfa5c0767810abf3147a4c28d8ef1e0733
--- /dev/null
+++ b/src/amdis/linear_algebra/istl/DirectRunner.hpp
@@ -0,0 +1,63 @@
+#pragma once
+
+#include <algorithm>
+#include <string>
+
+#include <amdis/Output.hpp>
+#include <amdis/linear_algebra/RunnerInterface.hpp>
+#include <amdis/linear_algebra/SolverInfo.hpp>
+
+namespace AMDiS
+{
+  /**
+   * \ingroup Solver
+   * \class AMDiS::DirectRunner
+   * \brief \implements RunnerInterface for the (external) direct solvers
+   */
+  template <class Solver, class Matrix, class VectorX, class VectorB>
+  class DirectRunner
+      : public RunnerInterface<Matrix, VectorX, VectorB>
+  {
+  protected:
+    using Super = RunnerInterface<Matrix, VectorX, VectorB>;
+    using BaseSolver  = Dune::InverseOperator<VectorX, VectorB>;
+
+  public:
+    /// Constructor.
+    DirectRunner(std::string prefix)
+      : solverCreator_(prefix)
+    {}
+
+    /// Implementes \ref RunnerInterface::init()
+    virtual void init(Matrix const& A) override
+    {
+      solver_ = solverCreator_.create(A);
+    }
+
+    /// Implementes \ref RunnerInterface::exit()
+    virtual void exit() override
+    {
+      solver_.reset();
+    }
+
+    /// Implementes \ref RunnerInterface::solve()
+    virtual int solve(Matrix const& A, VectorX& x, VectorB const& b,
+                      SolverInfo& solverInfo) override
+    {
+      // storing some statistics
+      Dune::InverseOperatorResult statistics;
+
+      // solve the linear system
+      VectorB _b = b;
+      solver_->apply(x, _b, statistics);
+
+      solverInfo.setRelResidual(statistics.reduction);
+
+      return 0;
+    }
+
+  private:
+    ISTLSolverCreator<Solver> solverCreator_;
+    std::shared_ptr<BaseSolver> solver_;
+  };
+}
diff --git a/src/amdis/linear_algebra/istl/Fwd.hpp b/src/amdis/linear_algebra/istl/Fwd.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e2ad22c295e1f6de8e8926c94be5b7bf36fc7f9c
--- /dev/null
+++ b/src/amdis/linear_algebra/istl/Fwd.hpp
@@ -0,0 +1,8 @@
+#pragma once
+
+namespace AMDiS
+{
+  template <class ISTLSolver, bool specialized = true>
+  struct ISTLSolverCreator;
+
+} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/ISTLRunner.hpp b/src/amdis/linear_algebra/istl/ISTLRunner.hpp
index 46d7790fed6db10f62cb510ea65b0e24f200e266..99c4927e876456f1bfffc9c2586467c8ec6fb69e 100644
--- a/src/amdis/linear_algebra/istl/ISTLRunner.hpp
+++ b/src/amdis/linear_algebra/istl/ISTLRunner.hpp
@@ -1,56 +1,45 @@
 #pragma once
 
-#include <amdis/linear_algebra/istl/ISTL_Solver.hpp>
-#include <amdis/linear_algebra/istl/Preconditioner.hpp>
+#include <dune/istl/preconditioner.hh>
+
+#include <amdis/linear_algebra/RunnerInterface.hpp>
+#include <amdis/linear_algebra/SolverInfo.hpp>
+#include <amdis/linear_algebra/istl/Fwd.hpp>
+#include <amdis/linear_algebra/istl/ISTL_Preconditioner.hpp>
 
 namespace AMDiS
 {
-  template <class Matrix, class VectorX, class VectorB>
+  template <class ISTLSolver, class Matrix, class VectorX, class VectorB>
   class ISTLRunner
-    : public RunnerInterface<Matrix, VectorX, VectorB>
+      : public RunnerInterface<Matrix, VectorX, VectorB>
   {
     using LinOperator = Dune::MatrixAdapter<Matrix, VectorX, VectorB>;
     using BaseSolver  = Dune::InverseOperator<VectorX, VectorB>;
-
-    using PreconRunner = ISTLPreconditioner<Matrix, VectorX, VectorB>;
-    using PreconBase   = PreconditionerInterface<Matrix, VectorX, VectorB>;
+    using BasePrecon  = Dune::Preconditioner<VectorX, VectorB>;
+    using ISTLPrecon  = ISTLPreconInterface<Matrix, VectorX, VectorB>;
 
   public:
-    ISTLRunner(std::string solverName, std::string prefix)
-      : solverName(solverName)
-      , prefix(prefix)
-      , solverCreator(prefix)
+    ISTLRunner(std::string prefix)
+      : solverCreator_(prefix)
     {
-      // get creator string for preconditioner
-      std::string preconName = "no";
-      Parameters::get(prefix + "->left precon", preconName);
-      if (preconName.empty() || preconName == "no")
-        Parameters::get(prefix + "->right precon", preconName);
-      if (preconName.empty() || preconName == "no")
-        Parameters::get(prefix + "->precon->name", preconName);
-
-      precon = std::make_shared<PreconRunner>(preconName, prefix + "->precon");
+      initPrecon(prefix);
     }
 
 
     /// Implementes \ref RunnerInterface::init()
     virtual void init(Matrix const& A) override
     {
-      precon->init(A);
-
-      Matrix& _A = const_cast<Matrix&>(A); // Unschoen!!!
-      linOperator = std::make_shared<LinOperator>(_A);
-      solver = solverCreator.create(solverName, *linOperator, *precon);
+      precon_ = preconCreator_->create(A);
+      linOperator_ = std::make_shared<LinOperator>(A);
+      solver_ = solverCreator_.create(*linOperator_, *precon_);
     }
 
     /// Implementes \ref RunnerInterface::exit()
     virtual void exit() override
     {
-      solver.reset();
-      linOperator.reset();
-
-      precon->exit();
-      precon.reset();
+      solver_.reset();
+      linOperator_.reset();
+      precon_.reset();
     }
 
     /// Implementes \ref RunnerInterface::solve()
@@ -62,26 +51,42 @@ namespace AMDiS
 
       // solve the linear system
       VectorB _b = b;
-      solver->apply(x, _b, statistics);
+      solver_->apply(x, _b, statistics);
 
       solverInfo.setRelResidual(statistics.reduction);
+
+      return statistics.converged ? 0 : 1;
     }
 
-    /// Implementes \ref RunnerInterface::getLeftPrecon()
-    virtual std::shared_ptr<PreconBase> getLeftPrecon() override
+  protected:
+    void initPrecon(std::string const& prefix)
     {
-      return preconAdapter;
+      // get creator string for preconditioner
+      std::string preconName = "no";
+      std::string initFileStr = "";
+      for (std::string postfix : {"left precon", "right precon", "precon->name"}) {
+        Parameters::get(prefix + "->" + postfix, preconName);
+        if (!preconName.empty() && preconName != "no") {
+          initFileStr = prefix + "->" + postfix;
+          break;
+        }
+      }
+
+      if (preconName.empty() || preconName == "no")
+        preconName = "default";
+
+      auto* creator = named(CreatorMap<ISTLPrecon>::getCreator(preconName, initFileStr));
+      preconCreator_ = creator->create(initFileStr);
+      assert(preconCreator_);
     }
 
   private:
-    std::string solverName;
-    std::string prefix;
-
-    ISTLSolverCreator<Matrix, VectorX, VectorB> solverCreator;
+    ISTLSolverCreator<ISTLSolver> solverCreator_;
+    std::shared_ptr<ISTLPrecon> preconCreator_;
 
-    std::shared_ptr<PreconRunner> precon;
-    std::shared_ptr<LinOperator>  linOperator;
-    std::shared_ptr<BaseSolver>   solver;
+    std::shared_ptr<BasePrecon>   precon_;
+    std::shared_ptr<LinOperator>  linOperator_;
+    std::shared_ptr<BaseSolver>   solver_;
   };
 
 } // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/ISTL_Preconditioner.hpp b/src/amdis/linear_algebra/istl/ISTL_Preconditioner.hpp
index 750b562b2e714b7a5d9720962017a65decbe23be..4a653fa796949ba28b3e3511da3628f4aed5cd1f 100644
--- a/src/amdis/linear_algebra/istl/ISTL_Preconditioner.hpp
+++ b/src/amdis/linear_algebra/istl/ISTL_Preconditioner.hpp
@@ -2,71 +2,116 @@
 
 #include <dune/istl/preconditioners.hh>
 
+#include <amdis/CreatorInterface.hpp>
+
 namespace AMDiS
 {
+  template <class Matrix, class VectorX, class VectorB>
+  struct ISTLPreconInterface
+  {
+    using PreconBase = Dune::Preconditioner<VectorX, VectorB>;
+    virtual std::unique_ptr<PreconBase> create(Matrix const& matrix) const = 0;
+  };
+
+  template <class> struct Type {};
 
   // creators for preconditioners, depending on matrix-type
   // ---------------------------------------------------------------------------
 
-  template <class Matrix, class VectorX, class VectorB = VectorX>
-  class ISTLPreconCreator;
-
-  template <class... MatrixRows, class VectorX, class VectorB>
-  class ISTLPreconCreator<Dune::MultiTypeBlockMatrix<MatrixRows...>, VectorX, VectorB>
+  template <class Precon, class Matrix>
+  struct ISTLPrecon
+      : public ISTLPreconInterface<Matrix, typename Precon::domain_type, typename Precon::range_type>
   {
-    using BasePreconditioner = Dune::Preconditioner<VectorX, VectorB>;
+    using VectorX = typename Precon::domain_type;
+    using VectorB = typename Precon::range_type;
+    using Super = ISTLPreconInterface<Matrix, VectorX, VectorB>;
+    using Self = ISTLPrecon;
 
-  public:
-    ISTLPreconCreator(std::string prefix)
-      : prefix(prefix)
-    {}
+    struct Creator : CreatorInterfaceName<Super>
+    {
+      virtual std::shared_ptr<Super> create(std::string prefix) override
+      {
+        return std::make_shared<Self>(prefix);
+      }
+    };
 
-    template <class Matrix>
-    std::unique_ptr<BasePreconditioner> create(std::string /*name*/, Matrix const& A)
+    ISTLPrecon(std::string prefix)
     {
-      double w = 1.0;
-      Parameters::get(prefix + "->relaxation", w);
+      Parameters::get(prefix + "->relaxation", w_);
+    }
 
-      warning("Only Richardson preconditioner available for MultiTypeBlockMatrix!");
-      return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w);
+    using PreconBase = Dune::Preconditioner<VectorX, VectorB>;
+    virtual std::unique_ptr<PreconBase> create(Matrix const& A) const override
+    {
+      return createImpl(A, Type<Precon>{});
     }
 
   private:
-    std::string prefix;
+    template <class P>
+    std::unique_ptr<P> createImpl(Matrix const& A, Type<P>) const
+    {
+      return std::make_unique<P>(A, 1, w_);
+    }
+
+    std::unique_ptr<Dune::SeqILU0<Matrix, VectorX, VectorB>>
+    createImpl(Matrix const& A, Type<Dune::SeqILU0<Matrix, VectorX, VectorB>>) const
+    {
+      return std::make_unique<Dune::SeqILU0<Matrix, VectorX, VectorB>>(A, w_);
+    }
+
+    std::unique_ptr<Dune::Richardson<VectorX, VectorB>>
+    createImpl(Matrix const& /*A*/, Type<Dune::Richardson<VectorX, VectorB>>) const
+    {
+      return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w_);
+    }
+
+    double w_ = 1.0;
   };
 
-  template <class Block, class Alloc, class VectorX, class VectorB>
-  class ISTLPreconCreator<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>
+
+  /// Adds default creators for linear solvers based on `Dune::BCRSMatrix`.
+  /**
+   * Adds creators for full-matrix aware solvers.
+   * - *cg*: conjugate gradient method, \see Dune::CGSolver
+   * - *bcgs*: stabilized bi-conjugate gradient method, \see Dune::BiCGSTABSolver
+   * - *minres*: Minimal residul method, \see Dune::MINRESSolver
+   * - *gmres*: Generalized minimal residula method, \see Dune::RestartedGMResSolver
+   * - *umfpack*: external UMFPACK solver, \see Dune::UMFPack
+   **/
+  template <class Matrix, class VectorX, class VectorB>
+  class DefaultCreators<ISTLPreconInterface<Matrix, VectorX, VectorB>>
   {
-    using BasePreconditioner = Dune::Preconditioner<VectorX, VectorB>;
+    template <class Precon>
+    using PreconCreator
+      = typename ISTLPrecon<Precon, Matrix>::Creator;
 
-  public:
-    ISTLPreconCreator(std::string prefix)
-      : prefix(prefix)
-    {}
+    using Map = CreatorMap<ISTLPreconInterface<Matrix, VectorX, VectorB>>;
 
-    template <class Matrix>
-    std::unique_ptr<BasePreconditioner> create(std::string name, Matrix const& A)
+  public:
+    static void init()
     {
-      double w = 1.0;
-      Parameters::get(prefix + "->relaxation", w);
-
-      if (name == "diag" || name == "jacobi")
-        return std::make_unique<Dune::SeqJac<Matrix, VectorX, VectorB>>(A, 1, w);
-      else if (name == "gs" || name == "gauss_seidel")
-        return std::make_unique<Dune::SeqGS<Matrix, VectorX, VectorB>>(A, 1, w);
-      else if (name == "sor")
-        return std::make_unique<Dune::SeqSOR<Matrix, VectorX, VectorB>>(A, 1, w);
-      else if (name == "ssor")
-        return std::make_unique<Dune::SeqSSOR<Matrix, VectorX, VectorB>>(A, 1, w);
-      else if (name == "ilu" || name == "ilu0")
-        return std::make_unique<Dune::SeqILU0<Matrix, VectorX, VectorB>>(A, w);
-      else
-        return std::make_unique<Dune::Richardson<VectorX, VectorB>>(w);
-    }
+      auto jacobi = new PreconCreator<Dune::SeqJac<Matrix, VectorX, VectorB>>;
+      Map::addCreator("diag", jacobi);
+      Map::addCreator("jacobi", jacobi);
 
-  private:
-    std::string prefix;
+      auto gs = new PreconCreator<Dune::SeqGS<Matrix, VectorX, VectorB>>;
+      Map::addCreator("gs", gs);
+      Map::addCreator("gauss_seidel", gs);
+
+      auto sor = new PreconCreator<Dune::SeqSOR<Matrix, VectorX, VectorB>>;
+      Map::addCreator("sor", sor);
+
+      auto ssor = new PreconCreator<Dune::SeqSSOR<Matrix, VectorX, VectorB>>;
+      Map::addCreator("ssor", ssor);
+
+      auto ilu = new PreconCreator<Dune::SeqILU0<Matrix, VectorX, VectorB>>;
+      Map::addCreator("ilu", ilu);
+      Map::addCreator("ilu0", ilu);
+
+      auto richardson = new PreconCreator<Dune::Richardson<VectorX, VectorB>>;
+      Map::addCreator("richardson", richardson);
+      Map::addCreator("default", richardson);
+    }
   };
 
 } // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/ISTL_Solver.hpp b/src/amdis/linear_algebra/istl/ISTL_Solver.hpp
index 558677dbdf152e31750adcbc4efefe9604381d15..e32d52fe3c726a57a951e40ad99819cc59d542a6 100644
--- a/src/amdis/linear_algebra/istl/ISTL_Solver.hpp
+++ b/src/amdis/linear_algebra/istl/ISTL_Solver.hpp
@@ -1,106 +1,196 @@
 #pragma once
 
+#include <memory>
+
 #include <dune/istl/solvers.hh>
 #include <dune/istl/umfpack.hh>
+#include <dune/istl/superlu.hh>
+
+#include <amdis/CreatorMap.hpp>
+#include <amdis/Initfile.hpp>
+
+#include <amdis/linear_algebra/istl/Fwd.hpp>
+#include <amdis/linear_algebra/istl/DirectRunner.hpp>
+#include <amdis/linear_algebra/istl/ISTLRunner.hpp>
 
 namespace AMDiS
 {
-  // creators for solvers, depending on matrix-type
-  // ---------------------------------------------------------------------------
+  template <class ISTLSolver, bool specialized>
+  struct ISTLSolverCreator
+  {
+    ISTLSolverCreator(std::string prefix)
+    {
+      Parameters::get(prefix + "->info", info_);
+      Parameters::get(prefix + "->max iteration", maxIter_);
+      Parameters::get(prefix + "->relative tolerance", rTol_);
+    }
 
-  template <class Matrix, class VectorX, class VectorB = VectorX>
-  class ISTLSolverCreator;
+    template <class LinOperator, class Precon>
+    std::unique_ptr<ISTLSolver> create(LinOperator& A, Precon& P) const
+    {
+      return std::make_unique<ISTLSolver>(A, P, rTol_, maxIter_, info_);
+    }
+
+    int info_ = 2;
+    std::size_t maxIter_ = 500;
+    double rTol_ = 1.e-6;
+  };
 
-  template <class... MatrixRows, class VectorX, class VectorB>
-  class ISTLSolverCreator<Dune::MultiTypeBlockMatrix<MatrixRows...>, VectorX, VectorB>
+  template <class VectorX>
+  struct ISTLSolverCreator<Dune::RestartedGMResSolver<VectorX>, true>
+      : public ISTLSolverCreator<Dune::RestartedGMResSolver<VectorX>, false>
   {
-    using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
+    using ISTLSolver = Dune::RestartedGMResSolver<VectorX>;
+    using Super = ISTLSolverCreator<ISTLSolver, false>;
 
-  public:
     ISTLSolverCreator(std::string prefix)
-      : prefix(prefix)
-      , info(2)
+      : Super(prefix)
     {
-      Parameters::get(prefix + "->info", info);
+      Parameters::get(prefix + "->restart", restart_);
     }
 
-    template <class Matrix, class Precon>
-    std::unique_ptr<BaseSolver> create(std::string name, Matrix& A, Precon& P)
+    template <class LinOperator, class Precon>
+    std::unique_ptr<ISTLSolver> create(LinOperator& A, Precon& P) const
     {
-      std::size_t maxIter = 500;
-      Parameters::get(prefix + "->max iteration", maxIter);
-
-      double rTol = 1.e-6;
-      Parameters::get(prefix + "->reduction", rTol);
-
-      if (name == "cg")
-        return std::make_unique<Dune::CGSolver<VectorX>>(A, P, rTol, maxIter, info);
-      else if (name == "bicgstab")
-        return std::make_unique<Dune::BiCGSTABSolver<VectorX>>(A, P, rTol, maxIter, info);
-      else if (name == "minres")
-        return std::make_unique<Dune::MINRESSolver<VectorX>>(A, P, rTol, maxIter, info);
-      else if (name == "gmres")
-      {
-        int restart = 30;
-        Parameters::get(prefix + "->restart", restart);
-        return std::make_unique<Dune::RestartedGMResSolver<VectorX>>(A, P, rTol, restart, maxIter, info);
-      }
-      else
-        error_exit("Unknown solver name!");
+      return std::make_unique<ISTLSolver>(A, P, this->rTol_, restart_, this->maxIter_, this->info_);
     }
 
-  private:
-    std::string prefix;
-    int         info;
+    int restart_ = 30;
   };
 
-
-  template <class Block, class Alloc, class VectorX, class VectorB>
-  class ISTLSolverCreator<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>
+  template <class VectorX>
+  struct ISTLSolverCreator<Dune::GeneralizedPCGSolver<VectorX>, true>
+      : public ISTLSolverCreator<Dune::GeneralizedPCGSolver<VectorX>, false>
   {
-    using BaseSolver = Dune::InverseOperator<VectorX, VectorB>;
-    using Matrix     = Dune::BCRSMatrix<Block,Alloc>;
+    using ISTLSolver = Dune::GeneralizedPCGSolver<VectorX>;
+    using Super = ISTLSolverCreator<ISTLSolver, false>;
 
-  public:
     ISTLSolverCreator(std::string prefix)
-      : prefix(prefix)
-      , info(2)
+      : Super(prefix)
     {
-      Parameters::get(prefix + "->info", info);
+      Parameters::get(prefix + "->restart", restart_);
     }
 
     template <class LinOperator, class Precon>
-    std::unique_ptr<BaseSolver> create(std::string name, LinOperator& A, Precon& P)
+    std::unique_ptr<ISTLSolver> create(LinOperator& A, Precon& P) const
     {
-      std::size_t maxIter = 500;
-      Parameters::get(prefix + "->max iteration", maxIter);
-
-      double rTol = 1.e-6;
-      Parameters::get(prefix + "->relative tolerance", rTol);
-
-      if (name == "cg")
-        return std::make_unique<Dune::CGSolver<VectorX>>(A, P, rTol, maxIter, info);
-      else if (name == "bicgstab")
-        return std::make_unique<Dune::BiCGSTABSolver<VectorX>>(A, P, rTol, maxIter, info);
-      else if (name == "minres")
-        return std::make_unique<Dune::MINRESSolver<VectorX>>(A, P, rTol, maxIter, info);
-      else if (name == "gmres")
-      {
-        int restart = 30;
-        Parameters::get(prefix + "->restart", restart);
-        return std::make_unique<Dune::RestartedGMResSolver<VectorX>>(A, P, rTol, restart, maxIter, info);
-      }
+      return std::make_unique<ISTLSolver>(A, P, this->rTol_, this->maxIter_, this->info_, restart_);
+    }
+
+    int restart_ = 30;
+  };
+
+
 #if HAVE_SUITESPARSE_UMFPACK
-      else if (name == "umfpack" || name == "direct")
-        return std::make_unique<Dune::UMFPack<Matrix>>(A, info);
+  template <class Matrix>
+  struct ISTLSolverCreator<Dune::UMFPack<Matrix>, true>
+  {
+    using Solver = Dune::UMFPack<Matrix>;
+    ISTLSolverCreator(std::string prefix)
+    {
+      Parameters::get(prefix + "->info", verbose_);
+    }
+
+    std::unique_ptr<Solver> create(Matrix const& A) const
+    {
+      return std::make_unique<Solver>(A, verbose_);
+    }
+
+  private:
+    int verbose_ = 2;
+  };
 #endif
-      else
-        error_exit("Unknown solver name!");
+
+#if HAVE_SUPERLU
+  template <class Matrix>
+  struct ISTLSolverCreator<Dune::SuperLU<Matrix>, true>
+  {
+    using Solver = Dune::SuperLU<Matrix>;
+    ISTLSolverCreator(std::string prefix)
+    {
+      int info = 2;
+      Parameters::get(prefix + "->info", info);
+      verbose_ = (info != 0);
+
+      Parameters::get(prefix + "->reuse vector", reuseVector_);
+    }
+
+    std::unique_ptr<Solver> create(Matrix const& A) const
+    {
+      return std::make_unique<Solver>(A, verbose_, reuseVector_);
     }
 
   private:
-    std::string prefix;
-    int         info;
+    bool verbose_ = true;
+    bool reuseVector_ = true;
+  };
+#endif
+
+  /// Adds default creators for linear solvers based on `Dune::BCRSMatrix`.
+  /**
+   * Adds creators for full-matrix aware solvers.
+   * - *cg*: conjugate gradient method, \see Dune::CGSolver
+   * - *bcgs*: stabilized bi-conjugate gradient method, \see Dune::BiCGSTABSolver
+   * - *minres*: Minimal residul method, \see Dune::MINRESSolver
+   * - *gmres*: Generalized minimal residual method, \see Dune::RestartedGMResSolver
+   * - *fcg*: Generalized preconditioned conjugate gradient solver, \see Dune::GeneralizedPCGSolver
+   * - *umfpack*: external UMFPACK solver, \see Dune::UMFPack
+   * - *superlu*: external SuperLU solver, \see Dune::SuperLU
+   **/
+  template <class Block, class Alloc, class VectorX, class VectorB>
+  class DefaultCreators<LinearSolverInterface<Dune::BCRSMatrix<Block,Alloc>, VectorX, VectorB>>
+  {
+    using Matrix = Dune::BCRSMatrix<Block,Alloc>;
+    using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;
+
+    template <class ISTLSolver>
+    using SolverCreator
+      = typename LinearSolver<Matrix, VectorX, VectorB,
+          ISTLRunner<ISTLSolver, Matrix, VectorX, VectorB>>::Creator;
+
+    template <class ISTLSolver>
+    using DirectSolverCreator
+      = typename LinearSolver<Matrix, VectorX, VectorB,
+          DirectRunner<ISTLSolver, Matrix, VectorX, VectorB>>::Creator;
+
+    using Map = CreatorMap<SolverBase>;
+
+  public:
+    static void init()
+    {
+      auto cg = new SolverCreator<Dune::CGSolver<VectorX>>;
+      Map::addCreator("cg", cg);
+
+      auto bicgstab = new SolverCreator<Dune::BiCGSTABSolver<VectorX>>;
+      Map::addCreator("bicgstab", bicgstab);
+      Map::addCreator("bcgs", bicgstab);
+
+      auto minres = new SolverCreator<Dune::MINRESSolver<VectorX>>;
+      Map::addCreator("minres", minres);
+
+      auto gmres = new SolverCreator<Dune::RestartedGMResSolver<VectorX>>;
+      Map::addCreator("gmres", gmres);
+
+      // Generalized preconditioned conjugate gradient solver.
+      auto fcg = new SolverCreator<Dune::GeneralizedPCGSolver<VectorX>>;
+      Map::addCreator("fcg", fcg);
+
+#if HAVE_SUITESPARSE_UMFPACK
+      auto umfpack = new DirectSolverCreator<Dune::UMFPack<Matrix>>;
+      Map::addCreator("umfpack", umfpack);
+#endif
+
+#if HAVE_SUPERLU
+      auto superlu = new DirectSolverCreator<Dune::SuperLU<Matrix>>;
+      Map::addCreator("superlu", superlu);
+#endif
+
+#if HAVE_SUITESPARSE_UMFPACK
+      Map::addCreator("direct", umfpack);
+#elif HAVE_SUPERLU
+      Map::addCreator("direct", superlu);
+#endif
+    }
   };
 
 } // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/LinearSolver.hpp b/src/amdis/linear_algebra/istl/LinearSolver.hpp
deleted file mode 100644
index 7799eaa29ae2c1d1407bf049d3f2d6d7f87043c8..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/LinearSolver.hpp
+++ /dev/null
@@ -1,70 +0,0 @@
-#pragma once
-
-#include <amdis/linear_algebra/LinearSolverInterface.hpp>
-#include <amdis/linear_algebra/istl/ISTLRunner.hpp>
-
-namespace AMDiS
-{
-  // A general solver class for istl solvers
-  // ---------------------------------------------------------------------------
-
-  template <class Matrix, class VectorX, class VectorB, class Runner>
-  class LinearSolver
-    : public LinearSolverInterface<Matrix, VectorX, VectorB>
-  {
-    using Super = LinearSolverInterface<Matrix, Vector, Vector>;
-    using RunnerBase = typename Super::RunnerBase;
-
-  public:
-    /// Constructor
-    LinearSolver(std::string name, std::string prefix)
-      : runner(std::make_shared<Runner>(name, prefix))
-    {}
-
-
-  private:
-    /// Implements \ref LinearSolverInterface::solveImpl()
-    virtual void solveImpl(Matrix const& A, VectorX& x, VectorB const& b,
-                           SolverInfo& solverInfo) override
-    {
-      if (solverInfo.doCreateMatrixData()) {
-        // init matrix or wrap block-matrix or ...
-        runner->init(A);
-      }
-
-      // solve the linear system
-      int error = runner->solve(A, x, b, solverInfo);
-      solverInfo.setError(error);
-
-      if (!solverInfo.doStoreMatrixData())
-        runner->exit();
-    }
-
-  private:
-    /// redirect the implementation to a runner. Implements a \ref RunnerInterface
-    std::shared_ptr<Runner> runner;
-  };
-
-
-  template <class Matrix, class VectorX, class VectorB>
-  class SolverCreator<Matrix, VectorX, VectorB>
-  {
-    using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;
-
-    template <class Runner>
-    using Solver = LinearSolver<Matrix, VectorX, VectorB, Runner>;
-
-    /// Instantiate a new linear solver
-    static std::unique_ptr<SolverBase> create(std::string name, std::string prefix)
-    {
-      using ISTL = ISTLRunner<Matrix, VectorX, VectorB>;
-
-      if (name == "cg" || name == "bicgstab" || name == "minres" || name == "gmres")
-        return std::make_unique<Solver<ISTL>>(name, prefix);
-      else
-        error_exit("Unknown solver name!"); // TODO: add direct solver interface here
-    }
-  };
-
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/Preconditioner.hpp b/src/amdis/linear_algebra/istl/Preconditioner.hpp
deleted file mode 100644
index 6c331d53e789198dee9dd064479ef48aca5dd97c..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/Preconditioner.hpp
+++ /dev/null
@@ -1,80 +0,0 @@
-#pragma once
-
-#include <type_traits>
-
-#include <dune/istl/preconditioner.hh>
-#include <dune/istl/solvercategory.hh>
-
-#include <amdis/linear_algebra/istl/ISTL_Preconditioner.hpp>
-
-namespace AMDiS {
-
-namespace detail
-{
-    template <class T, class = typename T::category >
-    std::integral_constant<bool, true> HasCategoryImpl(int);
-
-    template <class T>
-    std::integral_constant<bool, false> HasCategoryImpl(...);
-
-    template <class T>
-    using HasCategory = decltype(HasCategoryImpl<T>(int{}));
-
-} // end namespace detail
-
-
-template <class Matrix, class VectorX, class VectorB>
-class ISTLPreconditioner
-  : public Dune::Preconditioner<VectorX, VectorB>
-  , public PreconditionerInterface<Matrix, VectorX, VectorB>
-{
-public:
-    enum {category=Dune::SolverCategory::sequential};
-
-    ISTLPreconditioner(std::string preconName, std::string prefix)
-      : preconName(preconName)
-      , prefix(prefix)
-      , preconCreator(prefix)
-    {}
-
-    /// Implementation of \ref Dune::Preconditioner::pre()
-    virtual void pre(VectorX& x, VectorB& b) override { runner.pre(x, b); }
-
-    /// Implementation of \ref Dune::Preconditioner::apply()
-    virtual void apply(VectorX& v, const VectorB& d) override { runner.apply(v, d); }
-
-    /// Implementation of \ref Dune::Preconditioner::post()
-    virtual void post(VectorX& x) override { runner.post(x); }
-
-    // -------------------------------------------------------------------------
-
-    /// Implementation of \ref PreconditionerInterface::init()
-    virtual void init(Matrix const& A) override
-    {
-      Matrix& _A = const_cast<Matrix&>(A); // Unschoen!!!
-
-      runner = preconCreator.create(preconName, _A);
-    };
-
-    /// Implementation of \ref PreconditionerInterface::exit()
-    virtual void exit() override
-    {
-      runner.reset();
-    };
-
-    /// Implementation of \ref PreconditionerInterface::solve()
-    virtual void solve(VectorX const& x, VectorX& y) const override
-    {
-      runner.apply(v, d);
-    }
-
-private:
-    std::string preconName;
-    std::string prefix;
-
-    ISTLPreconCreator<Matrix, VectorX, VectorB> preconCreator;
-
-    std::shared_ptr<PreconRunner> runner;
-};
-
-} // end namespace AMDiS
\ No newline at end of file
diff --git a/src/amdis/linear_algebra/istl/SystemMatrix.cpp b/src/amdis/linear_algebra/istl/SystemMatrix.cpp
deleted file mode 100644
index 26d4d975b2d79c524632b40f846bbb0c2b84c0cc..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/SystemMatrix.cpp
+++ /dev/null
@@ -1,11 +0,0 @@
-#define AMDIS_NO_EXTERN_SYSTEMMATRIX
-#include <amdis/linear_algebra/istl/SystemMatrix.hpp>
-#undef AMDIS_NO_EXTERN_SYSTEMMATRIX
-
-
-namespace AMDiS
-{
-  // explicit template instatiation
-  template class SystemMatrix<typename TestTraits<2,1>::FeSpaces>;
-  template class SystemMatrix<typename TestTraits<2,1,2>::FeSpaces>;
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/SystemMatrix.hpp b/src/amdis/linear_algebra/istl/SystemMatrix.hpp
deleted file mode 100644
index fa3e50b219658eef03b2ea288132132eba7c9f9c..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/SystemMatrix.hpp
+++ /dev/null
@@ -1,318 +0,0 @@
-#pragma once
-
-#include <vector>
-#include <string>
-#include <memory>
-#include <tuple>
-
-#include <dune/istl/multitypeblockmatrix.hh>
-#include <dune/istl/multitypeblockvector.hh>
-
-#include <amdis/Output.hpp>
-#include <amdis/common/TupleUtility.hpp>
-#include <amdis/common/Utility.hpp>
-#include <amdis/linear_algebra/istl/DOFMatrix.hpp>
-
-// for explicit instantiation
-#include <amdis/ProblemStatTraits.hpp>
-
-namespace AMDiS
-{
-  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::get<R>(std::forward<MultiMatrix>(multiMatrix))
-          )...};
-      }
-
-      template <std::size_t R, std::size_t... C, class RowFeSpace, class ColFeSpaces, class MultiVector>
-      static DOFMatrixRow<R> make_row(index_t<R>,
-                                      Indices<C...>,
-                                      RowFeSpace&& rowFeSpace,
-                                      ColFeSpaces&& colFeSpace,
-                                      MultiVector&& multiVector)
-      {
-        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) + "]",
-          std::get<C>(std::forward<MultiVector>(multiVector))
-          )...};
-      }
-    };
-
-    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
-
-
-
-
-  /// \brief Container that repesents multiple data-Vectors of different value types.
-  /**
-    *  Represents a \ref Dune::MultiTypeBlockVector + a tuple of corresponding
-    *  feSpaces. This I'th \ref Dune::BlockVector compbined with the I'th FeSpace
-    *  builds a \ref DOFVector  and can be return by the \ref operator[].
-    *
-    *  By default, the ValueTypes are all \ref Dune::FieldMatrix of 1x1 double entry.
-    **/
-  template <class FeSpaces,
-            class FieldType  = double,
-            class Dimensions = Repeat_t<std::tuple_size<FeSpaces>::value, index_t<1>> >
-  class SystemMatrix
-  {
-    template <class RowFeSpace, class RowDim>
-    struct MultiMatrixRowGenerator
-    {
-      template <class ColDim>
-      using ValueTypeGenerator = Dune::FieldMatrix<FieldType, RowDim::value, ColDim::value>;
-      using ValueTypes = MakeTuple_t<ValueTypeGenerator, Dimensions>;
-
-      template <class ColFeSpace, class VT>
-      using DOFMatrixGenerator = DOFMatrix<RowFeSpace, ColFeSpace, VT>;
-      using DOFMatrices  = MakeTuple2_t<DOFMatrixGenerator, FeSpaces, ValueTypes>;
-
-      template <class VT>
-      using BaseMatrixGenerator = Dune::BCRSMatrix<VT, std::allocator<VT>>;
-      using BaseMatrices = MakeTuple_t<BaseMatrixGenerator, ValueTypes>;
-
-      using MultiVector = ExpandTuple_t<Dune::MultiTypeBlockVector, BaseMatrices>;
-    };
-
-    template <class RowFeSpace, class RowDim>
-    using MultiVectorGenerator = typename MultiMatrixRowGenerator<RowFeSpace, RowDim>::MultiVector;
-
-    template <class RowFeSpace, class RowDim>
-    using ValueTypeGenerator = typename MultiMatrixRowGenerator<RowFeSpace, RowDim>::ValueTypes;
-
-    template <class RowFeSpace, class RowDim>
-    using DOFMatrixGenerator = typename MultiMatrixRowGenerator<RowFeSpace, RowDim>::DOFMatrices;
-
-  public:
-    using MultiVectors = MakeTuple2_t<MultiVectorGenerator, FeSpaces, Dimensions>;
-    using ValueTypes   = MakeTuple2_t<ValueTypeGenerator, FeSpaces, Dimensions>;
-    using DOFMatrices  = MakeTuple2_t<DOFMatrixGenerator, FeSpaces, Dimensions>;
-    using MultiMatrix  = ExpandTuple_t<Dune::MultiTypeBlockMatrix, MultiVectors>;
-
-    static constexpr std::size_t nComponents = std::tuple_size<FeSpaces>::value;
-
-    AMDIS_STATIC_ASSERT( nComponents > 0 );
-    AMDIS_STATIC_ASSERT( nComponents == std::tuple_size<Dimensions>::value );
-
-    /// Constructor that constructs new base-vectors
-    SystemMatrix(FeSpaces const& feSpaces)
-      : feSpaces(feSpaces)
-      , matrix{}
-      , dofmatrices(Impl::buildDOFMatrices<DOFMatrices>(feSpaces, feSpaces, matrix))
-    {}
-
-
-    /// Return the number of vector entries
-    static constexpr std::size_t N()
-    {
-      return std::tuple_size<FeSpaces>::value;
-    }
-
-    static constexpr std::size_t M()
-    {
-      return std::tuple_size<FeSpaces>::value;
-    }
-
-    /// Return a shared pointer to the i'th underlaying base vector
-    template <std::size_t R, std::size_t C>
-    auto& getMatrix(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>(matrix));
-    }
-
-    /// Return a shared pointer to the i'th underlaying base vector
-    template <std::size_t R, std::size_t C>
-    auto const& getMatrix(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>(matrix));
-    }
-
-    /// Return the underlying multi vector
-    MultiMatrix& getMatrix()
-    {
-        return matrix;
-    }
-
-    MultiMatrix const& getMatrix() const
-    {
-      return matrix;
-    }
-
-    /// Return the I'th 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);
-    }
-
-    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);
-    }
-
-    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));
-    }
-
-    /// Create a DOFVector with references to the unerlaying data
-    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:
-    FeSpaces const& feSpaces;	///< a tuple of feSpaces
-    MultiMatrix 	matrix;		///< a tuple of data-vectors
-    DOFMatrices     dofmatrices;
-  };
-
-
-  // specialization for 1 component
-  // -------------------------------------------------------------------------
-
-  template <class FeSpace,
-            class FieldType,
-            int dim>
-  class SystemMatrix<std::tuple<FeSpace>, FieldType, std::tuple<index_t<dim>>>
-  {
-    using ValueType     = Dune::FieldMatrix<FieldType, dim, dim>;
-    using DOFMatrixType = DOFMatrix<FeSpace, FeSpace, ValueType>;
-
-  public:
-    using MultiMatrix   = typename DOFMatrixType::BaseMatrix;
-
-    static constexpr std::size_t nComponent = 1;
-
-
-    /// Constructor that constructs new base-vectors
-    template <class FeSpaces>
-    SystemMatrix(FeSpaces const& feSpaces)
-      : dofmatrix(std::get<0>(feSpaces), std::get<0>(feSpaces), "")
-    {}
-
-
-    /// Return the number of vector entries
-    static constexpr std::size_t N()
-    {
-      return 1;
-    }
-
-    static constexpr std::size_t M()
-    {
-      return 1;
-    }
-
-    /// Return a shared pointer to the i'th underlaying base vector
-    template <std::size_t R = 0, std::size_t C = 0>
-    MultiMatrix& getMatrix(const index_t<R> = {}, const index_t<C> = {})
-    {
-      static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
-      return dofmatrix.getMatrix();
-    }
-
-    /// Return a shared pointer to the i'th underlaying base vector
-    template <std::size_t R = 0, std::size_t C = 0>
-    MultiMatrix const& getMatrix(const index_t<R> = {}, const index_t<C> = {}) const
-    {
-      static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
-      return dofmatrix.getMatrix();
-    }
-
-    /// Return the underlying multi vector
-    MultiMatrix const& getMatrix() const
-    {
-      return dofmatrix.getMatrix();
-    }
-
-    MultiMatrix& getMatrix()
-    {
-      return dofmatrix.getMatrix();
-    }
-
-    /// Return the I'th finite element space
-    template <std::size_t I = 0>
-    FeSpace const& getRowFeSpace(const index_t<I> = {}) const
-    {
-      static_assert(I == 0, "Index out of range [0,1)");
-      return dofmatrix.getRowFeSpace();
-    }
-
-    template <std::size_t I = 0>
-    FeSpace const& getColFeSpace(const index_t<I> = {}) const
-    {
-      static_assert(I == 0, "Index out of range [0,1)");
-      return dofmatrix.getColFeSpace();
-    }
-
-    template <std::size_t R = 0, std::size_t C = 0>
-    DOFMatrixType& operator()(const index_t<R> = {}, const index_t<C> = {})
-    {
-      static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
-      return dofmatrix;
-    }
-
-    /// Create a DOFVector with references to the unerlaying data
-    template <std::size_t R = 0, std::size_t C = 0>
-    DOFMatrixType const& operator()(const index_t<R> = {}, const index_t<C> = {}) const
-    {
-      static_assert(R == 0 && C == 0, "Indices out of range [0,1)x[0,1)");
-      return dofmatrix;
-    }
-
-  private:
-    DOFMatrixType  dofmatrix;
-  };
-
-#ifndef AMDIS_NO_EXTERN_SYSTEMMATRIX
-    // explicit instantiation in SystemMatrix.cpp
-    extern template class SystemMatrix<typename TestTraits<2,1>::FeSpaces>;
-    extern template class SystemMatrix<typename TestTraits<2,1,2>::FeSpaces>;
-#endif
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/SystemVector.cpp b/src/amdis/linear_algebra/istl/SystemVector.cpp
deleted file mode 100644
index fc04f99bd3495ed19e973913814e6a57c97732c0..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/SystemVector.cpp
+++ /dev/null
@@ -1,11 +0,0 @@
-#define AMDIS_NO_EXTERN_SYSTEMVECTOR
-#include <amdis/linear_algebra/istl/SystemVector.hpp>
-#undef AMDIS_NO_EXTERN_SYSTEMVECTOR
-
-
-namespace AMDiS
-{
-  // explicit template instatiation
-  template class SystemVector<typename TestTraits<2,1>::FeSpaces>;
-  template class SystemVector<typename TestTraits<2,1,2>::FeSpaces>;
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/istl/SystemVector.hpp b/src/amdis/linear_algebra/istl/SystemVector.hpp
deleted file mode 100644
index 7dfcbb2bad000c96ee0a7fd6ee8e5daf012a08e4..0000000000000000000000000000000000000000
--- a/src/amdis/linear_algebra/istl/SystemVector.hpp
+++ /dev/null
@@ -1,383 +0,0 @@
-#pragma once
-
-#include <vector>
-#include <string>
-#include <memory>
-#include <tuple>
-
-#include <dune/istl/multitypeblockvector.hh>
-
-#include <amdis/Output.hpp>
-#include <amdis/common/Loops.hpp>
-#include <amdis/common/TupleUtility.hpp>
-#include <amdis/common/Utility.hpp>
-
-#include <amdis/linear_algebra/istl/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],
-          std::get<I>(std::forward<MultiVector>(multiVector))
-          )...};
-      }
-    };
-
-    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 \ref Dune::MultiTypeBlockVector + a tuple of corresponding
-    *  feSpaces. This I'th \ref Dune::BlockVector compbined with the I'th FeSpace
-    *  builds a \ref DOFVector  and can be return by the \ref operator[].
-    *
-    *  By default, the ValueTypes are all \ref Dune::FieldVector of 1 double entry.
-    **/
-  template <class FeSpaces,
-            class ValueTypes = Repeat_t<std::tuple_size<FeSpaces>::value, Dune::FieldVector<double,1>> >
-  class SystemVector
-  {
-    using Self = SystemVector;
-
-  public:
-    template <class T> using BlockVectorGenerator = Dune::BlockVector<T, std::allocator<T>>;
-
-    using BaseVectors = MakeTuple_t<BlockVectorGenerator, ValueTypes>;
-    using MultiVector = ExpandTuple_t<Dune::MultiTypeBlockVector, BaseVectors>;
-
-    using DOFVectors  = MakeTuple2_t<DOFVector, FeSpaces, ValueTypes>;
-
-    static constexpr std::size_t nComponents = std::tuple_size<FeSpaces>::value;
-
-    AMDIS_STATIC_ASSERT( nComponents > 0 );
-    AMDIS_STATIC_ASSERT( nComponents == std::tuple_size<ValueTypes>::value );
-
-    /// 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 vector entries
-    static constexpr std::size_t size()
-    {
-      return std::tuple_size<FeSpaces>::value;
-    }
-
-    /// Return number of elements
-    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 std::get<I>(vector);
-    }
-
-    /// 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 std::get<I>(vector);
-    }
-
-    /// Return the underlying multi vector
-    MultiVector& getVector()
-    {
-      return vector;
-    }
-
-    MultiVector const& getVector() 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];
-    }
-
-    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);
-    }
-
-    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);
-    }
-
-    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);
-    }
-
-    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();
-      });
-    }
-
-    void copy(Self const& that)
-    {
-      forEach(range_<0, size()>, [this, &that](const auto I) {
-        std::get<I>(dofvectors).copy(that[I]);
-      });
-    }
-
-    template <class Scalar>
-    std::enable_if_t< std::is_arithmetic<Scalar>::value, Self&>
-    operator*=(Scalar s)
-    {
-      forEach(range_<0, size()>, [this, s](const auto I) {
-        std::get<I>(vector) *= s;
-      });
-      return *this;
-    }
-
-    template <class Scalar>
-    std::enable_if_t< std::is_arithmetic<Scalar>::value, Self&>
-    operator=(Scalar s)
-    {
-      forEach(range_<0, size()>, [this, s](const auto I) {
-        std::get<I>(vector) = s;
-      });
-      return *this;
-    }
-
-  private:
-    FeSpaces const& feSpaces;	///< a tuple of feSpaces
-    std::vector<std::string> names;
-
-    MultiVector vector;		///< a tuple of data-vectors
-    DOFVectors	dofvectors;	///< a tuple of dofvectors referencing \ref vector
-  };
-
-
-  // specialization for the case of only 1 component.
-
-  template <class FeSpace, class ValueType>
-  class SystemVector<std::tuple<FeSpace>, std::tuple<ValueType>>
-  {
-    using Self = SystemVector;
-
-  public:
-    using MultiVector = Dune::BlockVector<ValueType>;
-    using DOFVectors  = DOFVector<FeSpace, ValueType>;
-
-    static constexpr std::size_t nComponent = 1;
-
-    /// Constructor that constructs new base-vectors
-    template <class FeSpaces>
-    SystemVector(FeSpaces const& feSpaces, std::vector<std::string> const& names)
-      : dofvector(std::get<0>(feSpaces), names.front())
-    {
-      compress();
-    }
-
-
-    /// Return the number of vector entries
-    static constexpr std::size_t size()
-    {
-      return 1;
-    }
-
-    /// Return number of elements
-    int count() const
-    {
-      return 1;
-    }
-
-    /// Return a shared pointer to the i'th underlaying base vector
-    template <std::size_t I = 0>
-    MultiVector& getVector(const index_t<I> = {})
-    {
-      static_assert(I == 0, "Index out of range [0,1)");
-      return dofvector.getVector();
-    }
-
-    /// Return a shared pointer to the i'th underlaying base vector
-    template <std::size_t I = 0>
-    MultiVector const& getVector(const index_t<I> = {}) const
-    {
-      static_assert(I == 0, "Index out of range [0,1)");
-      return dofvector.getVector();
-    }
-
-    /// Return the underlying multi vector
-    MultiVector& getVector()
-    {
-      return dofvector.getVector();
-    }
-
-    MultiVector const& getVector() const
-    {
-      return dofvector.getVector();
-    }
-
-    /// Return the I'th finite element space
-    template <std::size_t I = 0>
-    auto const& getFeSpace(const index_t<I> = {}) const
-    {
-      static_assert(I == 0, "Index out of range [0,1)");
-      return dofvector.getFeSpace();
-    }
-
-    /// Return the name of the i'th DOFVector
-    std::string getName(std::size_t i = 0) const
-    {
-      test_exit(i == 0, "Index out of range [0,1)");
-      return dofvector.getName();
-    }
-
-    template <std::size_t I>
-    auto& operator[](const index_t<I>)
-    {
-      static_assert(I == 0, "Index out of range [0,1)");
-      return dofvector;
-    }
-
-    template <std::size_t I>
-    auto const& operator[](const index_t<I>) const
-    {
-      static_assert(I == 0, "Index out of range [0,1)");
-      return dofvector;
-    }
-
-    template <std::size_t I>
-    auto& getDOFVector(const index_t<I> = {})
-    {
-      static_assert(I == 0, "Index out of range [0,1)");
-      return dofvector;
-    }
-
-    template <std::size_t I>
-    auto const& getDOFVector(const index_t<I> = {}) const
-    {
-      static_assert(I == 0, "Index out of range [0,1)");
-      return dofvector;
-    }
-
-    /// 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> = {})
-    {
-      static_assert(I == 0, "Index out of range [0,1)");
-      dofvector.compress();
-    }
-
-    /// Resize the \ref vectors to the sizes of the \ref feSpaces.
-    void compress()
-    {
-      dofvector.compress();
-    }
-
-    void copy(Self const& that)
-    {
-      dofvector.copy(that[_0]);
-    }
-
-    template <class Scalar>
-    std::enable_if_t< std::is_arithmetic<Scalar>::value, Self&>
-    operator*=(Scalar s)
-    {
-      dofvector.getVector() *= s;
-      return *this;
-    }
-
-    template <class Scalar>
-    std::enable_if_t< std::is_arithmetic<Scalar>::value, Self&>
-    operator=(Scalar s)
-    {
-      dofvector.getVector() = s;
-      return *this;
-    }
-
-  private:
-    DOFVectors      dofvector;  ///< a tuple of sofvectors referencing \ref vector
-  };
-
-
-
-
-#ifndef AMDIS_NO_EXTERN_SYSTEMVECTOR
-    // explicit instantiation in SystemVector.cpp
-    extern template class SystemVector<typename TestTraits<2,1>::FeSpaces>;
-    extern template class SystemVector<typename TestTraits<2,1,2>::FeSpaces>;
-#endif
-
-} // end namespace AMDiS
diff --git a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
index 7adbf1b879864e0688f4a48912797915a3e7c76e..67c702b32c0ed743f16fda44e65c5ffaaf59a35c 100644
--- a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
@@ -12,21 +12,14 @@
 
 #include <amdis/Output.hpp>
 #include <amdis/common/ClonablePtr.hpp>
+#include <amdis/linear_algebra/DOFMatrixBase.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,
-            class ColBasisType,
+  template <class RowBasisType, class ColBasisType,
             class ValueType = double>
   class DOFMatrix
   {
@@ -41,7 +34,7 @@ namespace AMDiS
     using BaseMatrix = mtl::compressed2D<ValueType>;
 
     /// The index/size - type
-    using size_type  = typename BaseMatrix::size_type;
+    using size_type = typename BaseMatrix::size_type;
 
     /// The type of the elements of the DOFMatrix
     using value_type = ValueType;
@@ -51,17 +44,14 @@ namespace AMDiS
 
   public:
     /// Constructor. Constructs new BaseMatrix.
-    DOFMatrix(RowBasis const& rowBasis,
-              ColBasis const& colBasis)
+    DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis)
       : rowBasis_(&rowBasis)
       , colBasis_(&colBasis)
       , matrix_(ClonablePtr<BaseMatrix>::make())
     {}
 
-    /// Constructor. Takes pointer of data-matrix.
-    DOFMatrix(RowBasis const& rowBasis,
-              ColBasis const& colBasis,
-              BaseMatrix& matrix)
+    /// Constructor. Takes a reference to a base matrix
+    DOFMatrix(RowBasis const& rowBasis, ColBasis const& colBasis, BaseMatrix& matrix)
       : rowBasis_(&rowBasis)
       , colBasis_(&colBasis)
       , matrix_(matrix)
diff --git a/src/amdis/linear_algebra/mtl/DOFVector.hpp b/src/amdis/linear_algebra/mtl/DOFVector.hpp
index f5ff128ce2b19f3f5163b89d22837854cf5fdf65..7287b807d33fe77f7770839cad744ffe95c22f48 100644
--- a/src/amdis/linear_algebra/mtl/DOFVector.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFVector.hpp
@@ -1,59 +1,50 @@
 #pragma once
 
-#include <algorithm>
-#include <memory>
-#include <string>
-
 #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/DOFVectorInterface.hpp>
+#include <amdis/linear_algebra/DOFVectorBase.hpp>
 #include <amdis/utility/MultiIndex.hpp>
 
 namespace AMDiS
 {
   /// The basic container that stores a base vector and a corresponding basis
-  template <class GlobalBasis, class ValueType = double>
+  template <class BasisType, class ValueType = double>
   class DOFVector
-      : public DOFVectorInterface
+      : public DOFVectorBase<BasisType>
   {
     using Self = DOFVector;
+    using Super = DOFVectorBase<BasisType>;
 
   public:
     /// The type of the \ref basis
-    using Basis = GlobalBasis;
+    using Basis = BasisType;
 
     /// The type of the base vector
     using BaseVector = mtl::vec::dense_vector<ValueType>;
 
     /// The index/size - type
-    using size_type  = typename GlobalBasis::size_type;
+    using size_type  = typename BasisType::size_type;
 
     /// The type of the elements of the DOFVector
     using value_type = ValueType;
 
     /// Constructor. Constructs new BaseVector.
-    DOFVector(GlobalBasis const& basis)
-      : basis_(&basis)
+    DOFVector(BasisType const& basis)
+      : Super(basis)
       , vector_(ClonablePtr<BaseVector>::make())
     {
       compress();
     }
 
     /// Constructor. Takes reference to existing BaseVector
-    DOFVector(GlobalBasis const& basis, BaseVector& vector)
-      : basis_(&basis)
+    DOFVector(BasisType const& basis, BaseVector& vector)
+      : Super(basis)
       , vector_(vector)
     {}
 
-    /// Return the basis \ref basis_ of the vector
-    Basis const& basis() const
-    {
-        return *basis_;
-    }
-
     /// Return the data-vector \ref vector_
     BaseVector const& vector() const
     {
@@ -66,17 +57,11 @@ namespace AMDiS
       return *vector_;
     }
 
-    /// Return the size of the \ref basis
-    size_type size() const
-    {
-      return basis_->dimension();
-    }
-
     /// Resize the \ref vector to the size of the \ref basis.
     virtual void compress() override
     {
-      if (num_rows(*vector_) != size()) {
-        resize(size());
+      if (mtl::vec::size(*vector_) != this->size()) {
+        resize(this->size());
         *vector_ = value_type(0);
       }
     }
@@ -93,8 +78,8 @@ namespace AMDiS
     value_type const& operator[](Index idx) const
     {
       size_type i = flatMultiIndex(idx);
-      test_exit_dbg(i < num_rows(*vector_),
-        "Index ", i, " out of range [0,", num_rows(*vector_), ")" );
+      test_exit_dbg(i < mtl::vec::size(*vector_),
+        "Index ", i, " out of range [0,", mtl::vec::size(*vector_), ")" );
       return (*vector_)[i];
     }
 
@@ -103,8 +88,8 @@ namespace AMDiS
     value_type& operator[](Index idx)
     {
       size_type i = flatMultiIndex(idx);
-      test_exit_dbg(i < num_rows(*vector_),
-        "Index ", i, " out of range [0,", num_rows(*vector_), ")" );
+      test_exit_dbg(i < mtl::vec::size(*vector_),
+        "Index ", i, " out of range [0,", mtl::vec::size(*vector_), ")" );
       return (*vector_)[i];
     }
 
@@ -133,36 +118,22 @@ namespace AMDiS
     }
 
   private:
-    /// The finite element space / basis associated with the data vector
-    Basis const* basis_;
-
     /// The data-vector (can hold a new BaseVector or a pointer to external data
     ClonablePtr<BaseVector> vector_;
   };
 
-#if DUNE_HAVE_CXX_CLASS_TEMPLATE_ARGUMENT_DEDUCTION
-  // Deduction rules
-  template <class GlobalBasis>
-  DOFVector(GlobalBasis const& basis)
-    -> DOFVector<GlobalBasis, double>;
-
-  template <class GlobalBasis, class ValueType>
-  DOFVector(GlobalBasis const& basis, mtl::vec::dense_vector<ValueType>& coefficients)
-    -> DOFVector<GlobalBasis, ValueType>;
-#endif
-
   /// Constructor a dofvector from given basis and name
-  template <class GlobalBasis, class ValueType = double>
-  DOFVector<GlobalBasis, ValueType>
-  makeDOFVector(GlobalBasis const& basis)
+  template <class Basis, class ValueType = double>
+  DOFVector<Basis, ValueType>
+  makeDOFVector(Basis const& basis)
   {
     return {basis};
   }
 
   /// Constructor a dofvector from given basis, name, and coefficient vector
-  template <class GlobalBasis, class ValueType>
-  DOFVector<GlobalBasis, ValueType>
-  makeDOFVector(GlobalBasis const& basis, mtl::vec::dense_vector<ValueType>& coefficients)
+  template <class Basis, class ValueType>
+  DOFVector<Basis, ValueType>
+  makeDOFVector(Basis const& basis, mtl::vec::dense_vector<ValueType>& coefficients)
   {
     return {basis, coefficients};
   }
diff --git a/src/amdis/linear_algebra/mtl/ITL_Solver.hpp b/src/amdis/linear_algebra/mtl/ITL_Solver.hpp
index 27a76f5c460dd2cf34af7f4db01359e52b9b8d5b..258e37500039cb67f5c59f617abb1f111c465aef 100644
--- a/src/amdis/linear_algebra/mtl/ITL_Solver.hpp
+++ b/src/amdis/linear_algebra/mtl/ITL_Solver.hpp
@@ -16,7 +16,7 @@
 
 #include <amdis/CreatorMap.hpp>
 #include <amdis/Initfile.hpp>
-#include <amdis/linear_algebra/mtl/LinearSolver.hpp>
+#include <amdis/linear_algebra/LinearSolver.hpp>
 #include <amdis/linear_algebra/mtl/KrylovRunner.hpp>
 #include <amdis/linear_algebra/mtl/UmfpackRunner.hpp>
 
@@ -393,12 +393,12 @@ namespace AMDiS
 
     template <class ITLSolver>
     using SolverCreator
-      = typename LinearSolver<Matrix, Vector,
+      = typename LinearSolver<Matrix, Vector, Vector,
           KrylovRunner<ITLSolver, Matrix, Vector>>::Creator;
 
 #ifdef HAVE_UMFPACK
     using UmfpackSolverCreator
-      = typename LinearSolver<Matrix, Vector,
+      = typename LinearSolver<Matrix, Vector, Vector,
           UmfpackRunner<Matrix, Vector>>::Creator;
 #endif
 
diff --git a/src/amdis/linear_algebra/mtl/Preconditioner.hpp b/src/amdis/linear_algebra/mtl/Preconditioner.hpp
index 6d9922ea9d3f945b270c603a7a5caca08fa9ba04..cc610a82cacf5a498ca9c67bcf55de638e4296bc 100644
--- a/src/amdis/linear_algebra/mtl/Preconditioner.hpp
+++ b/src/amdis/linear_algebra/mtl/Preconditioner.hpp
@@ -10,7 +10,8 @@ namespace AMDiS
    * \brief Wrapper for using ITL preconditioners in AMDiS.
    */
   template <class Matrix, class Vector, class PreconRunner>
-  class Preconditioner : public PreconditionerInterface<Matrix, Vector>
+  class Preconditioner
+      : public PreconditionerInterface<Matrix, Vector>
   {
     using Self = Preconditioner;
     using Super = PreconditionerInterface<Matrix, Vector>;