diff --git a/CMakeLists.txt b/CMakeLists.txt
index 942cea2a91bd86f8c6a971c0d7887db4fd3a94b0..f5e749f46116e7760c7acf6d8f09c18e947db40a 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -26,6 +26,7 @@ dune_project()
add_subdirectory(cmake/modules)
add_subdirectory(dune)
add_subdirectory(doc)
+add_subdirectory(bin)
# if Python bindings are enabled, include necessary sub directories.
if( DUNE_ENABLE_PYTHONBINDINGS )
diff --git a/bin/CMakeLists.txt b/bin/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..bdc9a18dff8ee9b75c88f27fa797699dd2991e35
--- /dev/null
+++ b/bin/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_executable(solve_mm EXCLUDE_FROM_ALL solve_mm.cc)
+add_dune_all_flags(solve_mm)
+dune_symlink_to_source_files(FILES solve_mm.ini)
diff --git a/bin/solve_mm.cc b/bin/solve_mm.cc
new file mode 100644
index 0000000000000000000000000000000000000000..c78e818d9050c0db2e3d3c87bb8d04c79e349501
--- /dev/null
+++ b/bin/solve_mm.cc
@@ -0,0 +1,161 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#include <config.h>
+
+// Uncomment the following line if you want to use direct solvers, but have MPI installed on you system
+//#undef HAVE_MPI
+
+#include <fenv.h>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/matrixmarket.hh>
+
+// include solvers and preconditioners for the solverfactory
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/cholmod.hh>
+#include <dune/istl/superlu.hh>
+#include <dune/istl/umfpack.hh>
+#include <dune/istl/ldl.hh>
+#include <dune/istl/spqr.hh>
+#include <dune/istl/paamg/amg.hh>
+
+
+#include "solve_mm.hh"
+
+using VectorFieldType = double;
+
+typedef Dune::BCRSMatrix<Dune::FieldMatrix<Dune::Simd::Scalar<VectorFieldType>,1,1>> Mat;
+typedef Dune::BlockVector<Dune::FieldVector<VectorFieldType,1>> Vec;
+
+void printHelp(){
+ std::cout << "solve_mm" << std::endl
+ << "Loads and solves a system from MatrixMarket" <<std::endl
+ << "format and stores the result in MatrixMarket format." << std::endl
+ << "This program is thought as test environment to play" << std::endl
+ << "around with different solvers and preconditioners." << std::endl
+ << "Parameters are read in a ParameterTree from solve_mm.ini" << std::endl
+ << "but can also be passed in by command line arguments." << std::endl
+ << std::endl
+ << "-matrix\t\tFilename of the MatrixMarket file of the operator matrix" << std::endl
+ << "-rhs\t\tFilename of the MatrixMarket file of the right-hand side" << std::endl
+ << "-verbose\tVerbosity (default: 1)" << std::endl
+ << "-random_rhs\tIf 1 generates a random RHS instead of reading it from a file (default: 0)" << std::endl
+ << "-distributed\tLoads a distributed MatrixMarket format (see matrixmarket.hh) (default: 0)" << std::endl
+ << "-redistribute\tRedistributes the matrix using ParMETIS (default: 0)" << std::endl
+ << "-ini\t\tFilename of the ini-file (default:solve_mm.ini)" << std::endl
+ << "-check_residual\tWhether to compute the defect at the end (default: 1)" << std::endl
+ << "-output\t\tFilename of the output filename in which the result is stored" << std::endl
+ << "-FP_EXCEPT\tenables floating point exceptions (default: 0)" << std::endl
+ << std::endl
+ << "The subtree 'solver' in the ParameterTree is passed to the SolverFactory" << std::endl;
+}
+
+int main(int argc, char** argv){
+ auto& mpihelper = Dune::MPIHelper::instance(argc, argv);
+
+ if(argc > 1 && argv[1][0]=='-' && argv[1][1] == 'h'){
+ printHelp();
+ return 0;
+ }
+
+ Dune::ParameterTreeParser::readOptions(argc, argv, config);
+ Dune::ParameterTreeParser::readINITree(config.get("ini","solve_mm.ini"), config, false);
+
+ if(mpihelper.size() > 1 && !config.get("distributed", false) && !config.get("redistribute", false)){
+ std::cerr << "To run this program in parallel you either need to load a distributed"
+ << " system (-distributed) or redistribute the system using parmetis"
+ << " (-redistribute)." << std::endl;
+ return 1;
+ }
+
+ if(config.get("FP_EXCEPT", false))
+ feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);// | FE_UNDERFLOW);
+
+ int verbose = config.get("verbose", 1);
+ if(mpihelper.rank() > 0)
+ verbose = 0;
+ std::shared_ptr<Vec> rhs = std::make_shared<Vec>();
+ std::shared_ptr<Mat> m = std::make_shared<Mat>();
+
+ if(verbose){
+ std::cout << "Processes: " << mpihelper.size() << std::endl;
+ std::cout << "Loading system... " << std::flush;
+ }
+ Dune::Timer t;
+#if HAVE_MPI
+ auto oocomm = loadSystem(m, rhs, config, mpihelper.getCommunication());
+ typedef Dune::OverlappingSchwarzOperator<Mat, Vec, Vec, OOCOMM> OP;
+ typedef Dune::OverlappingSchwarzScalarProduct<Vec, OOCOMM> SP;
+ std::shared_ptr<OP> op = std::make_shared<OP>(*m, *oocomm);
+ std::shared_ptr<SP> sp = std::make_shared<SP>(*oocomm);
+#else
+ loadSystem(m,rhs,config);
+ typedef Dune::MatrixAdapter<Mat, Vec, Vec> OP;
+ typedef Dune::ScalarProduct<Vec> SP;
+ std::shared_ptr<OP> op = std::make_shared<OP>(*m);
+ std::shared_ptr<SP> sp = std::make_shared<SP>();
+#endif
+ if(verbose)
+ std::cout << t.elapsed() << " s" << std::endl;
+
+ if(config.get("redistribute", false) && mpihelper.size() > 1){
+#if HAVE_PARMETIS && HAVE_MPI
+ Dune::Timer t;
+ if(verbose > 0)
+ std::cout << "Redistributing... " << std::flush;
+ redistribute(m, rhs, oocomm);
+ op = std::make_shared<OP>(*m, *oocomm);
+ sp = std::make_shared<SP>(*oocomm);
+ if(verbose > 0)
+ std::cout << t.elapsed() << " s" << std::endl;
+#else
+ std::cerr << "ParMETIS is necessary to redistribute the matrix." << std::endl;
+#endif
+ }
+
+ if(rhs->size() != m->N()){
+ std::cerr << "Dimensions of matrix and rhs do not match." << std::endl;
+ return 1;
+ }
+
+ Vec x(m->M());
+ typedef typename Vec::field_type field_type;
+ x=field_type(0.0);
+
+ // one output is enough!
+ if(mpihelper.rank()!=0)
+ config["solver.verbose"] = "0";
+
+ std::unique_ptr<Vec> pRHS;
+ bool check_residual = config.get("check_residual", true);
+ if(check_residual)
+ pRHS = std::make_unique<Vec>(*rhs);
+
+ solve(op, *rhs, x, config.sub("solver"),verbose);
+
+ if(check_residual){
+ op->applyscaleadd(-1.0, x, *pRHS);
+ auto defect = sp->norm(*pRHS);
+ if(verbose > 0)
+ std::cout << "Final defect: " << Dune::Simd::io(defect) << std::endl;
+ }
+
+ std::string outputfilename = config.get<std::string>("output", "");
+ if(outputfilename.size() != 0){
+#if HAVE_MPI
+ Dune::storeMatrixMarket(x, outputfilename, *oocomm, true);
+#else
+ std::ostringstream rfilename;
+ rfilename << outputfilename << ".mm";
+ Dune::storeMatrixMarket(x, rfilename.str());
+#endif
+ }
+ return 0;
+}
diff --git a/bin/solve_mm.hh b/bin/solve_mm.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3b6e2806d1877186f4687e266217d0eb8fd8065b
--- /dev/null
+++ b/bin/solve_mm.hh
@@ -0,0 +1,132 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+
+#ifndef BIN_SOLVE_MM_HH
+#define BIN_SOLVE_MM_HH
+
+#include <dune/common/parametertree.hh>
+#include <dune/common/timer.hh>
+#include <dune/istl/solverfactory.hh>
+
+Dune::ParameterTree config;
+
+template<class OP>
+void solve(const std::shared_ptr<OP>& op,
+ typename OP::range_type& rhs,
+ typename OP::domain_type& x,
+ const Dune::ParameterTree& config,
+ int verbose = 1){
+ Dune::Timer t;
+ if(verbose)
+ std::cout << "Initializing solver... " << std::flush;
+ Dune::initSolverFactories<OP>();
+ auto solver = Dune::getSolverFromFactory(op, config);
+ if(verbose){
+ std::cout << t.elapsed() << " s" << std::endl;
+ std::cout << "Solving system..." << std::flush;
+ }
+ t.reset();
+ Dune::InverseOperatorResult res;
+ solver->apply(x,rhs,res);
+ if(verbose)
+ std::cout << t.elapsed() << " s" << std::endl;
+}
+
+template<class Num>
+std::enable_if_t<Dune::IsNumber<Num>::value> fillRandom(Num& n){
+ using namespace Dune::Simd;
+ for(size_t l=0;l<lanes(n);++l){
+ lane(l,n) = (Scalar<Num>(std::rand()+1)/(Scalar<Num>(RAND_MAX)));
+ }
+}
+
+template<class Vec>
+std::enable_if_t<!Dune::IsNumber<Vec>::value> fillRandom(Vec& v){
+ for(auto& n : v){
+ fillRandom(n);
+ }
+}
+
+#if HAVE_MPI
+typedef Dune::OwnerOverlapCopyCommunication<long long> OOCOMM;
+template<class Mat, class Vec, class Comm>
+std::shared_ptr<OOCOMM> loadSystem(std::shared_ptr<Mat>& m,
+ std::shared_ptr<Vec>& rhs,
+ const Dune::ParameterTree& config,
+ Comm comm){
+ std::string matrixfilename = config.get<std::string>("matrix");
+ std::string rhsfilename;
+ if(!config.get("random_rhs", false))
+ rhsfilename = config.get<std::string>("rhs");
+ bool distributed = config.get("distributed", false);
+ std::shared_ptr<OOCOMM> oocomm;
+ if(distributed){
+ oocomm = std::make_shared<OOCOMM>(MPI_COMM_WORLD);
+ loadMatrixMarket(*m, matrixfilename, *oocomm);
+ if(config.get("random_rhs", false)){
+ rhs->resize(m->N());
+ srand(42);
+ fillRandom(*rhs);
+ }else{
+ loadMatrixMarket(*rhs, rhsfilename, *oocomm, false);
+ }
+ }else{
+ oocomm = std::make_shared<OOCOMM>(comm);
+ if(comm.rank()==0){
+ loadMatrixMarket(*m, matrixfilename);
+ if(config.get("random_rhs", false)){
+ rhs->resize(m->N());
+ fillRandom(*rhs);
+ }else{
+ loadMatrixMarket(*rhs, rhsfilename);
+ }
+ }
+ }
+ oocomm->remoteIndices().template rebuild<false>();
+ return oocomm;
+}
+#else
+template<class Mat, class Vec>
+void loadSystem(std::shared_ptr<Mat>& m,
+ std::shared_ptr<Vec>& rhs,
+ const Dune::ParameterTree& config){
+ std::string matrixfilename = config.get<std::string>("matrix");
+ std::string rhsfilename;
+ if(!config.get("random_rhs", false))
+ rhsfilename = config.get<std::string>("rhs");
+ loadMatrixMarket(*m, matrixfilename);
+ if(config.get("random_rhs", false)){
+ rhs->resize(m->N());
+ fillRandom(*rhs);
+ }else
+ loadMatrixMarket(*rhs, rhsfilename);
+}
+#endif
+
+#if HAVE_MPI
+template<class Mat, class Vec>
+void redistribute(std::shared_ptr<Mat>& m,
+ std::shared_ptr<Vec>& rhs,
+ std::shared_ptr<OOCOMM>& oocomm){
+ typedef typename Dune::Amg::MatrixGraph<Mat> MatrixGraph;
+ Dune::RedistributeInformation<OOCOMM> ri;
+ std::shared_ptr<Mat> m_redist = std::make_shared<Mat>();
+ std::shared_ptr<OOCOMM> oocomm_redist = std::make_shared<OOCOMM>(MPI_COMM_WORLD);
+ oocomm->remoteIndices().template rebuild<false>();
+ Dune::graphRepartition(MatrixGraph(*m), *oocomm,
+ oocomm->communicator().size(),
+ oocomm_redist,
+ ri.getInterface(), config.get("verbose", 1)>1);
+ ri.setSetup();
+ oocomm_redist->remoteIndices().template rebuild<false>();
+ redistributeMatrix(*m,*m_redist, *oocomm, *oocomm_redist, ri);
+ std::shared_ptr<Vec> rhs_redist = std::make_shared<Vec>(m_redist->N());
+ ri.redistribute(*rhs, *rhs_redist);
+ m = m_redist;
+ oocomm = oocomm_redist;
+ rhs = rhs_redist;
+ oocomm->copyOwnerToAll(*rhs, *rhs);
+}
+#endif
+
+#endif
diff --git a/bin/solve_mm.ini b/bin/solve_mm.ini
new file mode 100644
index 0000000000000000000000000000000000000000..dcefd4e380e5243a23354a4df183430ca0883028
--- /dev/null
+++ b/bin/solve_mm.ini
@@ -0,0 +1,15 @@
+random_rhs=1
+[solver]
+type=cgsolver
+verbose=2
+maxit=1000
+reduction=1e-7
+[solver.preconditioner]
+type=amg
+verbosity=2
+maxLevel=20
+coarsenTarget=100
+smoother=ssor
+smootherBlocks=4
+n=0
+relaxation=1.0