From dd48435bea3cea2c5f891035249d032b9ca325b0 Mon Sep 17 00:00:00 2001
From: Christian Engwer <christi@dune-project.org>
Date: Sat, 5 Dec 2015 23:42:30 +0100
Subject: [PATCH] [test] improve the tests for multiple rhs
---
dune/istl/solvers.hh | 28 +++++++++++------------
dune/istl/test/CMakeLists.txt | 37 ++++++++++++++++--------------
dune/istl/test/multirhstest.cc | 42 +++++++++++++++++++---------------
3 files changed, 58 insertions(+), 49 deletions(-)
diff --git a/dune/istl/solvers.hh b/dune/istl/solvers.hh
index ea7182256..5e9f207d8 100644
--- a/dune/istl/solvers.hh
+++ b/dune/istl/solvers.hh
@@ -20,7 +20,7 @@
#include "preconditioner.hh"
#include <dune/common/deprecated.hh>
#include <dune/common/exceptions.hh>
-#include <dune/common/simd.hh>
+#include <dune/common/conditional.hh>
#include <dune/common/rangeutilities.hh>
#include <dune/common/timer.hh>
#include <dune/common/ftraits.hh>
@@ -519,7 +519,7 @@ namespace Dune {
this->printOutput(std::cout,i,def);
_prec.post(x); // postprocess preconditioner
- res.iterations = i; // fill statistics
+ res.iterations = i; // fill statistics
res.reduction = static_cast<double>(max_value(max_value(def/def0)));
res.conv_rate = pow(res.reduction,1.0/i);
res.elapsed = watch.elapsed();
@@ -1093,7 +1093,7 @@ namespace Dune {
using std::abs;
using std::max;
using std::min;
- real_type eps = 1e-15;
+ const real_type eps = 1e-15;
real_type norm_dx = abs(dx);
real_type norm_dy = abs(dy);
real_type norm_max = max(norm_dx, norm_dy);
@@ -1102,22 +1102,22 @@ namespace Dune {
// we rewrite the code in a vectorizable fashion
cs = cond(norm_dy < eps,
field_type(1.0),
- cond(norm_dx < 1e-15,
+ cond(norm_dx < eps,
field_type(0.0),
cond(norm_dy > norm_dx,
- 1.0/sqrt(1.0 + temp*temp)*temp,
+ field_type(1.0)/sqrt(1.0 + temp*temp)*temp,
// dy is real in exact arithmetic
// so we don't need to conjugate here
- 1.0/sqrt(1.0 + temp*temp)*dx/norm_dx*dy/norm_dy)));
+ field_type(1.0)/sqrt(1.0 + temp*temp)*dx/norm_dx*dy/norm_dy)));
sn = cond(norm_dy < eps,
field_type(0.0),
cond(norm_dx < eps,
field_type(1.0),
cond(norm_dy > norm_dx,
- 1.0/sqrt(1.0 + temp*temp),
+ field_type(1.0)/sqrt(1.0 + temp*temp),
// dy and dx is real in exact arithmetic
// so we don't have to conjugate both of them
- 1.0/sqrt(1.0 + temp*temp)*dy/dx)));
+ field_type(1.0)/sqrt(1.0 + temp*temp)*dy/dx)));
}
SeqScalarProduct<X> ssp;
@@ -1446,7 +1446,7 @@ namespace Dune {
using std::abs;
using std::max;
using std::min;
- real_type eps = 1e-15;
+ const real_type eps = 1e-15;
real_type norm_dx = abs(dx);
real_type norm_dy = abs(dy);
real_type norm_max = max(norm_dx, norm_dy);
@@ -1455,18 +1455,18 @@ namespace Dune {
// we rewrite the code in a vectorizable fashion
cs = cond(norm_dy < eps,
field_type(1.0),
- cond(norm_dx < 1e-15,
+ cond(norm_dx < eps,
field_type(0.0),
cond(norm_dy > norm_dx,
- 1.0/sqrt(1.0 + temp*temp)*temp,
- 1.0/sqrt(1.0 + temp*temp)*dx/norm_dx*conjugate(dy)/norm_dy)));
+ field_type(1.0)/sqrt(1.0 + temp*temp)*temp,
+ field_type(1.0)/sqrt(1.0 + temp*temp)*dx/norm_dx*conjugate(dy)/norm_dy)));
sn = cond(norm_dy < eps,
field_type(0.0),
cond(norm_dx < eps,
field_type(1.0),
cond(norm_dy > norm_dx,
- 1.0/sqrt(1.0 + temp*temp),
- 1.0/sqrt(1.0 + temp*temp)*conjugate(dy/dx))));
+ field_type(1.0)/sqrt(1.0 + temp*temp),
+ field_type(1.0)/sqrt(1.0 + temp*temp)*conjugate(dy/dx))));
}
diff --git a/dune/istl/test/CMakeLists.txt b/dune/istl/test/CMakeLists.txt
index 665724947..a78ee0b83 100644
--- a/dune/istl/test/CMakeLists.txt
+++ b/dune/istl/test/CMakeLists.txt
@@ -20,7 +20,10 @@ dune_add_test(SOURCES matrixutilstest.cc)
dune_add_test(SOURCES matrixtest.cc)
+if(Vc_FOUND)
dune_add_test(SOURCES multirhstest.cc)
+ add_dune_vc_flags(multirhstest)
+endif(Vc_FOUND)
dune_add_test(SOURCES bvectortest.cc)
@@ -52,35 +55,35 @@ dune_add_test(SOURCES solvertest.cc)
dune_add_test(SOURCES solveraborttest.cc)
# Pardiso tests
-dune_add_test(SOURCES test_pardiso.cc)
+ dune_add_test(SOURCES test_pardiso.cc)
# SuperLU tests
-dune_add_test(NAME superlustest
- SOURCES superlutest.cc
- COMPILE_DEFINITIONS SUPERLU_NTYPE=0)
+ dune_add_test(NAME superlustest
+ SOURCES superlutest.cc
+ COMPILE_DEFINITIONS SUPERLU_NTYPE=0)
-dune_add_test(SOURCES superlutest.cc)
+ dune_add_test(SOURCES superlutest.cc)
-dune_add_test(NAME superluctest
- SOURCES superlutest.cc
- COMPILE_DEFINITIONS SUPERLU_NTYPE=2)
+ dune_add_test(NAME superluctest
+ SOURCES superlutest.cc
+ COMPILE_DEFINITIONS SUPERLU_NTYPE=2)
-dune_add_test(NAME superluztest
- SOURCES superlutest.cc
- COMPILE_DEFINITIONS SUPERLU_NTYPE=3)
+ dune_add_test(NAME superluztest
+ SOURCES superlutest.cc
+ COMPILE_DEFINITIONS SUPERLU_NTYPE=3)
-dune_add_test(SOURCES complexrhstest.cc
- COMPILE_DEFINITIONS SUPERLU_NTYPE=3)
+ dune_add_test(SOURCES complexrhstest.cc
+ COMPILE_DEFINITIONS SUPERLU_NTYPE=3)
# SuiteSparse tests
-dune_add_test(SOURCES umfpacktest.cc)
-set_directory_properties(PROPERTIES ADDITIONAL_MAKE_CLEAN_FILES "umfpack_decomp")
+ dune_add_test(SOURCES umfpacktest.cc)
+ set_directory_properties(PROPERTIES ADDITIONAL_MAKE_CLEAN_FILES "umfpack_decomp")
dune_add_test(SOURCES ldltest.cc)
dune_add_test(SOURCES spqrtest.cc)
-dune_add_test(SOURCES overlappingschwarztest.cc)
+ dune_add_test(SOURCES overlappingschwarztest.cc)
# MPI tests
dune_add_test(SOURCES matrixredisttest.cc
@@ -89,6 +92,6 @@ dune_add_test(SOURCES matrixredisttest.cc
dune_add_test(SOURCES vectorcommtest.cc
CMAKE_GUARD MPI_FOUND)
-dune_add_test(SOURCES matrixmarkettest.cc)
+ dune_add_test(SOURCES matrixmarkettest.cc)
exclude_from_headercheck(complexdata.hh)
diff --git a/dune/istl/test/multirhstest.cc b/dune/istl/test/multirhstest.cc
index b4c9a871f..bc0ecc95f 100644
--- a/dune/istl/test/multirhstest.cc
+++ b/dune/istl/test/multirhstest.cc
@@ -4,6 +4,8 @@
// start with including some headers
#include "config.h"
+// #undef HAVE_VC
+
#include <iostream> // for input/output to shell
#include <fstream> // for input/output to files
#include <vector> // STL vector class
@@ -35,6 +37,7 @@ struct Random {
}
};
+#if HAVE_VC
template<typename T, typename A>
struct Random<Vc::Vector<T,A>> {
static Vc::Vector<T,A> gen()
@@ -50,6 +53,7 @@ struct Random<Vc::SimdArray<T,N,V,M>> {
return Vc::SimdArray<T,N,V,M>::Random();
}
};
+#endif
template <typename V>
V detectVectorType(Dune::LinearOperator<V,V> &);
@@ -61,6 +65,8 @@ void run_test (std::string precName, std::string solverName, Operator & op, Solv
using FT = typename Vector::field_type;
Dune::Timer t;
+ std::cout << "Trying " << solverName << "(" << precName << ")"
+ << " with " << Dune::className<FT>() << std::endl;
for (unsigned int run = 0; run < Runs; run++) {
// set up system
Vector x(N),b(N);
@@ -74,8 +80,7 @@ void run_test (std::string precName, std::string solverName, Operator & op, Solv
Dune::InverseOperatorResult r;
solver.apply(x,b,r);
}
- std::cout << "Test " << Runs << " run(s) " << solverName << "(" << precName << ")"
- << " with " << Dune::className<FT>() << " took " << t.stop() << std::endl;
+ std::cout << Runs << " run(s) took " << t.stop() << std::endl;
}
template<typename Operator, typename Prec>
@@ -84,21 +89,21 @@ void test_all_solvers(std::string precName, Operator & op, Prec & prec, unsigned
using Vector = decltype(detectVectorType(op));
double reduction = 1e-4;
- int verb = 0;
+ int verb = 1;
Dune::LoopSolver<Vector> loop(op,prec,reduction,18000,verb);
Dune::CGSolver<Vector> cg(op,prec,reduction,8000,verb);
Dune::BiCGSTABSolver<Vector> bcgs(op,prec,reduction,8000,verb);
Dune::GradientSolver<Vector> grad(op,prec,reduction,18000,verb);
- // Dune::RestartedGMResSolver<Vector> gmres(op,prec,reduction,40,8000,verb);
- // Dune::MINRESSolver<Vector> minres(op,prec,reduction,8000,verb);
+ Dune::RestartedGMResSolver<Vector> gmres(op,prec,reduction,40,8000,verb);
+ Dune::MINRESSolver<Vector> minres(op,prec,reduction,8000,verb);
Dune::GeneralizedPCGSolver<Vector> gpcg(op,prec,reduction,8000,verb);
- run_test(precName, "Loop", op,loop,N,Runs);
+ // run_test(precName, "Loop", op,loop,N,Runs);
run_test(precName, "CG", op,cg,N,Runs);
- run_test(precName, "Gradient", op,bcgs,N,Runs);
- run_test(precName, "RestartedGMRes", op,grad,N,Runs);
- // run_test(precName, op,gmres,N,Runs);
- // run_test(precName, "MINRes", op,minres,N,Runs);
+ run_test(precName, "BiCGStab", op,bcgs,N,Runs);
+ run_test(precName, "Gradient", op,grad,N,Runs);
+ run_test(precName, "RestartedGMRes", op,gmres,N,Runs);
+ run_test(precName, "MINRes", op,minres,N,Runs);
run_test(precName, "GeneralizedPCG", op,gpcg,N,Runs);
}
@@ -121,12 +126,12 @@ void test_all(unsigned int Runs = 1)
Dune::MatrixAdapter<Matrix,Vector,Vector> op(A); // make linear operator from A
// create all preconditioners
- Dune::SeqJac<Matrix,Vector,Vector> jac(A,1,1); // Jacobi preconditioner
- Dune::SeqGS<Matrix,Vector,Vector> gs(A,1,0.5); // GS preconditioner
- Dune::SeqSOR<Matrix,Vector,Vector> sor(A,1,1.9520932); // SOR preconditioner
- Dune::SeqSSOR<Matrix,Vector,Vector> ssor(A,1,1.0); // SSOR preconditioner
- Dune::SeqILU0<Matrix,Vector,Vector> ilu0(A,1.0); // preconditioner object
- Dune::SeqILUn<Matrix,Vector,Vector> ilu1(A,1,0.92); // preconditioner object
+ Dune::SeqJac<Matrix,Vector,Vector> jac(A,1,0.1); // Jacobi preconditioner
+ Dune::SeqGS<Matrix,Vector,Vector> gs(A,1,0.1); // GS preconditioner
+ Dune::SeqSOR<Matrix,Vector,Vector> sor(A,1,0.1); // SOR preconditioner
+ Dune::SeqSSOR<Matrix,Vector,Vector> ssor(A,1,0.1); // SSOR preconditioner
+ Dune::SeqILU0<Matrix,Vector,Vector> ilu0(A,0.1); // preconditioner object
+ Dune::SeqILUn<Matrix,Vector,Vector> ilu1(A,1,0.1); // preconditioner object
// run the sub-tests
test_all_solvers("Jacobi", op,jac,N,Runs);
@@ -137,12 +142,13 @@ void test_all(unsigned int Runs = 1)
test_all_solvers("ILU1", op,ilu1,N,Runs);
}
-int main ()
+int main (int argc, char ** argv)
{
test_all<float>();
test_all<double>();
- test_all<Vc::double_v>();
#if HAVE_VC
+ test_all<Vc::float_v>();
+ test_all<Vc::double_v>();
test_all<Vc::Vector<double, Vc::VectorAbi::Scalar>>();
test_all<Vc::SimdArray<double,2>>();
test_all<Vc::SimdArray<double,2,Vc::Vector<double, Vc::VectorAbi::Scalar>,1>>();
--
GitLab