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_ISTL_MULTITYPEBLOCKVECTOR_HH
4 #define DUNE_ISTL_MULTITYPEBLOCKVECTOR_HH
5 
6 #include <cmath>
7 #include <iostream>
8 #include <tuple>
9 
10 #include <dune/common/dotproduct.hh>
11 #include <dune/common/ftraits.hh>
12 #include <dune/common/hybridutilities.hh>
13 #include <dune/common/typetraits.hh>
14 
15 #include "istlexception.hh"
16 
17 // forward declaration
18 namespace Dune {
19   template < typename... Args >
20   class MultiTypeBlockVector;
21 }
22 
23 #include "gsetc.hh"
24 
25 namespace Dune {
26 
27   /**
28       @addtogroup ISTL_SPMV
29       @{
30    */
31 
32 
33 
34 
35   /** @addtogroup DenseMatVec
36       @{
37    */
38   template <typename Arg0, typename... Args>
39   struct FieldTraits< MultiTypeBlockVector<Arg0, Args...> >
40   {
41     typedef typename FieldTraits<Arg0>::field_type field_type;
42     typedef typename FieldTraits<Arg0>::real_type real_type;
43   };
44   /**
45       @}
46    */
47 
48   /**
49       @brief A Vector class to support different block types
50 
51       This vector class combines elements of different types known at compile-time.
52    */
53   template < typename... Args >
54   class MultiTypeBlockVector
55   : public std::tuple<Args...>
56   {
57     /** \brief Helper type */
58     typedef std::tuple<Args...> TupleType;
59   public:
60 
61     /** \brief Type used for vector sizes */
62     using size_type = std::size_t;
63 
64     /**
65      * \brief Get the constructors from tuple
66      */
67     using TupleType::TupleType;
68 
69     /**
70      * own class' type
71      */
72     typedef MultiTypeBlockVector<Args...> type;
73 
74     /** \brief The type used for scalars
75      *
76      * The current code hardwires it to 'double', which is far from nice.
77      * On the other hand, it is not clear what the correct type is.  If the MultiTypeBlockVector class
78      * is instantiated with several vectors of different field_types, what should the resulting
79      * field_type be?
80      */
81     typedef double field_type;
82 
83     /** \brief Return the number of non-zero vector entries
84      *
85      * As this is a dense vector data structure, all entries are non-zero,
86      * and hence 'size' returns the same number as 'N'.
87      */
size()88     static constexpr size_type size()
89     {
90       return sizeof...(Args);
91     }
92 
93     /** \brief Number of elements
94      */
N()95     static constexpr size_type N()
96     {
97       return sizeof...(Args);
98     }
99 
100     /**
101      * number of elements
102      *
103      * \deprecated Use method <code>N</code> instead.
104      *             This will be removed after Dune 2.8.
105      */
106     [[deprecated("Use method 'N' instead")]]
count() const107     int count() const
108     {
109       return sizeof...(Args);
110     }
111 
112     /** \brief Number of scalar elements */
dim() const113     size_type dim() const
114     {
115       size_type result = 0;
116       Hybrid::forEach(std::make_index_sequence<N()>{},
117                       [&](auto i){result += std::get<i>(*this).dim();});
118 
119       return result;
120     }
121 
122     /** \brief Random-access operator
123      *
124      * This method mimicks the behavior of normal vector access with square brackets like, e.g., v[5] = 1.
125      * The problem is that the return type is different for each value of the argument in the brackets.
126      * Therefore we implement a trick using std::integral_constant.  To access the first entry of
127      * a MultiTypeBlockVector named v write
128      * \code
129      *  MultiTypeBlockVector<A,B,C,D> v;
130      *  std::integral_constant<std::size_t,0> _0;
131      *  v[_0] = ...
132      * \endcode
133      * The name '_0' used here as a static replacement of the integer number zero is arbitrary.
134      * Any other variable name can be used.  If you don't like the separate variable, you can writee
135      * \code
136      *  MultiTypeBlockVector<A,B,C,D> v;
137      *  v[std::integral_constant<std::size_t,0>()] = ...
138      * \endcode
139      */
140     template< size_type index >
141     typename std::tuple_element<index,TupleType>::type&
operator [](const std::integral_constant<size_type,index> indexVariable)142     operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable)
143     {
144       return std::get<index>(*this);
145     }
146 
147     /** \brief Const random-access operator
148      *
149      * This is the const version of the random-access operator.  See the non-const version for a full
150      * explanation of how to use it.
151      */
152     template< size_type index >
153     const typename std::tuple_element<index,TupleType>::type&
operator [](const std::integral_constant<size_type,index> indexVariable) const154     operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable) const
155     {
156       return std::get<index>(*this);
157     }
158 
159     /** \brief Assignment operator
160      */
161     template<typename T>
operator =(const T & newval)162     void operator= (const T& newval) {
163       Dune::Hybrid::forEach(*this, [&](auto&& entry) {
164         entry = newval;
165       });
166     }
167 
168     /**
169      * operator for MultiTypeBlockVector += MultiTypeBlockVector operations
170      */
operator +=(const type & newv)171     void operator+= (const type& newv) {
172       using namespace Dune::Hybrid;
173       forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
174         (*this)[i] += newv[i];
175       });
176     }
177 
178     /**
179      * operator for MultiTypeBlockVector -= MultiTypeBlockVector operations
180      */
operator -=(const type & newv)181     void operator-= (const type& newv) {
182       using namespace Dune::Hybrid;
183       forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
184         (*this)[i] -= newv[i];
185       });
186     }
187 
188     /** \brief Multiplication with a scalar */
189     template<class T,
190              std::enable_if_t< IsNumber<T>::value, int> = 0>
operator *=(const T & w)191     void operator*= (const T& w) {
192       Hybrid::forEach(*this, [&](auto&& entry) {
193         entry *= w;
194       });
195     }
196 
197     /** \brief Division by a scalar */
198     template<class T,
199              std::enable_if_t< IsNumber<T>::value, int> = 0>
operator /=(const T & w)200     void operator/= (const T& w) {
201       Hybrid::forEach(*this, [&](auto&& entry) {
202         entry /= w;
203       });
204     }
205 
operator *(const type & newv) const206     field_type operator* (const type& newv) const {
207       using namespace Dune::Hybrid;
208       return accumulate(integralRange(Hybrid::size(*this)), field_type(0), [&](auto&& a, auto&& i) {
209         return a + (*this)[i]*newv[i];
210       });
211     }
212 
dot(const type & newv) const213     field_type dot (const type& newv) const {
214       using namespace Dune::Hybrid;
215       return accumulate(integralRange(Hybrid::size(*this)), field_type(0), [&](auto&& a, auto&& i) {
216         return a + (*this)[i].dot(newv[i]);
217       });
218     }
219 
220     /** \brief Compute the 1-norm
221      */
one_norm() const222     auto one_norm() const {
223       using namespace Dune::Hybrid;
224       return accumulate(*this, typename FieldTraits<field_type>::real_type(0), [&](auto&& a, auto&& entry) {
225         return a + entry.one_norm();
226       });
227     }
228 
229     /** \brief Compute the simplified 1-norm (uses 1-norm also for complex values)
230      */
one_norm_real() const231     auto one_norm_real() const {
232       using namespace Dune::Hybrid;
233       return accumulate(*this, typename FieldTraits<field_type>::real_type(0), [&](auto&& a, auto&& entry) {
234         return a + entry.one_norm_real();
235       });
236     }
237 
238     /** \brief Compute the squared Euclidean norm
239      */
two_norm2() const240     typename FieldTraits<field_type>::real_type two_norm2() const {
241       using namespace Dune::Hybrid;
242       return accumulate(*this, typename FieldTraits<field_type>::real_type(0), [&](auto&& a, auto&& entry) {
243         return a + entry.two_norm2();
244       });
245     }
246 
247     /** \brief Compute the Euclidean norm
248      */
two_norm() const249     typename FieldTraits<field_type>::real_type two_norm() const {return sqrt(this->two_norm2());}
250 
251     /** \brief Compute the maximum norm
252      */
infinity_norm() const253     typename FieldTraits<field_type>::real_type infinity_norm() const
254     {
255       using namespace Dune::Hybrid;
256       using std::max;
257       using real_type = typename FieldTraits<field_type>::real_type;
258 
259       real_type result = 0.0;
260       // Compute max norm tracking appearing nan values
261       // if the field type supports nan.
262       if constexpr (HasNaN<field_type>()) {
263         // This variable will preserve any nan value
264         real_type nanTracker = 1.0;
265         using namespace Dune::Hybrid; // needed for icc, see issue #31
266         forEach(*this, [&](auto&& entry) {
267           real_type entryNorm = entry.infinity_norm();
268           result = max(entryNorm, result);
269           nanTracker += entryNorm;
270         });
271         // Incorporate possible nan value into result
272         result *= (nanTracker / nanTracker);
273       } else {
274         using namespace Dune::Hybrid; // needed for icc, see issue #31
275         forEach(*this, [&](auto&& entry) {
276           result = max(entry.infinity_norm(), result);
277         });
278       }
279       return result;
280     }
281 
282     /** \brief Compute the simplified maximum norm (uses 1-norm for complex values)
283      */
infinity_norm_real() const284     typename FieldTraits<field_type>::real_type infinity_norm_real() const
285     {
286       using namespace Dune::Hybrid;
287       using std::max;
288       using real_type = typename FieldTraits<field_type>::real_type;
289 
290       real_type result = 0.0;
291       // Compute max norm tracking appearing nan values
292       // if the field type supports nan.
293       if constexpr (HasNaN<field_type>()) {
294         // This variable will preserve any nan value
295         real_type nanTracker = 1.0;
296         using namespace Dune::Hybrid; // needed for icc, see issue #31
297         forEach(*this, [&](auto&& entry) {
298           real_type entryNorm = entry.infinity_norm_real();
299           result = max(entryNorm, result);
300           nanTracker += entryNorm;
301         });
302         // Incorporate possible nan value into result
303         result *= (nanTracker / nanTracker);
304       } else {
305         using namespace Dune::Hybrid; // needed for icc, see issue #31
306         forEach(*this, [&](auto&& entry) {
307           result = max(entry.infinity_norm_real(), result);
308         });
309       }
310       return result;
311     }
312 
313     /** \brief Axpy operation on this vector (*this += a * y)
314      *
315      * \tparam Ta Type of the scalar 'a'
316      */
317     template<typename Ta>
axpy(const Ta & a,const type & y)318     void axpy (const Ta& a, const type& y) {
319       using namespace Dune::Hybrid;
320       forEach(integralRange(Hybrid::size(*this)), [&](auto&& i) {
321         (*this)[i].axpy(a, y[i]);
322       });
323     }
324 
325   };
326 
327 
328 
329   /** \brief Send MultiTypeBlockVector to an outstream
330    */
331   template <typename... Args>
operator <<(std::ostream & s,const MultiTypeBlockVector<Args...> & v)332   std::ostream& operator<< (std::ostream& s, const MultiTypeBlockVector<Args...>& v) {
333     using namespace Dune::Hybrid;
334     forEach(integralRange(Dune::Hybrid::size(v)), [&](auto&& i) {
335       s << "\t(" << i << "):\n" << v[i] << "\n";
336     });
337     return s;
338   }
339 
340 
341 
342 } // end namespace
343 
344 #endif
345