diff --git a/dune/solvers/iterationsteps/multigridstep.cc b/dune/solvers/iterationsteps/multigridstep.cc
index 0361d6d34b1adadc0e1c7959a77e985cebb0f4c1..785b36df717004969d37cd2cf6f21f3e290a30d5 100644
--- a/dune/solvers/iterationsteps/multigridstep.cc
+++ b/dune/solvers/iterationsteps/multigridstep.cc
@@ -73,17 +73,14 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
 
     xHierarchy_.resize(numLevels());
     rhsHierarchy_.resize(numLevels());
-    matrixHierarchy_.resize(numLevels());
     ignoreNodesHierarchy_.resize(numLevels());
 
-    for (int i=0; i<int(matrixHierarchy_.size())-1; i++)
+    for (int i=0; i<int(xHierarchy_.size())-1; i++)
     {
-        matrixHierarchy_[i].reset();
         xHierarchy_[i].reset();
         ignoreNodesHierarchy_[i] = NULL;
     }
 
-    matrixHierarchy_.back() = this->mat_;
     xHierarchy_.back()      = Dune::stackobject_to_shared_ptr(*(this->x_));
     rhsHierarchy_.back()    = *(this->rhs_);
 
@@ -117,22 +114,35 @@ void MultigridStep<MatrixType, VectorType, BitVectorType>::preprocess()
     // /////////////////////////////////////////////////////
     //   Assemble the complete hierarchy of matrices
     // /////////////////////////////////////////////////////
-    for (int i=this->numLevels()-2; i>=0; i--)
+    if (not matrixHierarchyManuallySet_)
     {
-        this->matrixHierarchy_[i] = std::make_shared<MatrixType>();
-        this->xHierarchy_[i] = std::make_shared<VectorType>();
+        matrixHierarchy_.resize(numLevels());
+        for (int i=0; i<int(matrixHierarchy_.size())-1; i++)
+            matrixHierarchy_[i].reset();
+        matrixHierarchy_.back() = this->mat_;
+        for (int i=this->numLevels()-2; i>=0; i--)
+        {
+            this->matrixHierarchy_[i] = std::make_shared<MatrixType>();
 
-        // Compute which entries are present in the (sparse) coarse grid stiffness
-        // matrices.
-        this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+            // Compute which entries are present in the (sparse) coarse grid stiffness
+            // matrices.
+            this->mgTransfer_[i]->galerkinRestrictSetOccupation(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+
+            // setup matrix
+            this->mgTransfer_[i]->galerkinRestrict(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+        }
+        matrixHierarchyManuallySet_ = false;
+    }
 
-        // setup matrix
-        this->mgTransfer_[i]->galerkinRestrict(*this->matrixHierarchy_[i+1], *std::const_pointer_cast<MatrixType>(this->matrixHierarchy_[i]));
+    for (int i=this->numLevels()-2; i>=0; i--)
+    {
+        this->xHierarchy_[i] = std::make_shared<VectorType>();
 
         // Set solution vector sizes for the lower levels
         MatrixVector::resize(*(this->xHierarchy_[i]),*this->matrixHierarchy_[i]);
     }
 
+
     // /////////////////////////////////////////////////////
     //   Setup dirichlet bitfields
     // /////////////////////////////////////////////////////
diff --git a/dune/solvers/iterationsteps/multigridstep.hh b/dune/solvers/iterationsteps/multigridstep.hh
index 35ae17fce302781687c7f4e854b5c5f19282123d..1a03622be2fbbf15ccf8fae4f0cc029e463d3dce 100644
--- a/dune/solvers/iterationsteps/multigridstep.hh
+++ b/dune/solvers/iterationsteps/multigridstep.hh
@@ -208,6 +208,24 @@ namespace Dune {
             basesolver_ = wrap_own_share<Solver>(std::forward<BaseSolver>(baseSolver));
         }
 
+        /**
+         * \brief Set pre-coarsened matrix hierarchy
+         *
+         * This allows to set the coarsened matrix hierarchy manually.
+         * It is assumed, that the finest matrix coincides with
+         * the one passed to setProblem(). The latter is ignored.
+         * This has to be called before preprocess().
+         */
+        template<class M>
+        void setMatrixHirarchy(const std::vector<std::shared_ptr<M>>& matrixHierarchy)
+        {
+            matrixHierarchy_.resize(matrixHierarchy.size());
+            for(auto i : Dune::range(matrixHierarchy.size()))
+                matrixHierarchy_[i] = matrixHierarchy[i];
+            matrixHierarchyManuallySet_ = true;
+        }
+
+
     protected:
         /** \brief The presmoothers, one for each level */
         std::vector<std::shared_ptr<LinearStepType> > presmoother_;
@@ -226,6 +244,8 @@ namespace Dune {
 
         //! Number of coarse corrections
         int mu_;
+
+        bool matrixHierarchyManuallySet_= false;
     public:
         //! Variable used to store the current level
         int level_;