From eef32851cbb51eb0a0722321081e1009c6577e1d Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Mon, 10 Jun 2013 16:00:19 +0200
Subject: [PATCH] Allows copying AMG and using them independently.
This patch introduces shared_ptrs at various places,
namely for the matrix hierarchies and coarse solvers.
This allows copying an AMG and using the same matrix
hierarchies in both preconditioners. Thus one could
solve the same system for various right hand sides,
where each solve happens in a different thread.
---
dune/istl/paamg/amg.hh | 220 +++++++++++++------------
dune/istl/paamg/hierarchy.hh | 3 +
dune/istl/paamg/smoother.hh | 6 +
dune/istl/paamg/test/CMakeLists.txt | 8 +-
dune/istl/paamg/test/pthreadamgtest.cc | 209 +++++++++++++++++++++++
5 files changed, 338 insertions(+), 108 deletions(-)
create mode 100644 dune/istl/paamg/test/pthreadamgtest.cc
diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh
index 8fbdad3ec..d422b2731 100644
--- a/dune/istl/paamg/amg.hh
+++ b/dune/istl/paamg/amg.hh
@@ -214,6 +214,16 @@ namespace Dune
bool usesDirectCoarseLevelSolver() const;
private:
+ /**
+ * @brief Create matrix and smoother hierarchies.
+ * @param criterion The coarsening criterion.
+ * @param matrix The fine level matrix operator.
+ * @param pinfo The fine level parallel information.
+ */
+ template<class C>
+ void createHierarchies(C& criterion, Operator& matrix,
+ const PI& pinfo);
+
/** @brief Multigrid cycle on a level. */
void mgc();
@@ -246,13 +256,13 @@ namespace Dune
void initIteratorsWithFineLevel();
/** @brief The matrix we solve. */
- OperatorHierarchy* matrices_;
+ Dune::shared_ptr<OperatorHierarchy> matrices_;
/** @brief The arguments to construct the smoother */
SmootherArgs smootherArgs_;
/** @brief The hierarchy of the smoothers. */
- Hierarchy<Smoother,A> smoothers_;
+ Dune::shared_ptr<Hierarchy<Smoother,A> > smoothers_;
/** @brief The solver of the coarsest level. */
- CoarseSolver* solver_;
+ Dune::shared_ptr<CoarseSolver> solver_;
/** @brief The right hand side of our problem. */
Hierarchy<Range,A>* rhs_;
/** @brief The left approximate solution of our problem. */
@@ -263,8 +273,9 @@ namespace Dune
typedef Dune::ScalarProductChooser<X,PI,M::category> ScalarProductChooser;
/** @brief The type of the scalar product for the coarse solver. */
typedef typename ScalarProductChooser::ScalarProduct ScalarProduct;
+ typedef Dune::shared_ptr<ScalarProduct> ScalarProductPointer;
/** @brief Scalar product on the coarse level. */
- ScalarProduct* scalarProduct_;
+ ScalarProductPointer scalarProduct_;
/** @brief Gamma, 1 for V-cycle and 2 for W-cycle. */
std::size_t gamma_;
/** @brief The number of pre and postsmoothing steps. */
@@ -275,7 +286,7 @@ namespace Dune
bool buildHierarchy_;
bool additive;
bool coarsesolverconverged;
- Smoother *coarseSmoother_;
+ Dune::shared_ptr<Smoother> coarseSmoother_;
/** @brief The verbosity level. */
std::size_t verbosity_;
};
@@ -286,7 +297,7 @@ namespace Dune
std::size_t gamma, std::size_t preSmoothingSteps,
std::size_t postSmoothingSteps, bool additive_)
: matrices_(&matrices), smootherArgs_(smootherArgs),
- smoothers_(), solver_(&coarseSolver), scalarProduct_(0),
+ smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver), scalarProduct_(0),
gamma_(gamma), preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(false),
additive(additive_), coarsesolverconverged(true),
coarseSmoother_(), verbosity_(2)
@@ -294,7 +305,7 @@ namespace Dune
assert(matrices_->isBuilt());
// build the necessary smoother hierarchies
- matrices_->coarsenSmoother(smoothers_, smootherArgs_);
+ matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
}
template<class M, class X, class S, class PI, class A>
@@ -302,7 +313,7 @@ namespace Dune
const SmootherArgs& smootherArgs,
const Parameters& parms)
: matrices_(&matrices), smootherArgs_(smootherArgs),
- smoothers_(), solver_(&coarseSolver), scalarProduct_(0),
+ smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver), scalarProduct_(0),
gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
additive(parms.getAdditive()), coarsesolverconverged(true),
@@ -311,7 +322,7 @@ namespace Dune
assert(matrices_->isBuilt());
// build the necessary smoother hierarchies
- matrices_->coarsenSmoother(smoothers_, smootherArgs_);
+ matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
}
template<class M, class X, class S, class PI, class A>
@@ -324,7 +335,7 @@ namespace Dune
bool additive_,
const PI& pinfo)
: smootherArgs_(smootherArgs),
- smoothers_(), solver_(), scalarProduct_(0), gamma_(gamma),
+ smoothers_(new Hierarchy<Smoother,A>), solver_(), scalarProduct_(0), gamma_(gamma),
preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(true),
additive(additive_), coarsesolverconverged(true),
coarseSmoother_(), verbosity_(criterion.debugLevel())
@@ -334,16 +345,7 @@ namespace Dune
// TODO: reestablish compile time checks.
//dune_static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
// "Matrix and Solver must match in terms of category!");
- Timer watch;
- matrices_ = new OperatorHierarchy(const_cast<Operator&>(matrix), pinfo);
-
- matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
-
- // build the necessary smoother hierarchies
- matrices_->coarsenSmoother(smoothers_, smootherArgs_);
-
- if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
- std::cout<<"Building Hierarchy of "<<matrices_->maxlevels()<<" levels took "<<watch.elapsed()<<" seconds."<<std::endl;
+ createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
}
template<class M, class X, class S, class PI, class A>
@@ -353,7 +355,7 @@ namespace Dune
const SmootherArgs& smootherArgs,
const PI& pinfo)
: smootherArgs_(smootherArgs),
- smoothers_(), solver_(), scalarProduct_(0),
+ smoothers_(new Hierarchy<Smoother,A>), solver_(), scalarProduct_(0),
gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
additive(criterion.getAdditive()), coarsesolverconverged(true),
@@ -364,26 +366,91 @@ namespace Dune
// TODO: reestablish compile time checks.
//dune_static_assert(static_cast<int>(PI::category)==static_cast<int>(S::category),
// "Matrix and Solver must match in terms of category!");
+ createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
+ }
+
+ template<class M, class X, class S, class PI, class A>
+ AMG<M,X,S,PI,A>::~AMG()
+ {
+ if(buildHierarchy_) {
+ if(solver_)
+ solver_.reset();
+ if(coarseSmoother_)
+ coarseSmoother_.reset();
+ }
+ }
+
+ template<class M, class X, class S, class PI, class A>
+ template<class C>
+ void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, Operator& matrix,
+ const PI& pinfo)
+ {
Timer watch;
- matrices_ = new OperatorHierarchy(const_cast<Operator&>(matrix), pinfo);
+ matrices_.reset(new OperatorHierarchy(matrix, pinfo));
matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
// build the necessary smoother hierarchies
- matrices_->coarsenSmoother(smoothers_, smootherArgs_);
+ matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
std::cout<<"Building Hierarchy of "<<matrices_->maxlevels()<<" levels took "<<watch.elapsed()<<" seconds."<<std::endl;
- }
- template<class M, class X, class S, class PI, class A>
- AMG<M,X,S,PI,A>::~AMG()
- {
- if(buildHierarchy_) {
- delete matrices_;
+ if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
+ // We have the carsest level. Create the coarse Solver
+ SmootherArgs sargs(smootherArgs_);
+ sargs.iterations = 1;
+
+ typename ConstructionTraits<Smoother>::Arguments cargs;
+ cargs.setArgs(sargs);
+ if(matrices_->redistributeInformation().back().isSetup()) {
+ // Solve on the redistributed partitioning
+ cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
+ cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
+ }else{
+ cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
+ cargs.setComm(*matrices_->parallelInformation().coarsest());
+ }
+
+ coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
+ scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
+
+#if HAVE_SUPERLU
+ // Use superlu if we are purely sequential or with only one processor on the coarsest level.
+ if(is_same<ParallelInformation,SequentialInformation>::value // sequential mode
+ || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
+ || (matrices_->parallelInformation().coarsest().isRedistributed()
+ && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
+ && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) { // redistribute and 1 proc
+ if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
+ std::cout<<"Using superlu"<<std::endl;
+ if(matrices_->parallelInformation().coarsest().isRedistributed())
+ {
+ if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
+ // We are still participating on this level
+ solver_.reset(new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
+ else
+ solver_.reset();
+ }else
+ solver_.reset(new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(), false, false));
+ }else
+#endif
+ {
+ if(matrices_->parallelInformation().coarsest().isRedistributed())
+ {
+ if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
+ // We are still participating on this level
+ solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
+ *scalarProduct_,
+ *coarseSmoother_, 1E-2, 1000, 0));
+ else
+ solver_.reset();
+ }else
+ solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
+ *scalarProduct_,
+ *coarseSmoother_, 1E-2, 1000, 0));
+ }
}
- if(scalarProduct_)
- delete scalarProduct_;
}
@@ -419,8 +486,8 @@ namespace Dune
diagonal.solve(x[row.index()], b[row.index()]);
}
- if(smoothers_.levels()>0)
- smoothers_.finest()->pre(x,b);
+ if(smoothers_->levels()>0)
+ smoothers_->finest()->pre(x,b);
else
// No smoother to make x consistent! Do it by hand
matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
@@ -438,15 +505,15 @@ namespace Dune
typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
typedef typename Hierarchy<Range,A>::Iterator RIterator;
typedef typename Hierarchy<Domain,A>::Iterator DIterator;
- Iterator coarsest = smoothers_.coarsest();
- Iterator smoother = smoothers_.finest();
+ Iterator coarsest = smoothers_->coarsest();
+ Iterator smoother = smoothers_->finest();
RIterator rhs = rhs_->finest();
DIterator lhs = lhs_->finest();
- if(smoothers_.levels()>0) {
+ if(smoothers_->levels()>0) {
assert(lhs_->levels()==rhs_->levels());
- assert(smoothers_.levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
- assert(smoothers_.levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
+ assert(smoothers_->levels()==lhs_->levels() || matrices_->levels()==matrices_->maxlevels());
+ assert(smoothers_->levels()+1==lhs_->levels() || matrices_->levels()<matrices_->maxlevels());
if(smoother!=coarsest)
for(++smoother, ++lhs, ++rhs; smoother != coarsest; ++smoother, ++lhs, ++rhs)
@@ -460,63 +527,6 @@ namespace Dune
x = *lhs_->finest();
b = *rhs_->finest();
- if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
- // We have the carsest level. Create the coarse Solver
- SmootherArgs sargs(smootherArgs_);
- sargs.iterations = 1;
-
- typename ConstructionTraits<Smoother>::Arguments cargs;
- cargs.setArgs(sargs);
- if(matrices_->redistributeInformation().back().isSetup()) {
- // Solve on the redistributed partitioning
- cargs.setMatrix(matrices_->matrices().coarsest().getRedistributed().getmat());
- cargs.setComm(matrices_->parallelInformation().coarsest().getRedistributed());
-
- coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
- scalarProduct_ = ScalarProductChooser::construct(matrices_->parallelInformation().coarsest().getRedistributed());
- }else{
- cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
- cargs.setComm(*matrices_->parallelInformation().coarsest());
-
- coarseSmoother_ = ConstructionTraits<Smoother>::construct(cargs);
- scalarProduct_ = ScalarProductChooser::construct(*matrices_->parallelInformation().coarsest());
- }
-#if HAVE_SUPERLU
- // Use superlu if we are purely sequential or with only one processor on the coarsest level.
- if(is_same<ParallelInformation,SequentialInformation>::value // sequential mode
- || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
- || (matrices_->parallelInformation().coarsest().isRedistributed()
- && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
- && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) { // redistribute and 1 proc
- if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
- std::cout<<"Using superlu"<<std::endl;
- if(matrices_->parallelInformation().coarsest().isRedistributed())
- {
- if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
- // We are still participating on this level
- solver_ = new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat());
- else
- solver_ = 0;
- }else
- solver_ = new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat());
- }else
-#endif
- {
- if(matrices_->parallelInformation().coarsest().isRedistributed())
- {
- if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
- // We are still participating on this level
- solver_ = new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
- *scalarProduct_,
- *coarseSmoother_, 1E-2, 10000, 0);
- else
- solver_ = 0;
- }else
- solver_ = new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
- *scalarProduct_,
- *coarseSmoother_, 1E-2, 1000, 0);
- }
- }
}
template<class M, class X, class S, class PI, class A>
std::size_t AMG<M,X,S,PI,A>::levels()
@@ -560,7 +570,7 @@ namespace Dune
template<class M, class X, class S, class PI, class A>
void AMG<M,X,S,PI,A>::initIteratorsWithFineLevel()
{
- smoother = smoothers_.finest();
+ smoother = smoothers_->finest();
matrix = matrices_->matrices().finest();
pinfo = matrices_->parallelInformation().finest();
redist =
@@ -767,7 +777,7 @@ namespace Dune
//pinfo = matrices_->parallelInformation().coarsest
// calculate correction for all levels
lhs = lhs_->finest();
- typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_.finest();
+ typename Hierarchy<Smoother,A>::Iterator smoother = smoothers_->finest();
for(rhs=rhs_->finest(); rhs != rhs_->coarsest(); ++lhs, ++rhs, ++smoother) {
// presmoothing
@@ -801,20 +811,14 @@ namespace Dune
template<class M, class X, class S, class PI, class A>
void AMG<M,X,S,PI,A>::post(Domain& x)
{
- if(buildHierarchy_) {
- if(solver_)
- delete solver_;
- if(coarseSmoother_)
- ConstructionTraits<Smoother>::deconstruct(coarseSmoother_);
- }
// Postprocess all smoothers
typedef typename Hierarchy<Smoother,A>::Iterator Iterator;
typedef typename Hierarchy<Domain,A>::Iterator DIterator;
- Iterator coarsest = smoothers_.coarsest();
- Iterator smoother = smoothers_.finest();
+ Iterator coarsest = smoothers_->coarsest();
+ Iterator smoother = smoothers_->finest();
DIterator lhs = lhs_->finest();
- if(smoothers_.levels()>0) {
+ if(smoothers_->levels()>0) {
if(smoother != coarsest || matrices_->levels()<matrices_->maxlevels())
smoother->post(*lhs);
if(smoother!=coarsest)
@@ -822,13 +826,15 @@ namespace Dune
smoother->post(*lhs);
smoother->post(*lhs);
}
-
delete &(*lhs_->finest());
delete lhs_;
+ lhs_=nullptr;
delete &(*update_->finest());
delete update_;
+ update_=nullptr;
delete &(*rhs_->finest());
delete rhs_;
+ rhs_=nullptr;
}
template<class M, class X, class S, class PI, class A>
diff --git a/dune/istl/paamg/hierarchy.hh b/dune/istl/paamg/hierarchy.hh
index a1d54c047..dbf9ae32e 100644
--- a/dune/istl/paamg/hierarchy.hh
+++ b/dune/istl/paamg/hierarchy.hh
@@ -282,6 +282,9 @@ namespace Dune
~Hierarchy();
private:
+ /** @brief Deactivate copying */
+ Hierarchy(const Hierarchy&)
+ {}
/** @brief The finest element in the hierarchy. */
Element* finest_;
/** @brief The coarsest element in the hierarchy. */
diff --git a/dune/istl/paamg/smoother.hh b/dune/istl/paamg/smoother.hh
index 857d91075..e117014ea 100644
--- a/dune/istl/paamg/smoother.hh
+++ b/dune/istl/paamg/smoother.hh
@@ -119,6 +119,11 @@ namespace Dune
void setComm(T1& comm)
{}
+ const SequentialInformation& getComm()
+ {
+ return comm_;
+ }
+
const SmootherArgs getArgs() const
{
return *args_;
@@ -128,6 +133,7 @@ namespace Dune
const Matrix* matrix_;
private:
const SmootherArgs* args_;
+ SequentialInformation comm_;
};
template<class T>
diff --git a/dune/istl/paamg/test/CMakeLists.txt b/dune/istl/paamg/test/CMakeLists.txt
index f197caf57..b9dad7060 100644
--- a/dune/istl/paamg/test/CMakeLists.txt
+++ b/dune/istl/paamg/test/CMakeLists.txt
@@ -6,7 +6,7 @@ if(PARMETIS_FOUND)
set(PARMETISTESTS pamgtest pamg_comm_repart_test)
endif(PARMETIS_FOUND)
-set(NORMALTESTS amgtest fastamg graphtest kamgtest)
+set(NORMALTESTS amgtest fastamg graphtest kamgtest pthreadamgtest)
set(ALLTESTS ${MPITESTS} ${PARMETISTESTS} ${NORMALTESTS})
# We do not want want to build the tests during make all,
@@ -14,6 +14,12 @@ set(ALLTESTS ${MPITESTS} ${PARMETISTESTS} ${NORMALTESTS})
add_directory_test_target(_test_target)
add_dependencies(${_test_target} ${ALLTESTS})
+find_package(Threads)
+if(CMAKE_USE_PTHREADS_INIT)
+ add_executable(pthreadamgtest pthreadamgtest.cc)
+ target_link_libraries(pthreadamgtest "${CMAKE_THREAD_LIBS_INIT}")
+endif(CMAKE_USE_PTHREADS_INIT)
+
add_executable(amgtest "amgtest.cc")
add_dune_superlu_flags(amgtest)
add_executable(fastamg "fastamg.cc")
diff --git a/dune/istl/paamg/test/pthreadamgtest.cc b/dune/istl/paamg/test/pthreadamgtest.cc
new file mode 100644
index 000000000..f43c9e106
--- /dev/null
+++ b/dune/istl/paamg/test/pthreadamgtest.cc
@@ -0,0 +1,209 @@
+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#include "config.h"
+
+/*#include"xreal.h"
+
+ namespace std
+ {
+
+ HPA::xreal abs(const HPA::xreal& t)
+ {
+ return t>=0?t:-t;
+ }
+
+ };
+ */
+
+#include "anisotropic.hh"
+#include <dune/common/timer.hh>
+#include <dune/common/parallel/indexset.hh>
+#include <dune/common/parallel/collectivecommunication.hh>
+#include <dune/istl/paamg/amg.hh>
+#include <dune/istl/paamg/pinfo.hh>
+#include <dune/istl/solvers.hh>
+#include <cstdlib>
+#include <ctime>
+#include <pthread.h>
+#define NUM_THREADS 5
+
+typedef double XREAL;
+
+/*
+ typedef HPA::xreal XREAL;
+
+ namespace Dune
+ {
+ template<>
+ struct DoubleConverter<HPA::xreal>
+ {
+ static double toDouble(const HPA::xreal& t)
+ {
+ return t._2double();
+ }
+ };
+ }
+ */
+
+template<class M, class V>
+void randomize(const M& mat, V& b)
+{
+ V x=b;
+
+ srand((unsigned)std::clock());
+
+ typedef typename V::iterator iterator;
+ for(iterator i=x.begin(); i != x.end(); ++i)
+ *i=(rand() / (RAND_MAX + 1.0));
+
+ mat.mv(static_cast<const V&>(x), b);
+}
+
+typedef Dune::ParallelIndexSet<int,LocalIndex,512> ParallelIndexSet;
+typedef Dune::FieldMatrix<XREAL,1,1> MatrixBlock;
+typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
+typedef Dune::FieldVector<XREAL,1> VectorBlock;
+typedef Dune::BlockVector<VectorBlock> Vector;
+typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
+typedef Dune::CollectiveCommunication<void*> Comm;
+typedef Dune::SeqSSOR<BCRSMat,Vector,Vector> Smoother;
+typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
+typedef Dune::Amg::AMG<Operator,Vector,Smoother> AMG;
+
+struct thread_arg
+{
+ AMG *amg;
+ Vector *b;
+ Vector *x;
+ Operator *fop;
+};
+
+
+void *solve(void* arg)
+{
+ thread_arg *amgarg=(thread_arg*) arg;
+
+ Dune::GeneralizedPCGSolver<Vector> amgCG(*amgarg->fop,*amgarg->amg,1e-6,80,2);
+ //Dune::LoopSolver<Vector> amgCG(fop, amg, 1e-4, 10000, 2);
+ Dune::Timer watch;
+ Dune::InverseOperatorResult r;
+ amgCG.apply(*amgarg->x,*amgarg->b,r);
+
+ double solvetime = watch.elapsed();
+
+ std::cout<<"AMG solving took "<<solvetime<<" seconds"<<std::endl;
+
+ pthread_exit(NULL);
+}
+
+
+template <int BS>
+void testAMG(int N, int coarsenTarget, int ml)
+{
+
+ std::cout<<"N="<<N<<" coarsenTarget="<<coarsenTarget<<" maxlevel="<<ml<<std::endl;
+
+
+
+ ParallelIndexSet indices;
+ int n;
+
+ Comm c;
+ BCRSMat mat = setupAnisotropic2d<1,XREAL>(N, indices, c, &n, 1);
+
+ Vector b(mat.N()), x(mat.M());
+
+ b=0;
+ x=100;
+
+ setBoundary(x, b, N);
+
+ x=0;
+ randomize(mat, b);
+
+ std::vector<Vector> xs(NUM_THREADS, x);
+ std::vector<Vector> bs(NUM_THREADS, b);
+
+ if(N<6) {
+ Dune::printmatrix(std::cout, mat, "A", "row");
+ Dune::printvector(std::cout, x, "x", "row");
+ }
+
+ Dune::Timer watch;
+
+ watch.reset();
+ Operator fop(mat);
+
+ typedef Dune::Amg::CoarsenCriterion<Dune::Amg::UnSymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
+ Criterion;
+ //typedef Dune::SeqSOR<BCRSMat,Vector,Vector> Smoother;
+ //typedef Dune::SeqJac<BCRSMat,Vector,Vector> Smoother;
+ //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::MultiplicativeSchwarzMode> Smoother;
+ //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector,Dune::SymmetricMultiplicativeSchwarzMode> Smoother;
+ //typedef Dune::SeqOverlappingSchwarz<BCRSMat,Vector> Smoother;
+
+ SmootherArgs smootherArgs;
+
+ smootherArgs.iterations = 1;
+
+ //smootherArgs.overlap=SmootherArgs::vertex;
+ //smootherArgs.overlap=SmootherArgs::none;
+ //smootherArgs.overlap=SmootherArgs::aggregate;
+
+ smootherArgs.relaxationFactor = 1;
+
+ Criterion criterion(15,coarsenTarget);
+ criterion.setDefaultValuesIsotropic(2);
+ criterion.setAlpha(.67);
+ criterion.setBeta(1.0e-4);
+ criterion.setMaxLevel(ml);
+ criterion.setSkipIsolated(false);
+
+ Dune::SeqScalarProduct<Vector> sp;
+
+ Smoother smoother(mat,1,1);
+
+ AMG amg(fop, criterion, smootherArgs, 1, 1, 1, false);
+
+
+ double buildtime = watch.elapsed();
+
+ std::cout<<"Building hierarchy took "<<buildtime<<" seconds"<<std::endl;
+ std::vector<AMG> amgs(NUM_THREADS, amg);
+ std::vector<thread_arg> args(NUM_THREADS);
+ std::vector<pthread_t> threads(NUM_THREADS);
+ for(int i=0; i < NUM_THREADS; ++i)
+ {
+ args[i].amg=&amgs[i];
+ args[i].b=&bs[i];
+ args[i].x=&xs[i];
+ args[i].fop=&fop;
+ pthread_create(&threads[i], NULL, solve, (void*) &args[i]);
+ }
+ void* retval;
+
+ for(int i=0; i < NUM_THREADS; ++i)
+ pthread_join(threads[i], &retval);
+}
+
+
+int main(int argc, char** argv)
+{
+
+ int N=100;
+ int coarsenTarget=1200;
+ int ml=10;
+
+ if(argc>1)
+ N = atoi(argv[1]);
+
+ if(argc>2)
+ coarsenTarget = atoi(argv[2]);
+
+ if(argc>3)
+ ml = atoi(argv[3]);
+
+ testAMG<1>(N, coarsenTarget, ml);
+ //testAMG<2>(N, coarsenTarget, ml);
+
+}
--
GitLab