1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_AMGHIERARCHY_HH
4 #define DUNE_AMGHIERARCHY_HH
5 
6 #include <list>
7 #include <memory>
8 #include <limits>
9 #include <dune/common/stdstreams.hh>
10 #include <dune/common/timer.hh>
11 #include <dune/common/bigunsignedint.hh>
12 #include <dune/istl/paamg/construction.hh>
13 
14 namespace Dune
15 {
16   namespace Amg
17   {
18     /**
19      * @addtogroup ISTL_PAAMG
20      *
21      * @{
22      */
23 
24     /** @file
25      * @author Markus Blatt
26      * @brief Provides a classes representing the hierarchies in AMG.
27      */
28 
29     /**
30      * @brief A hierarchy of containers (e.g. matrices or vectors)
31      *
32      * Because sometimes a redistribution of the parallel data might be
33      * advisable one can add redistributed version of the container at
34      * each level.
35      */
36     template<typename T, typename A=std::allocator<T> >
37     class Hierarchy
38     {
39     public:
40       /**
41        * @brief The type of the container we store.
42        */
43       typedef T MemberType;
44 
45       template<typename T1, typename T2>
46       class LevelIterator;
47 
48     private:
49       /**
50        * @brief An element in the hierarchy.
51        */
52       struct Element
53       {
54         friend class LevelIterator<Hierarchy<T,A>, T>;
55         friend class LevelIterator<const Hierarchy<T,A>, const T>;
56 
57         /** @brief The next coarser element in the list. */
58         std::weak_ptr<Element> coarser_;
59 
60         /** @brief The next finer element in the list. */
61         std::shared_ptr<Element> finer_;
62 
63         /** @brief Pointer to the element. */
64         std::shared_ptr<MemberType> element_;
65 
66         /** @brief The redistributed version of the element. */
67         std::shared_ptr<MemberType> redistributed_;
68       };
69     public:
70 
71       /**
72        * @brief The allocator to use for the list elements.
73        */
74       using Allocator = typename std::allocator_traits<A>::template rebind_alloc<Element>;
75 
76       typedef typename ConstructionTraits<T>::Arguments Arguments;
77 
78       /**
79        * @brief Construct a new hierarchy.
80        * @param first std::shared_ptr to the first element in the hierarchy.
81        */
82       Hierarchy(const std::shared_ptr<MemberType> & first);
83 
84       /**
85        * @brief Construct an empty hierarchy.
86        */
Hierarchy()87       Hierarchy() : levels_(0)
88       {}
89 
90       /**
91        * @brief Copy constructor (deep copy!).
92        */
93       Hierarchy(const Hierarchy& other);
94 
95       /**
96        * @brief Add an element on a coarser level.
97        * @param args The arguments needed for the construction.
98        */
99       void addCoarser(Arguments& args);
100 
101       void addRedistributedOnCoarsest(Arguments& args);
102 
103       /**
104        * @brief Add an element on a finer level.
105        * @param args The arguments needed for the construction.
106        */
107       void addFiner(Arguments& args);
108 
109       /**
110        * @brief Iterator over the levels in the hierarchy.
111        *
112        * operator++() moves to the next coarser level in the hierarchy.
113        * while operator--() moves to the next finer level in the hierarchy.
114        */
115       template<class C, class T1>
116       class LevelIterator
117         : public BidirectionalIteratorFacade<LevelIterator<C,T1>,T1,T1&>
118       {
119         friend class LevelIterator<typename std::remove_const<C>::type,
120             typename std::remove_const<T1>::type >;
121         friend class LevelIterator<const typename std::remove_const<C>::type,
122             const typename std::remove_const<T1>::type >;
123 
124       public:
125         /** @brief Constructor. */
LevelIterator()126         LevelIterator()
127         {}
128 
LevelIterator(std::shared_ptr<Element> element)129         LevelIterator(std::shared_ptr<Element> element)
130           : element_(element)
131         {}
132 
133         /** @brief Copy constructor. */
LevelIterator(const LevelIterator<typename std::remove_const<C>::type,typename std::remove_const<T1>::type> & other)134         LevelIterator(const LevelIterator<typename std::remove_const<C>::type,
135                           typename std::remove_const<T1>::type>& other)
136           : element_(other.element_)
137         {}
138 
139         /** @brief Copy constructor. */
LevelIterator(const LevelIterator<const typename std::remove_const<C>::type,const typename std::remove_const<T1>::type> & other)140         LevelIterator(const LevelIterator<const typename std::remove_const<C>::type,
141                           const typename std::remove_const<T1>::type>& other)
142           : element_(other.element_)
143         {}
144 
145         /**
146          * @brief Equality check.
147          */
equals(const LevelIterator<typename std::remove_const<C>::type,typename std::remove_const<T1>::type> & other) const148         bool equals(const LevelIterator<typename std::remove_const<C>::type,
149                         typename std::remove_const<T1>::type>& other) const
150         {
151           return element_ == other.element_;
152         }
153 
154         /**
155          * @brief Equality check.
156          */
equals(const LevelIterator<const typename std::remove_const<C>::type,const typename std::remove_const<T1>::type> & other) const157         bool equals(const LevelIterator<const typename std::remove_const<C>::type,
158                         const typename std::remove_const<T1>::type>& other) const
159         {
160           return element_ == other.element_;
161         }
162 
163         /** @brief Dereference the iterator. */
dereference() const164         T1& dereference() const
165         {
166           return *(element_->element_);
167         }
168 
169         /** @brief Move to the next coarser level */
increment()170         void increment()
171         {
172           element_ = element_->coarser_.lock();
173         }
174 
175         /** @brief Move to the next fine level */
decrement()176         void decrement()
177         {
178           element_ = element_->finer_;
179         }
180 
181         /**
182          * @brief Check whether there was a redistribution at the current level.
183          * @return True if there is a redistributed version of the container at the current level.
184          */
isRedistributed() const185         bool isRedistributed() const
186         {
187           return (bool)element_->redistributed_;
188         }
189 
190         /**
191          * @brief Get the redistributed container.
192          * @return The redistributed container.
193          */
getRedistributed() const194         T1& getRedistributed() const
195         {
196           assert(element_->redistributed_);
197           return *element_->redistributed_;
198         }
addRedistributed(std::shared_ptr<T1> t)199         void addRedistributed(std::shared_ptr<T1> t)
200         {
201           element_->redistributed_ = t;
202         }
203 
deleteRedistributed()204         void deleteRedistributed()
205         {
206           element_->redistributed_ = nullptr;
207         }
208 
209       private:
210         std::shared_ptr<Element> element_;
211       };
212 
213       /** @brief Type of the mutable iterator. */
214       typedef LevelIterator<Hierarchy<T,A>,T> Iterator;
215 
216       /** @brief Type of the const iterator. */
217       typedef LevelIterator<const Hierarchy<T,A>, const T> ConstIterator;
218 
219       /**
220        * @brief Get an iterator positioned at the finest level.
221        * @return An iterator positioned at the finest level.
222        */
223       Iterator finest();
224 
225       /**
226        * @brief Get an iterator positioned at the coarsest level.
227        * @return An iterator positioned at the coarsest level.
228        */
229       Iterator coarsest();
230 
231 
232       /**
233        * @brief Get an iterator positioned at the finest level.
234        * @return An iterator positioned at the finest level.
235        */
236       ConstIterator finest() const;
237 
238       /**
239        * @brief Get an iterator positioned at the coarsest level.
240        * @return An iterator positioned at the coarsest level.
241        */
242       ConstIterator coarsest() const;
243 
244       /**
245        * @brief Get the number of levels in the hierarchy.
246        * @return The number of levels.
247        */
248       std::size_t levels() const;
249 
250     private:
251       /** @brief fix memory management of the finest element in the hierarchy
252 
253           This object is passed in from outside, so that we have to
254           deal with shared ownership.
255       */
256       std::shared_ptr<MemberType> originalFinest_;
257       /** @brief The finest element in the hierarchy. */
258       std::shared_ptr<Element> finest_;
259       /** @brief The coarsest element in the hierarchy. */
260       std::shared_ptr<Element> coarsest_;
261       /** @brief The allocator for the list elements. */
262       Allocator allocator_;
263       /** @brief The number of levels in the hierarchy. */
264       int levels_;
265     };
266 
267     template<class T, class A>
Hierarchy(const std::shared_ptr<MemberType> & first)268     Hierarchy<T,A>::Hierarchy(const std::shared_ptr<MemberType> & first)
269       : originalFinest_(first)
270     {
271       finest_ = std::allocate_shared<Element>(allocator_);
272       finest_->element_ = originalFinest_;
273       coarsest_ = finest_;
274       levels_ = 1;
275     }
276 
277     //! \brief deep copy of a given hierarchy
278     //TODO: do we actually want to support this? This might be very expensive?!
279     template<class T, class A>
Hierarchy(const Hierarchy & other)280     Hierarchy<T,A>::Hierarchy(const Hierarchy& other)
281     : allocator_(other.allocator_),
282       levels_(other.levels_)
283     {
284       if(!other.finest_)
285       {
286         finest_=coarsest_=nullptr;
287         return;
288       }
289       finest_ = std::allocate_shared<Element>(allocator_);
290       std::shared_ptr<Element> finer_;
291       std::shared_ptr<Element> current_ = finest_;
292       std::weak_ptr<Element> otherWeak_ = other.finest_;
293 
294       while(! otherWeak_.expired())
295       {
296         // create shared_ptr from weak_ptr, we just checked that this is safe
297         std::shared_ptr<Element> otherCurrent_ = std::shared_ptr<Element>(otherWeak_);
298         // clone current level
299         //TODO: should we use the allocator?
300         current_->element_ =
301           std::make_shared<MemberType>(*(otherCurrent_->element_));
302         current_->finer_=finer_;
303         if(otherCurrent_->redistributed_)
304           current_->redistributed_ =
305             std::make_shared<MemberType>(*(otherCurrent_->redistributed_));
306         finer_=current_;
307         if(not otherCurrent_->coarser_.expired())
308         {
309           auto c = std::allocate_shared<Element>(allocator_);
310           current_->coarser_ = c;
311           current_ = c;
312         }
313         // go to coarser level
314         otherWeak_ = otherCurrent_->coarser_;
315       }
316       coarsest_=current_;
317     }
318 
319     template<class T, class A>
levels() const320     std::size_t Hierarchy<T,A>::levels() const
321     {
322       return levels_;
323     }
324 
325     template<class T, class A>
addRedistributedOnCoarsest(Arguments & args)326     void Hierarchy<T,A>::addRedistributedOnCoarsest(Arguments& args)
327     {
328       coarsest_->redistributed_ = ConstructionTraits<MemberType>::construct(args);
329     }
330 
331     template<class T, class A>
addCoarser(Arguments & args)332     void Hierarchy<T,A>::addCoarser(Arguments& args)
333     {
334       if(!coarsest_) {
335         // we have no levels at all...
336         assert(!finest_);
337         // allocate into the shared_ptr
338         originalFinest_ = ConstructionTraits<MemberType>::construct(args);
339         coarsest_ = std::allocate_shared<Element>(allocator_);
340         coarsest_->element_ = originalFinest_;
341         finest_ = coarsest_;
342       }else{
343         auto old_coarsest = coarsest_;
344         coarsest_ = std::allocate_shared<Element>(allocator_);
345         coarsest_->finer_ = old_coarsest;
346         coarsest_->element_ = ConstructionTraits<MemberType>::construct(args);
347         old_coarsest->coarser_ = coarsest_;
348       }
349       ++levels_;
350     }
351 
352 
353     template<class T, class A>
addFiner(Arguments & args)354     void Hierarchy<T,A>::addFiner(Arguments& args)
355     {
356       //TODO: wouldn't it be better to do this in the constructor?'
357       if(!finest_) {
358         // we have no levels at all...
359         assert(!coarsest_);
360         // allocate into the shared_ptr
361         originalFinest_ = ConstructionTraits<MemberType>::construct(args);
362         finest_ = std::allocate_shared<Element>(allocator_);
363         finest_->element = originalFinest_;
364         coarsest_ = finest_;
365       }else{
366         finest_->finer_ = std::allocate_shared<Element>(allocator_);
367         finest_->finer_->coarser_ = finest_;
368         finest_ = finest_->finer_;
369         finest_->element = ConstructionTraits<T>::construct(args);
370       }
371       ++levels_;
372     }
373 
374     template<class T, class A>
finest()375     typename Hierarchy<T,A>::Iterator Hierarchy<T,A>::finest()
376     {
377       return Iterator(finest_);
378     }
379 
380     template<class T, class A>
coarsest()381     typename Hierarchy<T,A>::Iterator Hierarchy<T,A>::coarsest()
382     {
383       return Iterator(coarsest_);
384     }
385 
386     template<class T, class A>
finest() const387     typename Hierarchy<T,A>::ConstIterator Hierarchy<T,A>::finest() const
388     {
389       return ConstIterator(finest_);
390     }
391 
392     template<class T, class A>
coarsest() const393     typename Hierarchy<T,A>::ConstIterator Hierarchy<T,A>::coarsest() const
394     {
395       return ConstIterator(coarsest_);
396     }
397     /** @} */
398   } // namespace Amg
399 } // namespace Dune
400 
401 #endif
402