diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh
index d422b27317bda978867f8c33be6c27f7a1e39ef3..712078cf78cc79d0df0ae63a34e2cc5149e7e804 100644
--- a/dune/istl/paamg/amg.hh
+++ b/dune/istl/paamg/amg.hh
@@ -172,6 +172,11 @@ namespace Dune
const SmootherArgs& smootherArgs,
const ParallelInformation& pinfo=ParallelInformation());
+ /**
+ * @brief Copy constructor.
+ */
+ AMG(const AMG& amg);
+
~AMG();
/** \copydoc Preconditioner::pre */
@@ -291,13 +296,33 @@ namespace Dune
std::size_t verbosity_;
};
+ template<class M, class X, class S, class PI, class A>
+ inline AMG<M,X,S,PI,A>::AMG(const AMG& amg)
+ : matrices_(amg.matrices_), smootherArgs_(amg.smootherArgs_),
+ smoothers_(amg.smoothers_), solver_(amg.solver_),
+ rhs_(), lhs_(), update_(),
+ scalarProduct_(amg.scalarProduct_), gamma_(amg.gamma_),
+ preSteps_(amg.preSteps_), postSteps_(amg.postSteps_),
+ level(amg.level), buildHierarchy_(amg.buildHierarchy_),
+ additive(amg.additive), 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.update_)
+ update_=new Hierarchy<Domain,A>(*amg.update_);
+ }
+
template<class M, class X, class S, class PI, class A>
AMG<M,X,S,PI,A>::AMG(const OperatorHierarchy& matrices, CoarseSolver& coarseSolver,
const SmootherArgs& smootherArgs,
std::size_t gamma, std::size_t preSmoothingSteps,
std::size_t postSmoothingSteps, bool additive_)
: matrices_(&matrices), smootherArgs_(smootherArgs),
- smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver), scalarProduct_(0),
+ smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
+ rhs_(), lhs_(), update_(), scalarProduct_(0),
gamma_(gamma), preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(false),
additive(additive_), coarsesolverconverged(true),
coarseSmoother_(), verbosity_(2)
@@ -313,7 +338,8 @@ namespace Dune
const SmootherArgs& smootherArgs,
const Parameters& parms)
: matrices_(&matrices), smootherArgs_(smootherArgs),
- smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver), scalarProduct_(0),
+ smoothers_(new Hierarchy<Smoother,A>), solver_(&coarseSolver),
+ rhs_(), lhs_(), update_(), scalarProduct_(0),
gamma_(parms.getGamma()), preSteps_(parms.getNoPreSmoothSteps()),
postSteps_(parms.getNoPostSmoothSteps()), buildHierarchy_(false),
additive(parms.getAdditive()), coarsesolverconverged(true),
@@ -335,7 +361,8 @@ namespace Dune
bool additive_,
const PI& pinfo)
: smootherArgs_(smootherArgs),
- smoothers_(new Hierarchy<Smoother,A>), solver_(), scalarProduct_(0), gamma_(gamma),
+ smoothers_(new Hierarchy<Smoother,A>), solver_(),
+ rhs_(), lhs_(), update_(), scalarProduct_(0), gamma_(gamma),
preSteps_(preSmoothingSteps), postSteps_(postSmoothingSteps), buildHierarchy_(true),
additive(additive_), coarsesolverconverged(true),
coarseSmoother_(), verbosity_(criterion.debugLevel())
@@ -355,7 +382,8 @@ namespace Dune
const SmootherArgs& smootherArgs,
const PI& pinfo)
: smootherArgs_(smootherArgs),
- smoothers_(new Hierarchy<Smoother,A>), solver_(), scalarProduct_(0),
+ smoothers_(new Hierarchy<Smoother,A>), solver_(),
+ rhs_(), lhs_(), update_(), scalarProduct_(0),
gamma_(criterion.getGamma()), preSteps_(criterion.getNoPreSmoothSteps()),
postSteps_(criterion.getNoPostSmoothSteps()), buildHierarchy_(true),
additive(criterion.getAdditive()), coarsesolverconverged(true),
@@ -378,6 +406,15 @@ namespace Dune
if(coarseSmoother_)
coarseSmoother_.reset();
}
+ if(lhs_)
+ delete lhs_;
+ lhs_=nullptr;
+ if(update_)
+ delete update_;
+ update_=nullptr;
+ if(rhs_)
+ delete rhs_;
+ rhs_=nullptr;
}
template<class M, class X, class S, class PI, class A>
@@ -492,11 +529,17 @@ 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 Hierarchy<Range,A>(*copy);
+ if(rhs_)
+ delete rhs_;
+ rhs_ = new Hierarchy<Range,A>(copy);
Domain* dcopy = new Domain(x);
- lhs_ = new Hierarchy<Domain,A>(*dcopy);
+ if(lhs_)
+ delete lhs_;
+ lhs_ = new Hierarchy<Domain,A>(dcopy);
dcopy = new Domain(x);
- update_ = new Hierarchy<Domain,A>(*dcopy);
+ if(update_)
+ delete update_;
+ update_ = new Hierarchy<Domain,A>(dcopy);
matrices_->coarsenVector(*rhs_);
matrices_->coarsenVector(*lhs_);
matrices_->coarsenVector(*update_);
@@ -826,13 +869,13 @@ namespace Dune
smoother->post(*lhs);
smoother->post(*lhs);
}
- delete &(*lhs_->finest());
+ //delete &(*lhs_->finest());
delete lhs_;
lhs_=nullptr;
- delete &(*update_->finest());
+ //delete &(*update_->finest());
delete update_;
update_=nullptr;
- delete &(*rhs_->finest());
+ //delete &(*rhs_->finest());
delete rhs_;
rhs_=nullptr;
}
diff --git a/dune/istl/paamg/hierarchy.hh b/dune/istl/paamg/hierarchy.hh
index dbf9ae32ee94a396f19094fb3eed7c592717a234..0036300fd13d8172d33c0536039b0cc10f2d67d3 100644
--- a/dune/istl/paamg/hierarchy.hh
+++ b/dune/istl/paamg/hierarchy.hh
@@ -117,11 +117,23 @@ namespace Dune
*/
Hierarchy(MemberType& first);
+ /**
+ * @brief Construct a new hierarchy.
+ * @param first Pointer to the first element in the hierarchy.
+ * @warning Hierarchy will be responsible for the memory
+ * management of the pointer.
+ */
+ Hierarchy(MemberType* first);
+
/**
* @brief Construct a new empty hierarchy.
*/
Hierarchy();
+ /**
+ * @brief Copy constructor.
+ */
+ Hierarchy(const Hierarchy& other);
/**
* @brief Add an element on a coarser level.
* @param args The arguments needed for the construction.
@@ -234,7 +246,7 @@ namespace Dune
void deleteRedistributed()
{
- element_->redistributed_ = 0;
+ element_->redistributed_ = nullptr;
}
private:
@@ -282,9 +294,6 @@ 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. */
@@ -1233,13 +1242,25 @@ namespace Dune
{
finest_ = allocator_.allocate(1,0);
finest_->element_ = &first;
- finest_->redistributed_ = 0;
+ finest_->redistributed_ = nullptr;
nonAllocated_ = finest_;
coarsest_ = finest_;
- coarsest_->coarser_ = coarsest_->finer_ = 0;
+ coarsest_->coarser_ = coarsest_->finer_ = nullptr;
levels_ = 1;
}
+ template<class T, class A>
+ Hierarchy<T,A>::Hierarchy(MemberType* first)
+ : allocator_()
+ {
+ finest_ = allocator_.allocate(1,0);
+ finest_->element_ = first;
+ finest_->redistributed_ = nullptr;
+ nonAllocated_ = nullptr;
+ coarsest_ = finest_;
+ coarsest_->coarser_ = coarsest_->finer_ = nullptr;
+ levels_ = 1;
+ }
template<class T, class A>
Hierarchy<T,A>::~Hierarchy()
{
@@ -1252,8 +1273,45 @@ namespace Dune
ConstructionTraits<T>::deconstruct(current->element_);
}
allocator_.deallocate(current, 1);
- //coarsest_->coarser_ = 0;
+ current=nullptr;
+ //coarsest_->coarser_ = nullptr;
+ }
+ }
+
+ template<class T, class A>
+ Hierarchy<T,A>::Hierarchy(const Hierarchy& other)
+ : nonAllocated_(), allocator_(other.allocator_),
+ levels_(other.levels_)
+ {
+ if(!other.finest_)
+ {
+ finest_=coarsest_=nonAllocated_=nullptr;
+ return;
+ }
+ finest_=allocator_.allocate(1,0);
+ Element* finer_ = nullptr;
+ Element* current_ = finest_;
+ Element* otherCurrent_ = other.finest_;
+
+ while(otherCurrent_)
+ {
+ T* t=new T(*(otherCurrent_->element_));
+ current_->element_=t;
+ current_->finer_=finer_;
+ if(otherCurrent_->redistributed_)
+ current_->redistributed_ = new T(*otherCurrent_->redistributed_);
+ else
+ current_->redistributed_= nullptr;
+ finer_=current_;
+ if(otherCurrent_->coarser_)
+ {
+ current_->coarser_=allocator_.allocate(1,0);
+ current_=current_->coarser_;
+ }else
+ current_->coarser_=nullptr;
+ otherCurrent_=otherCurrent_->coarser_;
}
+ coarsest_=current_;
}
template<class T, class A>
@@ -1276,15 +1334,15 @@ namespace Dune
coarsest_ = allocator_.allocate(1,0);
coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
finest_ = coarsest_;
- coarsest_->finer_ = 0;
+ coarsest_->finer_ = nullptr;
}else{
coarsest_->coarser_ = allocator_.allocate(1,0);
coarsest_->coarser_->finer_ = coarsest_;
coarsest_ = coarsest_->coarser_;
coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
}
- coarsest_->redistributed_ = 0;
- coarsest_->coarser_=0;
+ coarsest_->redistributed_ = nullptr;
+ coarsest_->coarser_=nullptr;
++levels_;
}
@@ -1297,12 +1355,12 @@ namespace Dune
finest_ = allocator_.allocate(1,0);
finest_->element = ConstructionTraits<T>::construct(args);
coarsest_ = finest_;
- coarsest_->coarser_ = coarsest_->finer_ = 0;
+ coarsest_->coarser_ = coarsest_->finer_ = nullptr;
}else{
finest_->finer_ = allocator_.allocate(1,0);
finest_->finer_->coarser_ = finest_;
finest_ = finest_->finer_;
- finest_->finer = 0;
+ finest_->finer = nullptr;
finest_->element = ConstructionTraits<T>::construct(args);
}
++levels_;
diff --git a/dune/istl/paamg/test/pthreadamgtest.cc b/dune/istl/paamg/test/pthreadamgtest.cc
index f43c9e106cfe8be119c79fd8513571d030b7b707..be81e7cef2283cf14bcbe9f75da95e4b3275e9eb 100644
--- a/dune/istl/paamg/test/pthreadamgtest.cc
+++ b/dune/istl/paamg/test/pthreadamgtest.cc
@@ -25,7 +25,7 @@
#include <cstdlib>
#include <ctime>
#include <pthread.h>
-#define NUM_THREADS 5
+#define NUM_THREADS 1
typedef double XREAL;
@@ -96,6 +96,21 @@ void *solve(void* arg)
pthread_exit(NULL);
}
+void *solve1(void* arg)
+{
+ thread_arg *amgarg=(thread_arg*) arg;
+ (*amgarg->amg).apply(*amgarg->x,*amgarg->b);
+ (*amgarg->amg).post(*amgarg->x);
+
+}
+
+void *solve2(void* arg)
+{
+ thread_arg *amgarg=(thread_arg*) arg;
+ (*amgarg->amg).pre(*amgarg->x,*amgarg->b);
+ (*amgarg->amg).apply(*amgarg->x,*amgarg->b);
+ (*amgarg->amg).post(*amgarg->x);
+}
template <int BS>
void testAMG(int N, int coarsenTarget, int ml)
@@ -184,6 +199,37 @@ void testAMG(int N, int coarsenTarget, int ml)
for(int i=0; i < NUM_THREADS; ++i)
pthread_join(threads[i], &retval);
+
+ amgs.clear();
+ args.clear();
+ amg.pre(x, b);
+ amgs.resize(NUM_THREADS, amg);
+ 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, solve1, (void*) &args[i]);
+ }
+ for(int i=0; i < NUM_THREADS; ++i)
+ pthread_join(threads[i], &retval);
+
+ amgs.clear();
+ args.clear();
+ amg.pre(x, b);
+ amgs.resize(NUM_THREADS, amg);
+ 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, solve2, (void*) &args[i]);
+ }
+ for(int i=0; i < NUM_THREADS; ++i)
+ pthread_join(threads[i], &retval);
+
}