From 3da29f78907af49941719fb62a2a74cde37186ef Mon Sep 17 00:00:00 2001
From: Markus Blatt <markus@dr-blatt.de>
Date: Wed, 19 Jun 2013 12:38:00 +0200
Subject: [PATCH] Support copying of FastAMG and usage in different threads.
---
dune/istl/paamg/amg.hh | 2 +-
dune/istl/paamg/fastamg.hh | 210 +++++++++++++++----------
dune/istl/paamg/test/CMakeLists.txt | 12 +-
dune/istl/paamg/test/pthreadamgtest.cc | 16 +-
4 files changed, 148 insertions(+), 92 deletions(-)
diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh
index 5b347d798..4dc797628 100644
--- a/dune/istl/paamg/amg.hh
+++ b/dune/istl/paamg/amg.hh
@@ -169,7 +169,7 @@ namespace Dune
*/
template<class C>
AMG(const Operator& fineOperator, const C& criterion,
- const SmootherArgs& smootherArgs,
+ const SmootherArgs& smootherArgs=SmootherArgs(),
const ParallelInformation& pinfo=ParallelInformation());
/**
diff --git a/dune/istl/paamg/fastamg.hh b/dune/istl/paamg/fastamg.hh
index ea1f93be2..bd6309fb2 100644
--- a/dune/istl/paamg/fastamg.hh
+++ b/dune/istl/paamg/fastamg.hh
@@ -55,7 +55,7 @@ namespace Dune
* \tparam A An allocator for X
*/
template<class M, class X, class PI=SequentialInformation, class A=std::allocator<X> >
- class FastAMG : public Dune::Preconditioner<X,X>
+ class FastAMG : public Preconditioner<X,X>
{
public:
/** @brief The matrix operator type. */
@@ -68,7 +68,7 @@ namespace Dune
*/
typedef PI ParallelInformation;
/** @brief The operator hierarchy type. */
- typedef Dune::Amg::MatrixHierarchy<M, ParallelInformation, A> OperatorHierarchy;
+ typedef MatrixHierarchy<M, ParallelInformation, A> OperatorHierarchy;
/** @brief The parallal data distribution hierarchy type. */
typedef typename OperatorHierarchy::ParallelInformationHierarchy ParallelInformationHierarchy;
@@ -77,7 +77,7 @@ namespace Dune
/** @brief The range type. */
typedef X Range;
/** @brief the type of the coarse solver. */
- typedef Dune::InverseOperator<X,X> CoarseSolver;
+ typedef InverseOperator<X,X> CoarseSolver;
enum {
/** @brief The solver category. */
@@ -92,7 +92,7 @@ namespace Dune
* @param parms The parameters for the AMG.
*/
FastAMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
- const Dune::Amg::Parameters& parms,
+ const Parameters& parms,
bool symmetric=true);
/**
@@ -108,10 +108,15 @@ namespace Dune
*/
template<class C>
FastAMG(const Operator& fineOperator, const C& criterion,
- const Parameters& parms,
+ const Parameters& parms=Parameters(),
bool symmetric=true,
const ParallelInformation& pinfo=ParallelInformation());
+ /**
+ * @brief Copy constructor.
+ */
+ FastAMG(const FastAMG& amg),
+
~FastAMG();
/** \copydoc Preconditioner::pre */
@@ -144,7 +149,7 @@ namespace Dune
*/
void recalculateHierarchy()
{
- matrices_->recalculateGalerkin(Dune::NegateSet<typename PI::OwnerSet>());
+ matrices_->recalculateGalerkin(NegateSet<typename PI::OwnerSet>());
}
/**
@@ -154,6 +159,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(Domain& x, const Range& b);
@@ -161,10 +176,9 @@ namespace Dune
typename ParallelInformationHierarchy::Iterator pinfo;
typename OperatorHierarchy::RedistributeInfoList::const_iterator redist;
typename OperatorHierarchy::AggregatesMapList::const_iterator aggregates;
- typename Dune::Amg::Hierarchy<Domain,A>::Iterator lhs;
- typename Dune::Amg::Hierarchy<Domain,A>::Iterator residual;
- typename Dune::Amg::Hierarchy<Domain,A>::Iterator tmp;
- typename Dune::Amg::Hierarchy<Range,A>::Iterator rhs;
+ typename Hierarchy<Domain,A>::Iterator lhs;
+ typename Hierarchy<Domain,A>::Iterator residual;
+ typename Hierarchy<Range,A>::Iterator rhs;
/** @brief Apply pre smoothing on the current level. */
@@ -189,17 +203,16 @@ namespace Dune
/** @brief The solver of the coarsest level. */
Dune::shared_ptr<CoarseSolver> solver_;
/** @brief The right hand side of our problem. */
- Dune::Amg::Hierarchy<Range,A>* rhs_;
+ Hierarchy<Range,A>* rhs_;
/** @brief The left approximate solution of our problem. */
- Dune::Amg::Hierarchy<Domain,A>* lhs_;
+ Hierarchy<Domain,A>* lhs_;
/** @brief The current residual. */
- Dune::Amg::Hierarchy<Domain,A>* residual_;
+ Hierarchy<Domain,A>* residual_;
/** @brief The type of the chooser of the scalar product. */
- typedef Dune::ScalarProductChooser<X,PI,M::category> ScalarProductChooser;
+ typedef ScalarProductChooser<X,PI,M::category> ScalarProductChooserType;
/** @brief The type of the scalar product for the coarse solver. */
- typedef typename ScalarProductChooser::ScalarProduct ScalarProduct;
-
+ typedef typename ScalarProductChooserType::ScalarProduct ScalarProduct;
typedef Dune::shared_ptr<ScalarProduct> ScalarProductPointer;
/** @brief Scalar product on the coarse level. */
ScalarProductPointer scalarProduct_;
@@ -213,17 +226,34 @@ namespace Dune
bool buildHierarchy_;
bool symmetric;
bool coarsesolverconverged;
- typedef Dune::SeqSSOR<typename M::matrix_type,X,X> Smoother;
+ typedef SeqSSOR<typename M::matrix_type,X,X> Smoother;
typedef Dune::shared_ptr<Smoother> SmootherPointer;
SmootherPointer coarseSmoother_;
/** @brief The verbosity level. */
std::size_t verbosity_;
};
+ template<class M, class X, class PI, class A>
+ FastAMG<M,X,PI,A>::FastAMG(const FastAMG& amg)
+ : matrices_(amg.matrices_), solver_(amg.solver_),
+ rhs_(), lhs_(), residual_(), scalarProduct_(amg.scalarProduct_),
+ gamma_(amg.gamma_), preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
+ symmetric(amg.symmetric), coarsesolverconverged(amg.coarsesolverconverged),
+ coarseSmoother_(amg.coarseSmoother_), verbosity_(amg.verbosity_)
+ {
+ if(amg.rhs_)
+ rhs_=new Hierarchy<Range,A>(*amg.rhs_);
+ if(amg.lhs_)
+ lhs_=new Hierarchy<Domain,A>(*amg.lhs_);
+ if(amg.residual_)
+ residual_=new Hierarchy<Domain,A>(*amg.residual_);
+ }
+
template<class M, class X, class PI, class A>
FastAMG<M,X,PI,A>::FastAMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
- const Dune::Amg::Parameters& parms, bool symmetric_)
- : matrices_(&matrices), solver_(&coarseSolver), scalarProduct_(),
+ const Parameters& parms, bool symmetric_)
+ : matrices_(&matrices), solver_(&coarseSolver),
+ rhs_(), lhs_(), residual_(), scalarProduct_(),
gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
symmetric(symmetric_), coarsesolverconverged(true),
@@ -235,16 +265,16 @@ namespace Dune
preSteps_=postSteps_=0;
}
assert(matrices_->isBuilt());
- dune_static_assert((Dune::is_same<PI,Dune::Amg::SequentialInformation>::value), "Currently only sequential runs are supported");
+ dune_static_assert((is_same<PI,SequentialInformation>::value), "Currently only sequential runs are supported");
}
template<class M, class X, class PI, class A>
template<class C>
FastAMG<M,X,PI,A>::FastAMG(const Operator& matrix,
const C& criterion,
- const Dune::Amg::Parameters& parms,
+ const Parameters& parms,
bool symmetric_,
const PI& pinfo)
- : solver_(), scalarProduct_(), gamma_(parms.getGamma()),
+ : solver_(), rhs_(), lhs_(), residual_(), scalarProduct_(), gamma_(parms.getGamma()),
preSteps_(parms.getNoPreSmoothSteps()), postSteps_(parms.getNoPostSmoothSteps()),
buildHierarchy_(true),
symmetric(symmetric_), coarsesolverconverged(true),
@@ -255,42 +285,69 @@ namespace Dune
std::cerr<<"WARNING only one step of smoothing is supported!"<<std::endl;
preSteps_=postSteps_=1;
}
- dune_static_assert((Dune::is_same<PI,Dune::Amg::SequentialInformation>::value), "Currently only sequential runs are supported");
+ dune_static_assert((is_same<PI,SequentialInformation>::value), "Currently only sequential runs are supported");
// 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!");
- Dune::Timer watch;
- matrices_.reset(new OperatorHierarchy(const_cast<Operator&>(matrix), pinfo));
+ createHierarchies(criterion, const_cast<Operator&>(matrix), pinfo);
+ }
+
+ template<class M, class X, class PI, class A>
+ FastAMG<M,X,PI,A>::~FastAMG()
+ {
+ if(buildHierarchy_) {
+ if(solver_)
+ solver_.reset();
+ if(coarseSmoother_)
+ coarseSmoother_.reset();
+ }
+ if(lhs_)
+ delete lhs_;
+ lhs_=nullptr;
+ if(residual_)
+ delete residual_;
+ residual_=nullptr;
+ if(rhs_)
+ delete rhs_;
+ rhs_=nullptr;
+ }
- matrices_->template build<Dune::NegateSet<typename PI::OwnerSet> >(criterion);
+ template<class M, class X, class PI, class A>
+ template<class C>
+ void FastAMG<M,X,PI,A>::createHierarchies(C& criterion, Operator& matrix,
+ const PI& pinfo)
+ {
+ Timer watch;
+ matrices_.reset(new OperatorHierarchy(matrix, pinfo));
+
+ matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
+
+ if(verbosity_>0 && matrices_->parallelInformation().finest()->communicator().rank()==0)
+ std::cout<<"Building Hierarchy of "<<matrices_->maxlevels()<<" levels took "<<watch.elapsed()<<" seconds."<<std::endl;
if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
// We have the carsest level. Create the coarse Solver
- typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
+ typedef typename SmootherTraits<Smoother>::Arguments SmootherArgs;
SmootherArgs sargs;
sargs.iterations = 1;
- typename Dune::Amg::ConstructionTraits<Smoother>::Arguments cargs;
+ 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_ = SmootherPointer(Dune::Amg::ConstructionTraits<Smoother>::construct(cargs),
- Dune::Amg::ConstructionTraits<Smoother>::deconstruct);
- scalarProduct_ = ScalarProductPointer(ScalarProductChooser::construct(matrices_->parallelInformation().coarsest().getRedistributed()));
}else{
cargs.setMatrix(matrices_->matrices().coarsest()->getmat());
cargs.setComm(*matrices_->parallelInformation().coarsest());
-
- coarseSmoother_ = SmootherPointer(Dune::Amg::ConstructionTraits<Smoother>::construct(cargs),
- Dune::Amg::ConstructionTraits<Smoother>::deconstruct);
- scalarProduct_ = ScalarProductPointer(ScalarProductChooser::construct(*matrices_->parallelInformation().coarsest()));
}
+
+ coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
+ scalarProduct_.reset(ScalarProductChooserType::construct(cargs.getComm()));
+
#if HAVE_SUPERLU
// Use superlu if we are purely sequential or with only one processor on the coarsest level.
- if(Dune::is_same<ParallelInformation,Dune::Amg::SequentialInformation>::value // sequential mode
+ 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
@@ -301,9 +358,11 @@ namespace Dune
{
if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
// We are still participating on this level
- solver_ = Dune::shared_ptr<Dune::SuperLU<typename M::matrix_type> >( new Dune::SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat()));
+ solver_.reset(new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
+ else
+ solver_.reset();
}else
- solver_ = Dune::shared_ptr<Dune::SuperLU<typename M::matrix_type> >( new Dune::SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat()));
+ solver_.reset(new SuperLU<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(), false, false));
}else
#endif
{
@@ -311,13 +370,15 @@ namespace Dune
{
if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
// We are still participating on this level
- solver_ = Dune::shared_ptr<Dune::BiCGSTABSolver<X> >(new Dune::BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
- *scalarProduct_,
- *coarseSmoother_, 1E-2, 10000, 0));
+ solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(matrices_->matrices().coarsest().getRedistributed()),
+ *scalarProduct_,
+ *coarseSmoother_, 1E-2, 1000, 0));
+ else
+ solver_.reset();
}else
- solver_ = Dune::shared_ptr<Dune::BiCGSTABSolver<X> >(new Dune::BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
- *scalarProduct_,
- *coarseSmoother_, 1E-2, 1000, 0));
+ solver_.reset(new BiCGSTABSolver<X>(const_cast<M&>(*matrices_->matrices().coarsest()),
+ *scalarProduct_,
+ *coarseSmoother_, 1E-2, 1000, 0));
}
}
@@ -325,15 +386,11 @@ namespace Dune
std::cout<<"Building Hierarchy of "<<matrices_->maxlevels()<<" levels took "<<watch.elapsed()<<" seconds."<<std::endl;
}
- template<class M, class X, class PI, class A>
- FastAMG<M,X,PI,A>::~FastAMG()
- {}
-
template<class M, class X, class PI, class A>
void FastAMG<M,X,PI,A>::pre(Domain& x, Range& b)
{
- Dune::Timer watch, watch1;
+ Timer watch, watch1;
// Detect Matrix rows where all offdiagonal entries are
// zero and set x such that A_dd*x_d=b_d
// Thus users can be more careless when setting up their linear
@@ -367,11 +424,15 @@ namespace Dune
// No smoother to make x consistent! Do it by hand
matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
Range* copy = new Range(b);
- rhs_ = new Dune::Amg::Hierarchy<Range,A>(*copy);
+ if(rhs_)
+ delete rhs_;
+ rhs_ = new Hierarchy<Range,A>(copy);
Domain* dcopy = new Domain(x);
- lhs_ = new Dune::Amg::Hierarchy<Domain,A>(*dcopy);
+ if(lhs_)
+ delete lhs_;
+ lhs_ = new Hierarchy<Domain,A>(dcopy);
dcopy = new Domain(x);
- residual_ = new Dune::Amg::Hierarchy<Domain,A>(*dcopy);
+ residual_ = new Hierarchy<Domain,A>(dcopy);
matrices_->coarsenVector(*rhs_);
matrices_->coarsenVector(*lhs_);
matrices_->coarsenVector(*residual_);
@@ -380,10 +441,6 @@ namespace Dune
// copy the changes to the original vectors.
x = *lhs_->finest();
b = *rhs_->finest();
-
- std::cout<<" Preprocessing smoothers took "<<watch1.elapsed()<<std::endl;
- watch1.reset();
- std::cout<<"AMG::pre took "<<watch.elapsed()<<std::endl;
}
template<class M, class X, class PI, class A>
std::size_t FastAMG<M,X,PI,A>::levels()
@@ -406,20 +463,10 @@ namespace Dune
assert(v.two_norm()==0);
level=0;
- /*if(redist->isSetup())
- {
- *lhs=v;
- *rhs = d;
- mgc(*lhs, *rhs);
- if(postSteps_==0||matrices_->maxlevels()==1)
- pinfo->copyOwnerToAll(*lhs, *lhs);
- v=*lhs;
- }else{*/
+
mgc(v, d);
if(postSteps_==0||matrices_->maxlevels()==1)
pinfo->copyOwnerToAll(v, v);
- //}
-
}
template<class M, class X, class PI, class A>
@@ -448,14 +495,14 @@ namespace Dune
if(processNextLevel) {
//restrict defect to coarse level right hand side.
++pinfo;
- Dune::Amg::Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
+ Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(residual.getRedistributed()), *pinfo);
}
}else{
//restrict defect to coarse level right hand side.
- typename Dune::Amg::Hierarchy<Range,A>::Iterator fineRhs = rhs++;
+ typename Hierarchy<Range,A>::Iterator fineRhs = rhs++;
++pinfo;
- Dune::Amg::Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
+ Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::restrictVector(*(*aggregates), *rhs, static_cast<const Range&>(*residual), *pinfo);
}
@@ -495,15 +542,15 @@ namespace Dune
}
- typename Dune::Amg::Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
+ typename Hierarchy<Domain,A>::Iterator coarseLhs = lhs--;
if(redist->isSetup()) {
// Need to redistribute during prolongate
- Dune::Amg::Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
+ Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::prolongateVector(*(*aggregates), *coarseLhs, x, lhs.getRedistributed(), matrices_->getProlongationDampingFactor(),
*pinfo, *redist);
}else{
- Dune::Amg::Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
+ Transfer<typename OperatorHierarchy::AggregatesMap::AggregateDescriptor,Range,ParallelInformation>
::prolongateVector(*(*aggregates), *coarseLhs, x,
matrices_->getProlongationDampingFactor(), *pinfo);
@@ -540,7 +587,7 @@ namespace Dune
template<class M, class X, class PI, class A>
bool FastAMG<M,X,PI,A>::usesDirectCoarseLevelSolver() const
{
- return Dune::IsDirectSolver< CoarseSolver>::value;
+ return IsDirectSolver< CoarseSolver>::value;
}
template<class M, class X, class PI, class A>
@@ -548,7 +595,7 @@ namespace Dune
if(matrix == matrices_->matrices().coarsest() && levels()==maxlevels()) {
// Solve directly
- Dune::InverseOperatorResult res;
+ InverseOperatorResult res;
res.converged=true; // If we do not compute this flag will not get updated
if(redist->isSetup()) {
redist->redistribute(b, rhs.getRedistributed());
@@ -590,7 +637,7 @@ namespace Dune
if(matrix == matrices_->matrices().finest()) {
coarsesolverconverged = matrices_->parallelInformation().finest()->communicator().prod(coarsesolverconverged);
if(!coarsesolverconverged)
- DUNE_THROW(Dune::MathError, "Coarse solver did not converge");
+ DUNE_THROW(MathError, "Coarse solver did not converge");
}
// printvector(std::cout, *lhs, "update corrected", "u", 10, 10, 10);
@@ -606,13 +653,12 @@ namespace Dune
template<class M, class X, class PI, class A>
void FastAMG<M,X,PI,A>::post(Domain& x)
{
- delete &(*lhs_->finest());
delete lhs_;
- delete &(*residual_->finest());
- delete residual_;
-
- delete &(*rhs_->finest());
+ lhs_=nullptr;
delete rhs_;
+ rhs_=nullptr;
+ delete residual_;
+ residual_=nullptr;
}
template<class M, class X, class PI, class A>
diff --git a/dune/istl/paamg/test/CMakeLists.txt b/dune/istl/paamg/test/CMakeLists.txt
index b9dad7060..247301c0f 100644
--- a/dune/istl/paamg/test/CMakeLists.txt
+++ b/dune/istl/paamg/test/CMakeLists.txt
@@ -5,8 +5,10 @@ endif(MPI_FOUND)
if(PARMETIS_FOUND)
set(PARMETISTESTS pamgtest pamg_comm_repart_test)
endif(PARMETIS_FOUND)
-
-set(NORMALTESTS amgtest fastamg graphtest kamgtest pthreadamgtest)
+if(CMAKE_USE_PTHREADS_INIT)
+ set(PTHREADTESTS pthreadamgtest pthreadfastamgtest)
+endif(CMAKE_USE_PTHREADS_INIT)
+set(NORMALTESTS amgtest fastamg graphtest kamgtest ${PTHREADTESTS})
set(ALLTESTS ${MPITESTS} ${PARMETISTESTS} ${NORMALTESTS})
# We do not want want to build the tests during make all,
@@ -17,6 +19,12 @@ add_dependencies(${_test_target} ${ALLTESTS})
find_package(Threads)
if(CMAKE_USE_PTHREADS_INIT)
add_executable(pthreadamgtest pthreadamgtest.cc)
+ set_property(TARGET pthreadamgtest APPEND PROPERTY
+ COMPILE_DEFINITIONS "MYAMG=Dune::Amg::AMG<Operator,Vector,Smoother>")
+ target_link_libraries(pthreadamgtest "${CMAKE_THREAD_LIBS_INIT}")
+ add_executable(pthreadfastamgtest pthreadamgtest.cc)
+ set_property(TARGET pthreadfastamgtest APPEND PROPERTY
+ COMPILE_DEFINITIONS "MYAMG=Dune::Amg::FastAMG<Operator,Vector>")
target_link_libraries(pthreadamgtest "${CMAKE_THREAD_LIBS_INIT}")
endif(CMAKE_USE_PTHREADS_INIT)
diff --git a/dune/istl/paamg/test/pthreadamgtest.cc b/dune/istl/paamg/test/pthreadamgtest.cc
index be81e7cef..a4d3c1c20 100644
--- a/dune/istl/paamg/test/pthreadamgtest.cc
+++ b/dune/istl/paamg/test/pthreadamgtest.cc
@@ -20,12 +20,13 @@
#include <dune/common/parallel/indexset.hh>
#include <dune/common/parallel/collectivecommunication.hh>
#include <dune/istl/paamg/amg.hh>
+#include <dune/istl/paamg/fastamg.hh>
#include <dune/istl/paamg/pinfo.hh>
#include <dune/istl/solvers.hh>
#include <cstdlib>
#include <ctime>
#include <pthread.h>
-#define NUM_THREADS 1
+#define NUM_THREADS 3
typedef double XREAL;
@@ -68,11 +69,10 @@ 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;
+ MYAMG *amg;
Vector *b;
Vector *x;
Operator *fop;
@@ -99,6 +99,7 @@ void *solve(void* arg)
void *solve1(void* arg)
{
thread_arg *amgarg=(thread_arg*) arg;
+ *amgarg->x=0;
(*amgarg->amg).apply(*amgarg->x,*amgarg->b);
(*amgarg->amg).post(*amgarg->x);
@@ -107,12 +108,13 @@ void *solve1(void* arg)
void *solve2(void* arg)
{
thread_arg *amgarg=(thread_arg*) arg;
+ *amgarg->x=0;
(*amgarg->amg).pre(*amgarg->x,*amgarg->b);
(*amgarg->amg).apply(*amgarg->x,*amgarg->b);
(*amgarg->amg).post(*amgarg->x);
}
-template <int BS>
+template <int BS, typename AMG>
void testAMG(int N, int coarsenTarget, int ml)
{
@@ -177,8 +179,7 @@ void testAMG(int N, int coarsenTarget, int ml)
Dune::SeqScalarProduct<Vector> sp;
Smoother smoother(mat,1,1);
-
- AMG amg(fop, criterion, smootherArgs, 1, 1, 1, false);
+ AMG amg(fop, criterion);
double buildtime = watch.elapsed();
@@ -249,7 +250,8 @@ int main(int argc, char** argv)
if(argc>3)
ml = atoi(argv[3]);
- testAMG<1>(N, coarsenTarget, ml);
+ testAMG<1,MYAMG>(N, coarsenTarget, ml);
+
//testAMG<2>(N, coarsenTarget, ml);
}
--
GitLab