Skip to content
Snippets Groups Projects
Commit c0ccfcbd authored by Markus Blatt's avatar Markus Blatt
Browse files

Allow copying AMG even after calling pre()

parent 51445464
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/istl/paamg/amg.hh
+ 53
− 10
View file @ c0ccfcbd
......@@ -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;
}
......
dune/istl/paamg/hierarchy.hh
+ 70
− 12
View file @ c0ccfcbd
......@@ -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_;
......
dune/istl/paamg/test/pthreadamgtest.cc
+ 47
− 1
View file @ c0ccfcbd
......@@ -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);
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment