dune-common  2.9.0
densevector.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 // SPDX-FileCopyrightInfo: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
4 // SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
5 #ifndef DUNE_DENSEVECTOR_HH
6 #define DUNE_DENSEVECTOR_HH
7 
8 #include <algorithm>
9 #include <limits>
10 #include <type_traits>
11 
12 #include "genericiterator.hh"
13 #include "ftraits.hh"
14 #include "matvectraits.hh"
15 #include "promotiontraits.hh"
16 #include "dotproduct.hh"
17 #include "boundschecking.hh"
18 
19 namespace Dune {
20 
21  // forward declaration of template
22  template<typename V> class DenseVector;
23 
24  template<typename V>
25  struct FieldTraits< DenseVector<V> >
26  {
29  };
30 
40  namespace fvmeta
41  {
46  template<class K>
47  inline typename FieldTraits<K>::real_type absreal (const K& k)
48  {
49  using std::abs;
50  return abs(k);
51  }
52 
57  template<class K>
58  inline typename FieldTraits<K>::real_type absreal (const std::complex<K>& c)
59  {
60  using std::abs;
61  return abs(c.real()) + abs(c.imag());
62  }
63 
68  template<class K>
69  inline typename FieldTraits<K>::real_type abs2 (const K& k)
70  {
71  return k*k;
72  }
73 
78  template<class K>
79  inline typename FieldTraits<K>::real_type abs2 (const std::complex<K>& c)
80  {
81  return c.real()*c.real() + c.imag()*c.imag();
82  }
83 
88  template<class K, bool isInteger = std::numeric_limits<K>::is_integer>
89  struct Sqrt
90  {
91  static inline typename FieldTraits<K>::real_type sqrt (const K& k)
92  {
93  using std::sqrt;
94  return sqrt(k);
95  }
96  };
97 
102  template<class K>
103  struct Sqrt<K, true>
104  {
105  static inline typename FieldTraits<K>::real_type sqrt (const K& k)
106  {
107  using std::sqrt;
108  return typename FieldTraits<K>::real_type(sqrt(double(k)));
109  }
110  };
111 
116  template<class K>
117  inline typename FieldTraits<K>::real_type sqrt (const K& k)
118  {
119  return Sqrt<K>::sqrt(k);
120  }
121 
122  }
123 
128  template<class C, class T, class R =T&>
130  public Dune::RandomAccessIteratorFacade<DenseIterator<C,T,R>,T, R, std::ptrdiff_t>
131  {
132  friend class DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type, typename mutable_reference<R>::type >;
133  friend class DenseIterator<const typename std::remove_const<C>::type, const typename std::remove_const<T>::type, typename const_reference<R>::type >;
134 
135  typedef DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type, typename mutable_reference<R>::type > MutableIterator;
136  typedef DenseIterator<const typename std::remove_const<C>::type, const typename std::remove_const<T>::type, typename const_reference<R>::type > ConstIterator;
137  public:
138 
142  typedef std::ptrdiff_t DifferenceType;
143 
147  typedef typename C::size_type SizeType;
148 
149  // Constructors needed by the base iterators.
151  : container_(0), position_()
152  {}
153 
154  DenseIterator(C& cont, SizeType pos)
155  : container_(&cont), position_(pos)
156  {}
157 
159  : container_(other.container_), position_(other.position_)
160  {}
161 
163  : container_(other.container_), position_(other.position_)
164  {}
165 
166  // Methods needed by the forward iterator
167  bool equals(const MutableIterator &other) const
168  {
169  return position_ == other.position_ && container_ == other.container_;
170  }
171 
172 
173  bool equals(const ConstIterator & other) const
174  {
175  return position_ == other.position_ && container_ == other.container_;
176  }
177 
178  R dereference() const {
179  return container_->operator[](position_);
180  }
181 
182  void increment(){
183  ++position_;
184  }
185 
186  // Additional function needed by BidirectionalIterator
187  void decrement(){
188  --position_;
189  }
190 
191  // Additional function needed by RandomAccessIterator
193  return container_->operator[](position_+i);
194  }
195 
197  position_=position_+n;
198  }
199 
200  DifferenceType distanceTo(DenseIterator<const typename std::remove_const<C>::type,const typename std::remove_const<T>::type> other) const
201  {
202  assert(other.container_==container_);
203  return static_cast< DifferenceType >( other.position_ ) - static_cast< DifferenceType >( position_ );
204  }
205 
206  DifferenceType distanceTo(DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type> other) const
207  {
208  assert(other.container_==container_);
209  return static_cast< DifferenceType >( other.position_ ) - static_cast< DifferenceType >( position_ );
210  }
211 
213  SizeType index () const
214  {
215  return this->position_;
216  }
217 
218  private:
219  C *container_;
220  SizeType position_;
221  };
222 
227  template<typename V>
229  {
231  // typedef typename Traits::value_type K;
232 
233  // Curiously recurring template pattern
234  V & asImp() { return static_cast<V&>(*this); }
235  const V & asImp() const { return static_cast<const V&>(*this); }
236 
237  protected:
238  // construction allowed to derived classes only
239  constexpr DenseVector() = default;
240  // copying only allowed by derived classes
241  DenseVector(const DenseVector&) = default;
242 
243  public:
244  //===== type definitions and constants
245 
247  typedef typename Traits::derived_type derived_type;
248 
250  typedef typename Traits::value_type value_type;
251 
254 
256  typedef typename Traits::value_type block_type;
257 
259  typedef typename Traits::size_type size_type;
260 
262  constexpr static int blocklevel = 1;
263 
264  //===== assignment from scalar
267  {
268  for (size_type i=0; i<size(); i++)
269  asImp()[i] = k;
270  return asImp();
271  }
272 
273  //===== assignment from other DenseVectors
274  protected:
276  DenseVector& operator=(const DenseVector&) = default;
277 
278  public:
279 
281  template <typename W,
282  std::enable_if_t<
283  std::is_assignable<value_type&, typename DenseVector<W>::value_type>::value, int> = 0>
285  {
286  assert(other.size() == size());
287  for (size_type i=0; i<size(); i++)
288  asImp()[i] = other[i];
289  return asImp();
290  }
291 
292  //===== access to components
293 
296  {
297  return asImp()[i];
298  }
299 
301  {
302  return asImp()[i];
303  }
304 
307  {
308  return asImp()[0];
309  }
310 
312  const value_type& front() const
313  {
314  return asImp()[0];
315  }
316 
319  {
320  return asImp()[size()-1];
321  }
322 
324  const value_type& back() const
325  {
326  return asImp()[size()-1];
327  }
328 
330  bool empty() const
331  {
332  return size() == 0;
333  }
334 
336  size_type size() const
337  {
338  return asImp().size();
339  }
340 
345 
348  {
349  return Iterator(*this,0);
350  }
351 
354  {
355  return Iterator(*this,size());
356  }
357 
361  {
362  return Iterator(*this,size()-1);
363  }
364 
368  {
369  return Iterator(*this,-1);
370  }
371 
374  {
375  return Iterator(*this,std::min(i,size()));
376  }
377 
382 
385  {
386  return ConstIterator(*this,0);
387  }
388 
391  {
392  return ConstIterator(*this,size());
393  }
394 
398  {
399  return ConstIterator(*this,size()-1);
400  }
401 
405  {
406  return ConstIterator(*this,-1);
407  }
408 
411  {
412  return ConstIterator(*this,std::min(i,size()));
413  }
414 
415  //===== vector space arithmetic
416 
418  template <class Other>
420  {
421  DUNE_ASSERT_BOUNDS(x.size() == size());
422  for (size_type i=0; i<size(); i++)
423  (*this)[i] += x[i];
424  return asImp();
425  }
426 
428  template <class Other>
430  {
431  DUNE_ASSERT_BOUNDS(x.size() == size());
432  for (size_type i=0; i<size(); i++)
433  (*this)[i] -= x[i];
434  return asImp();
435  }
436 
438  template <class Other>
440  {
441  derived_type z = asImp();
442  return (z+=b);
443  }
444 
446  template <class Other>
448  {
449  derived_type z = asImp();
450  return (z-=b);
451  }
452 
455  {
456  V result;
457  using idx_type = typename decltype(result)::size_type;
458 
459  for (idx_type i = 0; i < size(); ++i)
460  result[i] = -asImp()[i];
461 
462  return result;
463  }
464 
466 
474  template <typename ValueType>
475  typename std::enable_if<
476  std::is_convertible<ValueType, value_type>::value,
478  >::type&
479  operator+= (const ValueType& kk)
480  {
481  const value_type& k = kk;
482  for (size_type i=0; i<size(); i++)
483  (*this)[i] += k;
484  return asImp();
485  }
486 
488 
496  template <typename ValueType>
497  typename std::enable_if<
498  std::is_convertible<ValueType, value_type>::value,
500  >::type&
501  operator-= (const ValueType& kk)
502  {
503  const value_type& k = kk;
504  for (size_type i=0; i<size(); i++)
505  (*this)[i] -= k;
506  return asImp();
507  }
508 
510 
518  template <typename FieldType>
519  typename std::enable_if<
520  std::is_convertible<FieldType, field_type>::value,
522  >::type&
523  operator*= (const FieldType& kk)
524  {
525  const field_type& k = kk;
526  for (size_type i=0; i<size(); i++)
527  (*this)[i] *= k;
528  return asImp();
529  }
530 
532 
540  template <typename FieldType>
541  typename std::enable_if<
542  std::is_convertible<FieldType, field_type>::value,
544  >::type&
545  operator/= (const FieldType& kk)
546  {
547  const field_type& k = kk;
548  for (size_type i=0; i<size(); i++)
549  (*this)[i] /= k;
550  return asImp();
551  }
552 
554  template <class Other>
555  bool operator== (const DenseVector<Other>& x) const
556  {
557  DUNE_ASSERT_BOUNDS(x.size() == size());
558  for (size_type i=0; i<size(); i++)
559  if ((*this)[i]!=x[i])
560  return false;
561 
562  return true;
563  }
564 
566  template <class Other>
567  bool operator!= (const DenseVector<Other>& x) const
568  {
569  return !operator==(x);
570  }
571 
572 
574  template <class Other>
576  {
577  DUNE_ASSERT_BOUNDS(x.size() == size());
578  for (size_type i=0; i<size(); i++)
579  (*this)[i] += a*x[i];
580  return asImp();
581  }
582 
590  template<class Other>
592  typedef typename PromotionTraits<field_type, typename DenseVector<Other>::field_type>::PromotedType PromotedType;
593  PromotedType result(0);
594  assert(x.size() == size());
595  for (size_type i=0; i<size(); i++) {
596  result += PromotedType((*this)[i]*x[i]);
597  }
598  return result;
599  }
600 
608  template<class Other>
610  typedef typename PromotionTraits<field_type, typename DenseVector<Other>::field_type>::PromotedType PromotedType;
611  PromotedType result(0);
612  assert(x.size() == size());
613  for (size_type i=0; i<size(); i++) {
614  result += Dune::dot((*this)[i],x[i]);
615  }
616  return result;
617  }
618 
619  //===== norms
620 
623  using std::abs;
624  typename FieldTraits<value_type>::real_type result( 0 );
625  for (size_type i=0; i<size(); i++)
626  result += abs((*this)[i]);
627  return result;
628  }
629 
630 
633  {
634  typename FieldTraits<value_type>::real_type result( 0 );
635  for (size_type i=0; i<size(); i++)
636  result += fvmeta::absreal((*this)[i]);
637  return result;
638  }
639 
642  {
643  typename FieldTraits<value_type>::real_type result( 0 );
644  for (size_type i=0; i<size(); i++)
645  result += fvmeta::abs2((*this)[i]);
646  return fvmeta::sqrt(result);
647  }
648 
651  {
652  typename FieldTraits<value_type>::real_type result( 0 );
653  for (size_type i=0; i<size(); i++)
654  result += fvmeta::abs2((*this)[i]);
655  return result;
656  }
657 
659  template <typename vt = value_type,
660  typename std::enable_if<!HasNaN<vt>::value, int>::type = 0>
662  using real_type = typename FieldTraits<vt>::real_type;
663  using std::abs;
664  using std::max;
665 
666  real_type norm = 0;
667  for (auto const &x : *this) {
668  real_type const a = abs(x);
669  norm = max(a, norm);
670  }
671  return norm;
672  }
673 
675  template <typename vt = value_type,
676  typename std::enable_if<!HasNaN<vt>::value, int>::type = 0>
678  using real_type = typename FieldTraits<vt>::real_type;
679  using std::max;
680 
681  real_type norm = 0;
682  for (auto const &x : *this) {
683  real_type const a = fvmeta::absreal(x);
684  norm = max(a, norm);
685  }
686  return norm;
687  }
688 
690  template <typename vt = value_type,
691  typename std::enable_if<HasNaN<vt>::value, int>::type = 0>
693  using real_type = typename FieldTraits<vt>::real_type;
694  using std::abs;
695  using std::max;
696 
697  real_type norm = 0;
698  real_type isNaN = 1;
699  for (auto const &x : *this) {
700  real_type const a = abs(x);
701  norm = max(a, norm);
702  isNaN += a;
703  }
704  return norm * (isNaN / isNaN);
705  }
706 
708  template <typename vt = value_type,
709  typename std::enable_if<HasNaN<vt>::value, int>::type = 0>
711  using real_type = typename FieldTraits<vt>::real_type;
712  using std::max;
713 
714  real_type norm = 0;
715  real_type isNaN = 1;
716  for (auto const &x : *this) {
717  real_type const a = fvmeta::absreal(x);
718  norm = max(a, norm);
719  isNaN += a;
720  }
721  return norm * (isNaN / isNaN);
722  }
723 
724  //===== sizes
725 
727  size_type N () const
728  {
729  return size();
730  }
731 
733  size_type dim () const
734  {
735  return size();
736  }
737 
738  };
739 
748  template<typename V>
749  std::ostream& operator<< (std::ostream& s, const DenseVector<V>& v)
750  {
751  for (typename DenseVector<V>::size_type i=0; i<v.size(); i++)
752  s << ((i>0) ? " " : "") << v[i];
753  return s;
754  }
755 
758 } // end namespace
759 
760 #endif // DUNE_DENSEVECTOR_HH
Macro for wrapping boundary checks.
Provides the functions dot(a,b) := and dotT(a,b) := .
Type traits to determine the type of reals (when working with complex numbers)
Implements a generic iterator class for writing stl conformant iterators.
Documentation of the traits classes you need to write for each implementation of DenseVector or Dense...
Compute type of the result of an arithmetic operation involving two different number types.
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition: boundschecking.hh:30
auto dot(const A &a, const B &b) -> typename std::enable_if< IsNumber< A >::value &&!IsVector< A >::value &&!std::is_same< typename FieldTraits< A >::field_type, typename FieldTraits< A >::real_type > ::value, decltype(conj(a) *b)>::type
computes the dot product for fundamental data types according to Petsc's VectDot function: dot(a,...
Definition: dotproduct.hh:42
std::ostream & operator<<(std::ostream &s, const bigunsignedint< k > &x)
Definition: bigunsignedint.hh:278
Dune namespace.
Definition: alignedallocator.hh:13
auto min(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:447
auto max(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:425
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Definition: fvector.hh:604
Interface for a class of dense vectors over a given field.
Definition: densevector.hh:229
Traits::value_type value_type
export the type representing the field
Definition: densevector.hh:250
value_type & back()
return reference to last element
Definition: densevector.hh:318
ConstIterator const_iterator
typedef for stl compliant access
Definition: densevector.hh:381
Iterator iterator
typedef for stl compliant access
Definition: densevector.hh:344
std::enable_if< std::is_convertible< FieldType, field_type >::value, derived_type >::type & operator*=(const FieldType &kk)
vector space multiplication with scalar
Definition: densevector.hh:523
ConstIterator find(size_type i) const
return iterator to given element or end()
Definition: densevector.hh:410
ConstIterator end() const
end ConstIterator
Definition: densevector.hh:390
ConstIterator beforeBegin() const
Definition: densevector.hh:404
bool operator==(const DenseVector< Other > &x) const
Binary vector comparison.
Definition: densevector.hh:555
FieldTraits< value_type >::real_type two_norm2() const
square of two norm (sum over squared values of entries), need for block recursion
Definition: densevector.hh:650
Iterator begin()
begin iterator
Definition: densevector.hh:347
Iterator beforeBegin()
Definition: densevector.hh:367
DenseIterator< const DenseVector, const value_type > ConstIterator
ConstIterator class for sequential access.
Definition: densevector.hh:379
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:575
Traits::derived_type derived_type
type of derived vector class
Definition: densevector.hh:247
derived_type operator+(const DenseVector< Other > &b) const
Binary vector addition.
Definition: densevector.hh:439
size_type size() const
size method
Definition: densevector.hh:336
const value_type & front() const
return reference to first element
Definition: densevector.hh:312
PromotionTraits< field_type, typename DenseVector< Other >::field_type >::PromotedType operator*(const DenseVector< Other > &x) const
indefinite vector dot product which corresponds to Petsc's VecTDot
Definition: densevector.hh:591
size_type dim() const
dimension of the vector space
Definition: densevector.hh:733
constexpr static int blocklevel
The number of block levels we contain. This is the leaf, that is, 1.
Definition: densevector.hh:262
ConstIterator beforeEnd() const
Definition: densevector.hh:397
FieldTraits< value_type >::real_type one_norm() const
one norm (sum over absolute values of entries)
Definition: densevector.hh:622
Iterator end()
end iterator
Definition: densevector.hh:353
std::enable_if< std::is_convertible< FieldType, field_type >::value, derived_type >::type & operator/=(const FieldType &kk)
vector space division by scalar
Definition: densevector.hh:545
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densevector.hh:259
DenseIterator< DenseVector, value_type > Iterator
Iterator class for sequential access.
Definition: densevector.hh:342
const value_type & back() const
return reference to last element
Definition: densevector.hh:324
DenseVector(const DenseVector &)=default
Iterator beforeEnd()
Definition: densevector.hh:360
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:641
derived_type & operator+=(const DenseVector< Other > &x)
vector space addition
Definition: densevector.hh:419
bool operator!=(const DenseVector< Other > &x) const
Binary vector incomparison.
Definition: densevector.hh:567
value_type & operator[](size_type i)
random access
Definition: densevector.hh:295
ConstIterator begin() const
begin ConstIterator
Definition: densevector.hh:384
derived_type & operator=(const value_type &k)
Assignment operator for scalar.
Definition: densevector.hh:266
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densevector.hh:677
constexpr DenseVector()=default
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (maximum of absolute values of entries)
Definition: densevector.hh:661
FieldTraits< value_type >::field_type field_type
export the type representing the field
Definition: densevector.hh:253
Traits::value_type block_type
export the type representing the components
Definition: densevector.hh:256
PromotionTraits< field_type, typename DenseVector< Other >::field_type >::PromotedType dot(const DenseVector< Other > &x) const
vector dot product which corresponds to Petsc's VecDot
Definition: densevector.hh:609
DenseVector & operator=(const DenseVector &)=default
Assignment operator for other DenseVector of same type.
derived_type operator-() const
Vector negation.
Definition: densevector.hh:454
derived_type & operator-=(const DenseVector< Other > &x)
vector space subtraction
Definition: densevector.hh:429
Iterator find(size_type i)
return iterator to given element or end()
Definition: densevector.hh:373
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition: densevector.hh:727
bool empty() const
checks whether the container is empty
Definition: densevector.hh:330
FieldTraits< value_type >::real_type one_norm_real() const
simplified one norm (uses Manhattan norm for complex values)
Definition: densevector.hh:632
value_type & front()
return reference to first element
Definition: densevector.hh:306
FieldTraits< typename DenseMatVecTraits< V >::value_type >::field_type field_type
Definition: densevector.hh:27
FieldTraits< typename DenseMatVecTraits< V >::value_type >::real_type real_type
Definition: densevector.hh:28
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:131
void increment()
Definition: densevector.hh:182
SizeType index() const
return index
Definition: densevector.hh:213
bool equals(const MutableIterator &other) const
Definition: densevector.hh:167
DenseIterator(const MutableIterator &other)
Definition: densevector.hh:158
bool equals(const ConstIterator &other) const
Definition: densevector.hh:173
R elementAt(DifferenceType i) const
Definition: densevector.hh:192
DifferenceType distanceTo(DenseIterator< const typename std::remove_const< C >::type, const typename std::remove_const< T >::type > other) const
Definition: densevector.hh:200
void decrement()
Definition: densevector.hh:187
DenseIterator(const ConstIterator &other)
Definition: densevector.hh:162
DifferenceType distanceTo(DenseIterator< typename std::remove_const< C >::type, typename std::remove_const< T >::type > other) const
Definition: densevector.hh:206
DenseIterator(C &cont, SizeType pos)
Definition: densevector.hh:154
std::ptrdiff_t DifferenceType
The type of the difference between two positions.
Definition: densevector.hh:142
R dereference() const
Definition: densevector.hh:178
void advance(DifferenceType n)
Definition: densevector.hh:196
C::size_type SizeType
The type to index the underlying container.
Definition: densevector.hh:147
Definition: ftraits.hh:26
T field_type
export the type representing the field
Definition: ftraits.hh:28
T real_type
export the type representing the real type of the field
Definition: ftraits.hh:30
get the 'mutable' version of a reference to a const object
Definition: genericiterator.hh:116
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:434
Definition: matvectraits.hh:31
Compute type of the result of an arithmetic operation involving two different number types.
Definition: promotiontraits.hh:27