Skip to content
Snippets Groups Projects
Commit c8a41d3b authored by Elias Pipping's avatar Elias Pipping
Browse files

Tests: Cleanup

Whitespace, typedefs, remove #undef left over from debugging
parent 4b984142
No related branches found
No related tags found
No related merge requests found
......@@ -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;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment