Skip to content
Snippets Groups Projects
Commit 13fb32f7 authored by Christian Engwer's avatar Christian Engwer
Browse files

[amg] update Hierarchy to internally manage the finest level vis shared_ptr

this change became necessary to correctly manage data created in the
solver factory. We implement the following changes:

* we introduce a new constructor taking a shared_ptr<MemberType>
* all old constructors are deprecated
* we update the test and the AMG implementations to the shared_ptr and
  pass it to Hierarchy
parent 41c7002a
No related branches found
No related tags found
No related merge requests found
Pipeline #
......@@ -208,7 +208,8 @@ namespace Dune
* @param pinfo The fine level parallel information.
*/
template<class C>
void createHierarchies(C& criterion, Operator& matrix,
void createHierarchies(C& criterion,
const std::shared_ptr<const Operator>& matrixptr,
const PI& pinfo);
/**
* @brief A struct that holds the context of the current level.
......@@ -383,12 +384,12 @@ namespace Dune
// TODO: reestablish compile time checks.
//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);
auto matrixptr = stackobject_to_shared_ptr(matrix);
createHierarchies(criterion, matrixptr, pinfo);
}
template<class M, class X, class S, class PI, class A>
AMG<M,X,S,PI,A>::AMG(std::shared_ptr<const Operator> matrix,
AMG<M,X,S,PI,A>::AMG(std::shared_ptr<const Operator> matrixptr,
const ParameterTree& configuration,
const ParallelInformation& pinfo) :
smoothers_(new Hierarchy<Smoother,A>),
......@@ -438,7 +439,7 @@ namespace Dune
verbosity_ = configuration.get("verbosity",2);
criterion.setDebugLevel (verbosity_);
createHierarchies(criterion, const_cast<Operator&>(*matrix), pinfo);
createHierarchies(criterion, matrixptr, pinfo);
}
template<class M, class X, class S, class PI, class A>
......@@ -528,11 +529,14 @@ namespace Dune
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)
void AMG<M,X,S,PI,A>::createHierarchies(C& criterion,
const std::shared_ptr<const Operator>& matrixptr,
const PI& pinfo)
{
Timer watch;
matrices_.reset(new OperatorHierarchy(matrix, pinfo));
matrices_ = std::make_shared<OperatorHierarchy>(
std::const_pointer_cast<Operator>(matrixptr),
stackobject_to_shared_ptr(const_cast<PI&>(pinfo)));
matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
......@@ -648,18 +652,15 @@ namespace Dune
else
// No smoother to make x consistent! Do it by hand
matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
Range* copy = new Range(b);
if(rhs_)
delete rhs_;
rhs_ = new Hierarchy<Range,A>(copy);
Domain* dcopy = new Domain(x);
rhs_ = new Hierarchy<Range,A>(std::make_shared<Range>(b));
if(lhs_)
delete lhs_;
lhs_ = new Hierarchy<Domain,A>(dcopy);
dcopy = new Domain(x);
lhs_ = new Hierarchy<Domain,A>(std::make_shared<Domain>(x));
if(update_)
delete update_;
update_ = new Hierarchy<Domain,A>(dcopy);
update_ = new Hierarchy<Domain,A>(std::make_shared<Domain>(x));
matrices_->coarsenVector(*rhs_);
matrices_->coarsenVector(*lhs_);
matrices_->coarsenVector(*update_);
......
......@@ -167,7 +167,8 @@ namespace Dune
* @param pinfo The fine level parallel information.
*/
template<class C>
void createHierarchies(C& criterion, Operator& matrix,
void createHierarchies(C& criterion,
const std::shared_ptr<const Operator>& matrixptr,
const PI& pinfo);
/**
......@@ -342,7 +343,8 @@ namespace Dune
// TODO: reestablish compile time checks.
//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);
auto matrixptr = stackobject_to_shared_ptr(matrix);
createHierarchies(criterion, matrixptr, pinfo);
}
template<class M, class X, class PI, class A>
......@@ -367,11 +369,14 @@ namespace Dune
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)
void FastAMG<M,X,PI,A>::createHierarchies(C& criterion,
const std::shared_ptr<const Operator>& matrixptr,
const PI& pinfo)
{
Timer watch;
matrices_.reset(new OperatorHierarchy(matrix, pinfo));
matrices_ = std::make_shared<OperatorHierarchy>(
std::const_pointer_cast<Operator>(matrixptr),
stackobject_to_shared_ptr(const_cast<PI&>(pinfo)));
matrices_->template build<NegateSet<typename PI::OwnerSet> >(criterion);
......@@ -483,16 +488,13 @@ namespace Dune
watch1.reset();
// No smoother to make x consistent! Do it by hand
matrices_->parallelInformation().coarsest()->copyOwnerToAll(x,x);
Range* copy = new Range(b);
if(rhs_)
delete rhs_;
rhs_ = new Hierarchy<Range,A>(copy);
Domain* dcopy = new Domain(x);
rhs_ = new Hierarchy<Range,A>(std::make_shared<Range>(b));
if(lhs_)
delete lhs_;
lhs_ = new Hierarchy<Domain,A>(dcopy);
dcopy = new Domain(x);
residual_ = new Hierarchy<Domain,A>(dcopy);
lhs_ = new Hierarchy<Domain,A>(std::make_shared<Domain>(x));
residual_ = new Hierarchy<Domain,A>(std::make_shared<Domain>(x));
matrices_->coarsenVector(*rhs_);
matrices_->coarsenVector(*lhs_);
matrices_->coarsenVector(*residual_);
......
......@@ -111,19 +111,27 @@ namespace Dune
typedef typename ConstructionTraits<T>::Arguments Arguments;
/**
* @brief Construct a new hierarchy.
* @param first std::shared_ptr to the first element in the hierarchy.
*/
Hierarchy(const std::shared_ptr<MemberType> & first);
/**
* @brief Construct a new hierarchy.
* @param first The first element in the hierarchy.
* @todo remove in 2.6
*/
Hierarchy(MemberType& first);
Hierarchy(MemberType& first) DUNE_DEPRECATED_MSG("Use the constructor taking a shared_ptr");
/**
* @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.
* @todo remove in 2.6
*/
Hierarchy(MemberType* first);
Hierarchy(MemberType* first) DUNE_DEPRECATED_MSG("Use the constructor taking a shared_ptr");
/**
* @brief Construct a new empty hierarchy.
......@@ -294,12 +302,16 @@ namespace Dune
~Hierarchy();
private:
/** @brief fix memory management of the finest element in the hierarchy
This object is passed in from outside, so that we have to
deal with shared ownership.
*/
std::shared_ptr<MemberType> originalFinest_;
/** @brief The finest element in the hierarchy. */
Element* finest_;
/** @brief The coarsest element in the hierarchy. */
Element* coarsest_;
/** @brief Whether the first element was not allocated by us. */
Element* nonAllocated_;
/** @brief The allocator for the list elements. */
Allocator allocator_;
/** @brief The number of levels in the hierarchy. */
......@@ -357,9 +369,8 @@ namespace Dune
* @param fineMatrix The matrix to coarsen.
* @param pinfo The information about the parallel data decomposition at the first level.
*/
MatrixHierarchy(const MatrixOperator& fineMatrix,
const ParallelInformation& pinfo=ParallelInformation());
MatrixHierarchy(std::shared_ptr<MatrixOperator> fineMatrix,
std::shared_ptr<ParallelInformation> pinfo = std::make_shared<ParallelInformation>());
~MatrixHierarchy();
......@@ -657,10 +668,10 @@ namespace Dune
}
template<class M, class IS, class A>
MatrixHierarchy<M,IS,A>::MatrixHierarchy(const MatrixOperator& fineOperator,
const ParallelInformation& pinfo)
: matrices_(const_cast<MatrixOperator&>(fineOperator)),
parallelInformation_(const_cast<ParallelInformation&>(pinfo))
MatrixHierarchy<M,IS,A>::MatrixHierarchy(std::shared_ptr<MatrixOperator> fineMatrix,
std::shared_ptr<ParallelInformation> pinfo)
: matrices_(fineMatrix),
parallelInformation_(pinfo)
{
#warning reenable category checks
// static_assert((static_cast<int>(MatrixOperator::category) ==
......@@ -1242,17 +1253,28 @@ namespace Dune
template<class T, class A>
Hierarchy<T,A>::Hierarchy()
: finest_(0), coarsest_(0), nonAllocated_(0), allocator_(), levels_(0)
: finest_(0), coarsest_(0), allocator_(), levels_(0)
{}
template<class T, class A>
Hierarchy<T,A>::Hierarchy(const std::shared_ptr<MemberType> & first)
: originalFinest_(first), allocator_()
{
finest_ = allocator_.allocate(1,0);
finest_->element_ = originalFinest_.get();
finest_->redistributed_ = nullptr;
coarsest_ = finest_;
coarsest_->coarser_ = coarsest_->finer_ = nullptr;
levels_ = 1;
}
template<class T, class A>
Hierarchy<T,A>::Hierarchy(MemberType& first)
: allocator_()
: originalFinest_(stackobject_to_shared_ptr(first)), allocator_()
{
finest_ = allocator_.allocate(1,0);
finest_->element_ = &first;
finest_->element_ = originalFinest_.get();
finest_->redistributed_ = nullptr;
nonAllocated_ = finest_;
coarsest_ = finest_;
coarsest_->coarser_ = coarsest_->finer_ = nullptr;
levels_ = 1;
......@@ -1260,12 +1282,11 @@ namespace Dune
template<class T, class A>
Hierarchy<T,A>::Hierarchy(MemberType* first)
: allocator_()
: originalFinest_(first), allocator_()
{
finest_ = allocator_.allocate(1,0);
finest_->element_ = first;
finest_->element_ = originalFinest_.get();
finest_->redistributed_ = nullptr;
nonAllocated_ = nullptr;
coarsest_ = finest_;
coarsest_->coarser_ = coarsest_->finer_ = nullptr;
levels_ = 1;
......@@ -1276,7 +1297,9 @@ namespace Dune
while(coarsest_) {
Element* current = coarsest_;
coarsest_ = coarsest_->finer_;
if(current != nonAllocated_) {
// we changed the internal behaviour
// now the finest level is _always_ managedd by a shared_ptr
if(current != finest_) {
if(current->redistributed_)
ConstructionTraits<T>::deconstruct(current->redistributed_);
ConstructionTraits<T>::deconstruct(current->element_);
......@@ -1289,12 +1312,12 @@ namespace Dune
template<class T, class A>
Hierarchy<T,A>::Hierarchy(const Hierarchy& other)
: nonAllocated_(), allocator_(other.allocator_),
: allocator_(other.allocator_),
levels_(other.levels_)
{
if(!other.finest_)
{
finest_=coarsest_=nonAllocated_=nullptr;
finest_=coarsest_=nullptr;
return;
}
finest_=allocator_.allocate(1,0);
......@@ -1339,9 +1362,14 @@ namespace Dune
void Hierarchy<T,A>::addCoarser(Arguments& args)
{
if(!coarsest_) {
// we have no levels at all...
assert(!finest_);
// allocate into the shared_ptr
originalFinest_ =
std::shared_ptr<MemberType>(
ConstructionTraits<MemberType>::construct(args));
coarsest_ = allocator_.allocate(1,0);
coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
coarsest_->element_ = originalFinest_.get();
finest_ = coarsest_;
coarsest_->finer_ = nullptr;
}else{
......@@ -1360,9 +1388,14 @@ namespace Dune
void Hierarchy<T,A>::addFiner(Arguments& args)
{
if(!finest_) {
// we have no levels at all...
assert(!coarsest_);
// allocate into the shared_ptr
originalFinest_ =
std::shared_ptr<MemberType>(
ConstructionTraits<MemberType>::construct(args));
finest_ = allocator_.allocate(1,0);
finest_->element = ConstructionTraits<T>::construct(args);
finest_->element = originalFinest_.get();
coarsest_ = finest_;
coarsest_->coarser_ = coarsest_->finer_ = nullptr;
}else{
......
......@@ -43,8 +43,11 @@ void testHierarchy(int N)
typedef Dune::Amg::Hierarchy<Vector> VHierarchy;
Operator op(mat, pinfo);
Hierarchy hierarchy(op, pinfo);
VHierarchy vh(b);
Hierarchy hierarchy(
Dune::stackobject_to_shared_ptr(op),
Dune::stackobject_to_shared_ptr(pinfo));
VHierarchy vh(
Dune::stackobject_to_shared_ptr(b));
typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<BCRSMat,Dune::Amg::FirstDiagonal> >
Criterion;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment