diff --git a/dune/matrix-vector/test/triangularsolvetest.cc b/dune/matrix-vector/test/triangularsolvetest.cc index fec0232568d2999b5b7b38dc1f249be14ae7be89..2e9f24cc57c34f5c7a2d3d2d5c84825328b2c94d 100644 --- a/dune/matrix-vector/test/triangularsolvetest.cc +++ b/dune/matrix-vector/test/triangularsolvetest.cc @@ -2,8 +2,6 @@ #include "config.h" #endif -#undef NDEBUG - #include <dune/common/bitsetvector.hh> #include <dune/common/fmatrix.hh> #include <dune/common/fvector.hh> @@ -21,17 +19,25 @@ double const reg = 1e-6; template <int n, bool lower> bool test() { + using ft = double; + using LocalMatrix = Dune::FieldMatrix<ft, 1, 1>; + using Matrix = Dune::BCRSMatrix<LocalMatrix>; + using LocalVector = Dune::FieldVector<ft, 1>; + using Vector = Dune::BlockVector<LocalVector>; + bool passed = true; + std::random_device randomDevice; std::mt19937 generator(randomDevice()); - Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>> M; + Matrix M; Dune::MatrixIndexSet indices(n, n); 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); + std::uniform_int_distribution<int> integerDistribution(1, + finalJ - initialJ); auto aroundTen = [&]() { auto x = integerDistribution(generator); return x <= 10; @@ -43,7 +49,7 @@ bool test() { } indices.exportIdx(M); - std::uniform_real_distribution<> unconstrainedDistribution(-1, 1); + std::uniform_real_distribution<ft> 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) { @@ -55,17 +61,17 @@ bool test() { } } - Dune::BlockVector<Dune::FieldVector<double, 1>> sol(n); + Vector sol(n); for (auto &x : sol) x = unconstrainedDistribution(generator); - Dune::BlockVector<Dune::FieldVector<double, 1>> b(n); + Vector b(n); M.mv(sol, b); Dune::BitSetVector<1> *ignore = nullptr; { auto x = lower ? Dune::MatrixVector::lowerTriangularSolve(M, b, ignore) : Dune::MatrixVector::upperTriangularSolve(M, b, ignore); - Dune::BlockVector<Dune::FieldVector<double, 1>> y(n); + Vector y(n); M.mv(x, y); auto diff = diffDune(b, y); std::cout << "|x - A^{-1}b| = " << diff << std::endl; @@ -75,19 +81,19 @@ bool test() { auto x = lower ? Dune::MatrixVector::upperTriangularSolve(M, b, ignore, true) : Dune::MatrixVector::lowerTriangularSolve(M, b, ignore, true); - Dune::BlockVector<Dune::FieldVector<double, 1>> y(n); + Vector y(n); M.mtv(x, y); auto diff = diffDune(b, y); - std::cout << "|x - A^{-T}b| = " << diff << " (manual transpose)" << std::endl; + std::cout << "|x - A^{-T}b| = " << diff << " (manual transpose)" + << std::endl; passed &= diff < tol; } { - Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>> Mt; + Matrix 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); + auto x = lower ? Dune::MatrixVector::upperTriangularSolve(Mt, b, ignore) + : Dune::MatrixVector::lowerTriangularSolve(Mt, b, ignore); + Vector y(n); M.mtv(x, y); auto diff = diffDune(b, y); std::cout << "|x - A^{-T}b| = " << diff << std::endl;