diff --git a/dune/solvers/test/common.hh b/dune/solvers/test/common.hh
index 7015536debd45974e770bb404e2b8fed4f9c5275..203b1e028d4bce9f7e0413d2b3f1dd45f736aac1 100644
--- a/dune/solvers/test/common.hh
+++ b/dune/solvers/test/common.hh
@@ -26,6 +26,9 @@
 #include <dune/alugrid/grid.hh>
 #endif
 
+#include <dune/matrix-vector/addtodiagonal.hh>
+
+#include <dune/solvers/common/defaultbitvector.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/norms/twonorm.hh>
 
@@ -64,13 +67,6 @@ void constructPQ1Pattern(const GridView& gridView, Matrix& matrix)
     indices.exportIdx(matrix);
 }
 
-template<class T, int n, int m, class F>
-void addToDiagonal(Dune::FieldMatrix<T, n, m>& M, const F& x)
-{
-    for(int i=0; i<n; ++i)
-        M[i][i] += x;
-}
-
 template<class GridView, class Matrix>
 void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
 {
@@ -131,11 +127,11 @@ void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
                     int jGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(j).subEntity(), dim);
 
                     double zij = (gradients[i] * gradients[j]) * z;
-                    addToDiagonal(matrix[iGlobal][jGlobal], zij);
-                    addToDiagonal(matrix[jGlobal][iGlobal], zij);
+                    Dune::MatrixVector::addToDiagonal(matrix[iGlobal][jGlobal], zij);
+                    Dune::MatrixVector::addToDiagonal(matrix[jGlobal][iGlobal], zij);
                 }
                 double zii = (gradients[i] * gradients[i]) * z;
-                addToDiagonal(matrix[iGlobal][iGlobal], zii);
+                Dune::MatrixVector::addToDiagonal(matrix[iGlobal][iGlobal], zii);
             }
         }
     }
@@ -314,15 +310,19 @@ void markBoundaryDOFs(const GridView& gridView, BitVector& isBoundary)
     }
 }
 
+/**
+ *
+ * \tparam blocksize If this is zero, then double will be used instead of FieldVector/FieldMatrix
+ */
 template <size_t blocksize, class GridView, bool trivialDirichletOnly = true,
           bool makePositiveDefinitive = not trivialDirichletOnly>
 class SymmetricSampleProblem {
 public:
-  typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
-  typedef typename Dune::FieldVector<double, blocksize> VectorBlock;
+  using MatrixBlock = std::conditional_t<blocksize==0, double, Dune::FieldMatrix<double, blocksize, blocksize> >;
+  using VectorBlock = std::conditional_t<blocksize==0, double, Dune::FieldVector<double, blocksize> >;
   typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
   typedef typename Dune::BlockVector<VectorBlock> Vector;
-  typedef typename Dune::BitSetVector<blocksize> BitVector;
+  using BitVector = Dune::Solvers::DefaultBitVector_t<Vector>;
   typedef EnergyNorm<Matrix, Vector> Norm;
 
   SymmetricSampleProblem(GridView const& gridView) {
@@ -333,15 +333,21 @@ public:
     energyNorm.setMatrix(&A);
 
     ignore.resize(A.N());
-    ignore.unsetAll();
+    if constexpr (std::is_same_v<decltype(ignore),Dune::BitSetVector<blocksize> >)
+      ignore.unsetAll();
+    else
+      std::fill(ignore.begin(), ignore.end(), false);
     if (trivialDirichletOnly)
         markBoundaryDOFs(gridView, ignore);
     else {
         // Mark the first component only
-        Dune::BitSetVector<blocksize> boundaryNodes(A.N(), false);
+        BitVector boundaryNodes(A.N(), false);
         markBoundaryDOFs(gridView, boundaryNodes);
         for (size_t i=0; i < boundaryNodes.size(); ++i)
-            ignore[i][0] = boundaryNodes[i][0];
+            if constexpr (std::is_same_v<decltype(ignore),Dune::BitSetVector<blocksize> >)
+                ignore[i][0] = boundaryNodes[i][0];
+            else
+                ignore[i] = boundaryNodes[i];
     }
 
     u.resize(A.N());
@@ -357,10 +363,17 @@ public:
     // Set up a random 'solution'
     u = 0;
     for (size_t i = 0; i < u.size(); ++i)
-      for (size_t j = 0; j < blocksize; ++j) {
+      if constexpr (blocksize==0)
+      {
         if (makePositiveDefinitive)
-          A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]);
-        u_ex[i][j] = distribution(generator);
+          A[i][i] += 0.5*std::abs(A[0][0]);
+        u_ex[i] = distribution(generator);
+      } else {
+        for (size_t j = 0; j < blocksize; ++j) {
+          if (makePositiveDefinitive)
+            A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]);
+          u_ex[i][j] = distribution(generator);
+        }
       }
 
     // Construct right hand side corresponding to the 'solution'
@@ -368,9 +381,15 @@ public:
 
     // set dirichlet values
     for (size_t i = 0; i < u.size(); ++i)
-      for (size_t j = 0; j < blocksize; ++j)
-        if (ignore[i][j])
-          u[i][j] = u_ex[i][j];
+      if constexpr (blocksize==0)
+      {
+        if (ignore[i])
+          u[i] = u_ex[i];
+      } else {
+        for (size_t j = 0; j < blocksize; ++j)
+          if (ignore[i][j])
+            u[i][j] = u_ex[i][j];
+      }
   }
 
   Matrix A;