diff --git a/dune/matrix-vector/test/triangularsolvetest.cc b/dune/matrix-vector/test/triangularsolvetest.cc
index bca54c400deb2486e80885bb947e50d84c7f5055..fec0232568d2999b5b7b38dc1f249be14ae7be89 100644
--- a/dune/matrix-vector/test/triangularsolvetest.cc
+++ b/dune/matrix-vector/test/triangularsolvetest.cc
@@ -11,6 +11,7 @@
 #include <dune/istl/bvector.hh>
 #include <dune/istl/matrixindexset.hh>
 
+#include "../transpose.hh"
 #include "../triangularsolve.hh"
 
 #include "common.hh"
@@ -42,18 +43,18 @@ bool test() {
   }
   indices.exportIdx(M);
 
-  std::uniform_real_distribution<> lowDistribution(reg, 0.4);
+  std::uniform_real_distribution<> unconstrainedDistribution(-1, 1);
   for (auto it = M.begin(); it != M.end(); ++it) {
     size_t const i = it.index();
     for (auto cIt = it->begin(); cIt != it->end(); ++cIt) {
       size_t const j = cIt.index();
-      *cIt = lowDistribution(generator);
+      *cIt = unconstrainedDistribution(generator);
+      // Make sure the diagonal is bounded away from zero
       if (i == j)
-        *cIt += 0.6;
+        *cIt += (*cIt > 0 ? 1 : -1) + reg;
     }
   }
 
-  std::uniform_real_distribution<> unconstrainedDistribution(reg, 0.4);
   Dune::BlockVector<Dune::FieldVector<double, 1>> sol(n);
   for (auto &x : sol)
     x = unconstrainedDistribution(generator);
@@ -67,7 +68,7 @@ bool test() {
     Dune::BlockVector<Dune::FieldVector<double, 1>> y(n);
     M.mv(x, y);
     auto diff = diffDune(b, y);
-    std::cout << "Difference: " << diff << std::endl;
+    std::cout << "|x - A^{-1}b| = " << diff << std::endl;
     passed &= diff < tol;
   }
   {
@@ -77,7 +78,19 @@ bool test() {
     Dune::BlockVector<Dune::FieldVector<double, 1>> y(n);
     M.mtv(x, y);
     auto diff = diffDune(b, y);
-    std::cout << "Difference: " << diff << std::endl;
+    std::cout << "|x - A^{-T}b| = " << diff << " (manual transpose)" << std::endl;
+    passed &= diff < tol;
+  }
+  {
+    Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>> Mt;
+    Dune::MatrixVector::transpose(M, Mt);
+    auto x = lower
+                 ? Dune::MatrixVector::upperTriangularSolve(Mt, b, ignore)
+                 : Dune::MatrixVector::lowerTriangularSolve(Mt, b, ignore);
+    Dune::BlockVector<Dune::FieldVector<double, 1>> y(n);
+    M.mtv(x, y);
+    auto diff = diffDune(b, y);
+    std::cout << "|x - A^{-T}b| = " << diff << std::endl;
     passed &= diff < tol;
   }
   return passed;
@@ -85,7 +98,9 @@ bool test() {
 
 int main() {
   bool passed = true;
+  std::cout << "Testing random upper triangular matrix" << std::endl;
   passed &= test<5000, false>();
+  std::cout << "Testing lower triangular matrix" << std::endl;
   passed &= test<5000, true>();
   return passed ? 0 : 1;
 }