From bd05ab5a26142c1bd7b9e0b2050dc6992714cb3e Mon Sep 17 00:00:00 2001
From: Christian Engwer <christi@dune-project.org>
Date: Thu, 24 Jan 2019 15:42:05 +0100
Subject: [PATCH] [amg] update Hierarchy to internally manage the finest level
 via shared_ptr

this is a first step to cleanup memory management in the AMG and it will
facilitate the construction of an AMG using a factory.

* 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
---
 dune/istl/paamg/amg.hh                | 24 ++++----
 dune/istl/paamg/fastamg.hh            | 24 ++++----
 dune/istl/paamg/hierarchy.hh          | 81 +++++++++++++++++++--------
 dune/istl/paamg/test/hierarchytest.cc |  7 ++-
 4 files changed, 88 insertions(+), 48 deletions(-)

diff --git a/dune/istl/paamg/amg.hh b/dune/istl/paamg/amg.hh
index 249c0e093..12d98420c 100644
--- a/dune/istl/paamg/amg.hh
+++ b/dune/istl/paamg/amg.hh
@@ -181,7 +181,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.
@@ -356,7 +357,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);
     }
 
 
@@ -447,11 +449,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);
 
@@ -600,18 +605,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_);
diff --git a/dune/istl/paamg/fastamg.hh b/dune/istl/paamg/fastamg.hh
index 5b712656d..a89b7e6d0 100644
--- a/dune/istl/paamg/fastamg.hh
+++ b/dune/istl/paamg/fastamg.hh
@@ -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_);
diff --git a/dune/istl/paamg/hierarchy.hh b/dune/istl/paamg/hierarchy.hh
index ff59ff0e5..0499e0f18 100644
--- a/dune/istl/paamg/hierarchy.hh
+++ b/dune/istl/paamg/hierarchy.hh
@@ -104,19 +104,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.
@@ -287,12 +295,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. */
@@ -350,9 +362,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();
 
@@ -649,12 +660,12 @@ 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)
     {
-      if (SolverCategory::category(fineOperator) != SolverCategory::category(pinfo))
+      if (SolverCategory::category(*fineMatrix) != SolverCategory::category(*pinfo))
         DUNE_THROW(ISTLError, "MatrixOperator and ParallelInformation must belong to the same category!");
     }
 
@@ -1226,17 +1237,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;
@@ -1244,12 +1266,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;
@@ -1260,7 +1281,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_);
@@ -1273,12 +1296,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);
@@ -1323,9 +1346,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{
@@ -1344,9 +1372,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{
diff --git a/dune/istl/paamg/test/hierarchytest.cc b/dune/istl/paamg/test/hierarchytest.cc
index d7de15259..9e1619dc0 100644
--- a/dune/istl/paamg/test/hierarchytest.cc
+++ b/dune/istl/paamg/test/hierarchytest.cc
@@ -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;
-- 
GitLab