Tests: Use a sparse matrix

Otherwise, rounding errors dominate once the problem becomes larger

dune/matrix-vector/test/triangularsolvetest.cc

@@ -16,34 +16,47 @@
 #include "common.hh"
 
 double const tol = 1e-10;
+double const reg = 1e-6;
 
 template <int n, bool lower>
 bool test() {
   bool passed = true;
+  std::random_device randomDevice;
+  std::mt19937 generator(randomDevice());
+
   Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>> M;
   Dune::MatrixIndexSet indices(n, n);
-  for (size_t i = 0; i < n; ++i)
-    for (size_t j = lower ? 0 : i; j < (lower ? i + 1 : n); ++j)
-      indices.add(i, j);
+  for (size_t i = 0; i < n; ++i) {
+    size_t initialJ = lower ? 0 : i;
+    size_t finalJ = lower ? i + 1 : n;
+    for (size_t j = initialJ; j < finalJ; ++j) {
+      std::uniform_int_distribution<int> integerDistribution(1, finalJ - initialJ);
+      auto aroundTen = [&]() {
+        auto x = integerDistribution(generator);
+        return x <= 10;
+      };
+      // Add a diagonal entry and approximately ten non-zero entries
+      if (i == j or aroundTen())
+        indices.add(i, j);
+    }
+  }
   indices.exportIdx(M);
 
-  std::random_device randomDevice;
-  std::mt19937 generator(randomDevice());
-  std::uniform_real_distribution<> distribution(0.1, 0.4);
-
+  std::uniform_real_distribution<> lowDistribution(reg, 0.4);
   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 = distribution(generator);
+      *cIt = lowDistribution(generator);
       if (i == j)
         *cIt += 0.6;
     }
   }
 
+  std::uniform_real_distribution<> unconstrainedDistribution(reg, 0.4);
   Dune::BlockVector<Dune::FieldVector<double, 1>> sol(n);
   for (auto &x : sol)
-    x = distribution(generator);
+    x = unconstrainedDistribution(generator);
   Dune::BlockVector<Dune::FieldVector<double, 1>> b(n);
   M.mv(sol, b);
 
@@ -72,7 +85,7 @@ bool test() {
 
 int main() {
   bool passed = true;
-  passed &= test<1000, false>();
-  passed &= test<1000, true>();
+  passed &= test<5000, false>();
+  passed &= test<5000, true>();
   return passed ? 0 : 1;
 }