diff --git a/CHANGELOG.md b/CHANGELOG.md
index dca526b52608dabdf36911f9fe410b2614a20b15..f5dfb64e0d57cb277ea99aeb2e6be4f435876bb3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,8 @@
 # Master (will become release 2.10)
 
+- `UMFPackSolver` accepts each correctly formed blocked matrix. The user has to make sure that the vector types of `x` and `rhs` are compatible to the matrix.
+  The main advantage is that it is now possible to use `MultiTypeBlockMatrix`.
+
 - A new solver `ProximalNewtonSolver` is added which solves non-smooth minimization problems.
 
 # Release 2.9
diff --git a/dune/solvers/solvers/umfpacksolver.hh b/dune/solvers/solvers/umfpacksolver.hh
index 195d632cdca65ef41374d9e7d8ce018db4693d29..1c0743888a8b33826bcb6249327f958a150c5d6f 100644
--- a/dune/solvers/solvers/umfpacksolver.hh
+++ b/dune/solvers/solvers/umfpacksolver.hh
@@ -61,8 +61,6 @@ public:
 
   void solve() override
   {
-    // We may use the original rhs, but ISTL modifies it, so we need a non-const type here
-    VectorType mutableRhs = *rhs_;
 
     if (not this->hasIgnore())
     {
@@ -72,83 +70,86 @@ public:
       Dune::InverseOperatorResult statistics;
       Dune::UMFPack<MatrixType> solver(*matrix_);
       solver.setOption(UMFPACK_PRL, 0);   // no output at all
+
+      // We may use the original rhs, but ISTL modifies it, so we need a non-const type here
+      VectorType mutableRhs = *rhs_;
       solver.apply(*x_, mutableRhs, statistics);
 
     }
     else
     {
-      ///////////////////////////////////////////////////////////////////////////////////////////
-      //  Extract the set of matrix rows that do not correspond to ignored degrees of freedom.
-      //  Unfortunately, not all cases are handled by the ISTL UMFPack solver.  Currently,
-      //  you can only remove complete block rows.  If a block is only partially ignored,
-      //  the dune-istl UMFPack solver cannot do it, and the code here will throw an exception.
-      //  All this can be fixed, but it needs going into the istl UMFPack code.
-      ///////////////////////////////////////////////////////////////////////////////////////////
-      std::set<std::size_t> nonIgnoreRows;
-      for (size_t i=0; i<matrix_->N(); i++)
-      {
-        auto const &ignore = this->ignore()[i];
-        if constexpr (std::is_convertible<decltype(ignore), bool>::value)
-        {
-          if (!ignore)
-            nonIgnoreRows.insert(i);
-        } else {
-          if (ignore.none())
-            nonIgnoreRows.insert(i);
-          else if (not ignore.all())
-            DUNE_THROW(Dune::NotImplemented, "Individual blocks must be either ignored completely, or not at all");
-        }
-      }
-
       // Construct the solver
       Dune::InverseOperatorResult statistics;
       Dune::UMFPack<MatrixType> solver;
       solver.setOption(UMFPACK_PRL, 0);   // no output at all
       // We eliminate all rows and columns(!) from the matrix that correspond to ignored degrees of freedom.
       // Here is where the sparse LU decomposition is happenening.
-      solver.setSubMatrix(*matrix_,nonIgnoreRows);
+      solver.setMatrix(*matrix_,this->ignore());
+
+
+      // total number of dofs
+      auto N = flatVectorForEach(*rhs_, [](auto&&, auto&&){});
 
       // We need to modify the rhs vector by static condensation.
-      VectorType shortRhs(nonIgnoreRows.size());
-      int shortRowCount=0;
-      for (auto it = nonIgnoreRows.begin(); it!=nonIgnoreRows.end(); ++shortRowCount, ++it)
-        shortRhs[shortRowCount] = mutableRhs[*it];
+      std::vector<bool> flatIgnore(N);
+      size_t numberOfIgnoredDofs = 0;
+
+      flatVectorForEach(this->ignore(), [&](auto&& entry, auto&& offset)
+      {
+        flatIgnore[offset] = entry;
+        if ( entry )
+        {
+          numberOfIgnoredDofs++;
+        }
+      });
+
 
-      shortRowCount = 0;
-      for (size_t i=0; i<matrix_->N(); i++)
+      using field_type = typename MatrixType::field_type;
+      std::vector<field_type> shortRhs(N-numberOfIgnoredDofs);
+
+      // mapping of long indices to short indices
+      std::vector<size_t> subIndices(N,std::numeric_limits<size_t>::max());
+
+      int shortRhsCount=0;
+      flatVectorForEach(*rhs_, [&](auto&& entry, auto&& offset)
       {
-        if constexpr (std::is_convertible<decltype(this->ignore()[i]), const bool>::value) {
-          if (this->ignore()[i])
-            continue;
-        } else {
-          if (this->ignore()[i].all())
-            continue;
+        if ( not flatIgnore[offset] )
+        {
+          shortRhs[shortRhsCount] = entry;
+          subIndices[offset] = shortRhsCount;
+          shortRhsCount++;
         }
+      });
 
-        auto cIt    = (*matrix_)[i].begin();
-        auto cEndIt = (*matrix_)[i].end();
+      std::vector<field_type> flatX(N);
+      flatVectorForEach(*x_, [&](auto&& entry, auto&& offset)
+      {
+        flatX[offset] = entry;
+      });
 
-        for (; cIt!=cEndIt; ++cIt)
-          if constexpr (std::is_convertible<decltype(this->ignore()[cIt.index()]), const bool>::value) {
-            if (this->ignore()[cIt.index()])
-              Dune::Impl::asMatrix(*cIt).mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
-          } else {
-            if (this->ignore()[cIt.index()].all())
-              Dune::Impl::asMatrix(*cIt).mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
-          }
 
-        shortRowCount++;
-      }
+      flatMatrixForEach(*matrix_, [&](auto&& entry, auto&& row, auto&& col)
+      {
+        if ( flatIgnore[col] and not flatIgnore[row] )
+        {
+          shortRhs[ subIndices[ row ] ] -= entry * flatX[ col ];
+        }
+      });
 
       // Solve the reduced system
-      VectorType shortX(nonIgnoreRows.size());
-      solver.apply(shortX, shortRhs, statistics);
+      std::vector<field_type> shortX(N-numberOfIgnoredDofs);
 
-      // Blow up the solution vector
-      shortRowCount=0;
-      for (auto it = nonIgnoreRows.begin(); it!=nonIgnoreRows.end(); ++shortRowCount, ++it)
-        (*x_)[*it] = shortX[shortRowCount];
+      // call the raw-pointer variant of the ISTL UMPACK solver
+      solver.apply(shortX.data(), shortRhs.data());
 
+      // Blow up the solution vector
+      flatVectorForEach(*x_, [&](auto&& entry, auto&& offset)
+      {
+        if ( not flatIgnore[ offset ] )
+        {
+          entry = shortX[ subIndices[ offset ] ];
+        }
+      });
     }
   }