dune-common  2.8.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 #ifndef DUNE_DENSEVECTOR_HH
4 #define DUNE_DENSEVECTOR_HH
5 
6 #include <algorithm>
7 #include <limits>
8 #include <type_traits>
9 
10 #include "genericiterator.hh"
11 #include "ftraits.hh"
12 #include "matvectraits.hh"
13 #include "promotiontraits.hh"
14 #include "dotproduct.hh"
15 #include "boundschecking.hh"
16 
17 namespace Dune {
18 
19  // forward declaration of template
20  template<typename V> class DenseVector;
21 
22  template<typename V>
23  struct FieldTraits< DenseVector<V> >
24  {
27  };
28 
38  namespace fvmeta
39  {
44  template<class K>
45  inline typename FieldTraits<K>::real_type absreal (const K& k)
46  {
47  using std::abs;
48  return abs(k);
49  }
50 
55  template<class K>
56  inline typename FieldTraits<K>::real_type absreal (const std::complex<K>& c)
57  {
58  using std::abs;
59  return abs(c.real()) + abs(c.imag());
60  }
61 
66  template<class K>
67  inline typename FieldTraits<K>::real_type abs2 (const K& k)
68  {
69  return k*k;
70  }
71 
76  template<class K>
77  inline typename FieldTraits<K>::real_type abs2 (const std::complex<K>& c)
78  {
79  return c.real()*c.real() + c.imag()*c.imag();
80  }
81 
86  template<class K, bool isInteger = std::numeric_limits<K>::is_integer>
87  struct Sqrt
88  {
89  static inline typename FieldTraits<K>::real_type sqrt (const K& k)
90  {
91  using std::sqrt;
92  return sqrt(k);
93  }
94  };
95 
100  template<class K>
101  struct Sqrt<K, true>
102  {
103  static inline typename FieldTraits<K>::real_type sqrt (const K& k)
104  {
105  using std::sqrt;
106  return typename FieldTraits<K>::real_type(sqrt(double(k)));
107  }
108  };
109 
114  template<class K>
115  inline typename FieldTraits<K>::real_type sqrt (const K& k)
116  {
117  return Sqrt<K>::sqrt(k);
118  }
119 
120  }
121 
126  template<class C, class T, class R =T&>
128  public Dune::RandomAccessIteratorFacade<DenseIterator<C,T,R>,T, R, std::ptrdiff_t>
129  {
130  friend class DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type, typename mutable_reference<R>::type >;
131  friend class DenseIterator<const typename std::remove_const<C>::type, const typename std::remove_const<T>::type, typename const_reference<R>::type >;
132 
133  typedef DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type, typename mutable_reference<R>::type > MutableIterator;
134  typedef DenseIterator<const typename std::remove_const<C>::type, const typename std::remove_const<T>::type, typename const_reference<R>::type > ConstIterator;
135  public:
136 
140  typedef std::ptrdiff_t DifferenceType;
141 
145  typedef typename C::size_type SizeType;
146 
147  // Constructors needed by the base iterators.
149  : container_(0), position_()
150  {}
151 
152  DenseIterator(C& cont, SizeType pos)
153  : container_(&cont), position_(pos)
154  {}
155 
157  : container_(other.container_), position_(other.position_)
158  {}
159 
161  : container_(other.container_), position_(other.position_)
162  {}
163 
164  // Methods needed by the forward iterator
165  bool equals(const MutableIterator &other) const
166  {
167  return position_ == other.position_ && container_ == other.container_;
168  }
169 
170 
171  bool equals(const ConstIterator & other) const
172  {
173  return position_ == other.position_ && container_ == other.container_;
174  }
175 
176  R dereference() const {
177  return container_->operator[](position_);
178  }
179 
180  void increment(){
181  ++position_;
182  }
183 
184  // Additional function needed by BidirectionalIterator
185  void decrement(){
186  --position_;
187  }
188 
189  // Additional function needed by RandomAccessIterator
191  return container_->operator[](position_+i);
192  }
193 
195  position_=position_+n;
196  }
197 
198  DifferenceType distanceTo(DenseIterator<const typename std::remove_const<C>::type,const typename std::remove_const<T>::type> other) const
199  {
200  assert(other.container_==container_);
201  return static_cast< DifferenceType >( other.position_ ) - static_cast< DifferenceType >( position_ );
202  }
203 
204  DifferenceType distanceTo(DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type> other) const
205  {
206  assert(other.container_==container_);
207  return static_cast< DifferenceType >( other.position_ ) - static_cast< DifferenceType >( position_ );
208  }
209 
211  SizeType index () const
212  {
213  return this->position_;
214  }
215 
216  private:
217  C *container_;
218  SizeType position_;
219  };
220 
225  template<typename V>
227  {
229  // typedef typename Traits::value_type K;
230 
231  // Curiously recurring template pattern
232  V & asImp() { return static_cast<V&>(*this); }
233  const V & asImp() const { return static_cast<const V&>(*this); }
234 
235  protected:
236  // construction allowed to derived classes only
237  constexpr DenseVector() = default;
238  // copying only allowed by derived classes
239  DenseVector(const DenseVector&) = default;
240 
241  public:
242  //===== type definitions and constants
243 
245  typedef typename Traits::derived_type derived_type;
246 
248  typedef typename Traits::value_type value_type;
249 
252 
254  typedef typename Traits::value_type block_type;
255 
257  typedef typename Traits::size_type size_type;
258 
260  enum {
262  blocklevel = 1
263  };
264 
265  //===== assignment from scalar
268  {
269  for (size_type i=0; i<size(); i++)
270  asImp()[i] = k;
271  return asImp();
272  }
273 
274  //===== assignment from other DenseVectors
275  protected:
277  DenseVector& operator=(const DenseVector&) = default;
278 
279  public:
280 
282  template <typename W,
283  std::enable_if_t<
284  std::is_assignable<value_type&, typename DenseVector<W>::value_type>::value, int> = 0>
286  {
287  assert(other.size() == size());
288  for (size_type i=0; i<size(); i++)
289  asImp()[i] = other[i];
290  return asImp();
291  }
292 
293  //===== access to components
294 
297  {
298  return asImp()[i];
299  }
300 
302  {
303  return asImp()[i];
304  }
305 
308  {
309  return asImp()[0];
310  }
311 
313  const value_type& front() const
314  {
315  return asImp()[0];
316  }
317 
320  {
321  return asImp()[size()-1];
322  }
323 
325  const value_type& back() const
326  {
327  return asImp()[size()-1];
328  }
329 
331  bool empty() const
332  {
333  return size() == 0;
334  }
335 
337  size_type size() const
338  {
339  return asImp().size();
340  }
341 
346 
349  {
350  return Iterator(*this,0);
351  }
352 
355  {
356  return Iterator(*this,size());
357  }
358 
362  {
363  return Iterator(*this,size()-1);
364  }
365 
369  {
370  return Iterator(*this,-1);
371  }
372 
375  {
376  return Iterator(*this,std::min(i,size()));
377  }
378 
383 
386  {
387  return ConstIterator(*this,0);
388  }
389 
392  {
393  return ConstIterator(*this,size());
394  }
395 
399  {
400  return ConstIterator(*this,size()-1);
401  }
402 
406  {
407  return ConstIterator(*this,-1);
408  }
409 
412  {
413  return ConstIterator(*this,std::min(i,size()));
414  }
415 
416  //===== vector space arithmetic
417 
419  template <class Other>
421  {
422  DUNE_ASSERT_BOUNDS(x.size() == size());
423  for (size_type i=0; i<size(); i++)
424  (*this)[i] += x[i];
425  return asImp();
426  }
427 
429  template <class Other>
431  {
432  DUNE_ASSERT_BOUNDS(x.size() == size());
433  for (size_type i=0; i<size(); i++)
434  (*this)[i] -= x[i];
435  return asImp();
436  }
437 
439  template <class Other>
441  {
442  derived_type z = asImp();
443  return (z+=b);
444  }
445 
447  template <class Other>
449  {
450  derived_type z = asImp();
451  return (z-=b);
452  }
453 
456  {
457  V result;
458  using idx_type = typename decltype(result)::size_type;
459 
460  for (idx_type i = 0; i < size(); ++i)
461  result[i] = -asImp()[i];
462 
463  return result;
464  }
465 
467 
475  template <typename ValueType>
476  typename std::enable_if<
477  std::is_convertible<ValueType, value_type>::value,
479  >::type&
480  operator+= (const ValueType& kk)
481  {
482  const value_type& k = kk;
483  for (size_type i=0; i<size(); i++)
484  (*this)[i] += k;
485  return asImp();
486  }
487 
489 
497  template <typename ValueType>
498  typename std::enable_if<
499  std::is_convertible<ValueType, value_type>::value,
501  >::type&
502  operator-= (const ValueType& kk)
503  {
504  const value_type& k = kk;
505  for (size_type i=0; i<size(); i++)
506  (*this)[i] -= k;
507  return asImp();
508  }
509 
511 
519  template <typename FieldType>
520  typename std::enable_if<
521  std::is_convertible<FieldType, field_type>::value,
523  >::type&
524  operator*= (const FieldType& kk)
525  {
526  const field_type& k = kk;
527  for (size_type i=0; i<size(); i++)
528  (*this)[i] *= k;
529  return asImp();
530  }
531 
533 
541  template <typename FieldType>
542  typename std::enable_if<
543  std::is_convertible<FieldType, field_type>::value,
545  >::type&
546  operator/= (const FieldType& kk)
547  {
548  const field_type& k = kk;
549  for (size_type i=0; i<size(); i++)
550  (*this)[i] /= k;
551  return asImp();
552  }
553 
555  template <class Other>
556  bool operator== (const DenseVector<Other>& x) const
557  {
558  DUNE_ASSERT_BOUNDS(x.size() == size());
559  for (size_type i=0; i<size(); i++)
560  if ((*this)[i]!=x[i])
561  return false;
562 
563  return true;
564  }
565 
567  template <class Other>
568  bool operator!= (const DenseVector<Other>& x) const
569  {
570  return !operator==(x);
571  }
572 
573 
575  template <class Other>
577  {
578  DUNE_ASSERT_BOUNDS(x.size() == size());
579  for (size_type i=0; i<size(); i++)
580  (*this)[i] += a*x[i];
581  return asImp();
582  }
583 
591  template<class Other>
593  typedef typename PromotionTraits<field_type, typename DenseVector<Other>::field_type>::PromotedType PromotedType;
594  PromotedType result(0);
595  assert(x.size() == size());
596  for (size_type i=0; i<size(); i++) {
597  result += PromotedType((*this)[i]*x[i]);
598  }
599  return result;
600  }
601 
609  template<class Other>
611  typedef typename PromotionTraits<field_type, typename DenseVector<Other>::field_type>::PromotedType PromotedType;
612  PromotedType result(0);
613  assert(x.size() == size());
614  for (size_type i=0; i<size(); i++) {
615  result += Dune::dot((*this)[i],x[i]);
616  }
617  return result;
618  }
619 
620  //===== norms
621 
624  using std::abs;
625  typename FieldTraits<value_type>::real_type result( 0 );
626  for (size_type i=0; i<size(); i++)
627  result += abs((*this)[i]);
628  return result;
629  }
630 
631 
634  {
635  typename FieldTraits<value_type>::real_type result( 0 );
636  for (size_type i=0; i<size(); i++)
637  result += fvmeta::absreal((*this)[i]);
638  return result;
639  }
640 
643  {
644  typename FieldTraits<value_type>::real_type result( 0 );
645  for (size_type i=0; i<size(); i++)
646  result += fvmeta::abs2((*this)[i]);
647  return fvmeta::sqrt(result);
648  }
649 
652  {
653  typename FieldTraits<value_type>::real_type result( 0 );
654  for (size_type i=0; i<size(); i++)
655  result += fvmeta::abs2((*this)[i]);
656  return result;
657  }
658 
660  template <typename vt = value_type,
661  typename std::enable_if<!HasNaN<vt>::value, int>::type = 0>
663  using real_type = typename FieldTraits<vt>::real_type;
664  using std::abs;
665  using std::max;
666 
667  real_type norm = 0;
668  for (auto const &x : *this) {
669  real_type const a = abs(x);
670  norm = max(a, norm);
671  }
672  return norm;
673  }
674 
676  template <typename vt = value_type,
677  typename std::enable_if<!HasNaN<vt>::value, int>::type = 0>
679  using real_type = typename FieldTraits<vt>::real_type;
680  using std::max;
681 
682  real_type norm = 0;
683  for (auto const &x : *this) {
684  real_type const a = fvmeta::absreal(x);
685  norm = max(a, norm);
686  }
687  return norm;
688  }
689 
691  template <typename vt = value_type,
692  typename std::enable_if<HasNaN<vt>::value, int>::type = 0>
694  using real_type = typename FieldTraits<vt>::real_type;
695  using std::abs;
696  using std::max;
697 
698  real_type norm = 0;
699  real_type isNaN = 1;
700  for (auto const &x : *this) {
701  real_type const a = abs(x);
702  norm = max(a, norm);
703  isNaN += a;
704  }
705  return norm * (isNaN / isNaN);
706  }
707 
709  template <typename vt = value_type,
710  typename std::enable_if<HasNaN<vt>::value, int>::type = 0>
712  using real_type = typename FieldTraits<vt>::real_type;
713  using std::max;
714 
715  real_type norm = 0;
716  real_type isNaN = 1;
717  for (auto const &x : *this) {
718  real_type const a = fvmeta::absreal(x);
719  norm = max(a, norm);
720  isNaN += a;
721  }
722  return norm * (isNaN / isNaN);
723  }
724 
725  //===== sizes
726 
728  size_type N () const
729  {
730  return size();
731  }
732 
734  size_type dim () const
735  {
736  return size();
737  }
738 
739  };
740 
749  template<typename V>
750  std::ostream& operator<< (std::ostream& s, const DenseVector<V>& v)
751  {
752  for (typename DenseVector<V>::size_type i=0; i<v.size(); i++)
753  s << ((i>0) ? " " : "") << v[i];
754  return s;
755  }
756 
759 } // end namespace
760 
761 #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.
auto dot(const A &a, const B &b) -> typename std::enable_if<!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:40
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition: boundschecking.hh:28
std::ostream & operator<<(std::ostream &s, const bigunsignedint< k > &x)
Definition: bigunsignedint.hh:273
Dune namespace.
Definition: alignedallocator.hh:11
auto min(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:434
auto max(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:412
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Definition: fvector.hh:610
Interface for a class of dense vectors over a given field.
Definition: densevector.hh:227
Traits::value_type value_type
export the type representing the field
Definition: densevector.hh:248
value_type & back()
return reference to last element
Definition: densevector.hh:319
ConstIterator const_iterator
typedef for stl compliant access
Definition: densevector.hh:382
Iterator iterator
typedef for stl compliant access
Definition: densevector.hh:345
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:524
ConstIterator find(size_type i) const
return iterator to given element or end()
Definition: densevector.hh:411
ConstIterator end() const
end ConstIterator
Definition: densevector.hh:391
ConstIterator beforeBegin() const
Definition: densevector.hh:405
bool operator==(const DenseVector< Other > &x) const
Binary vector comparison.
Definition: densevector.hh:556
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:651
Iterator begin()
begin iterator
Definition: densevector.hh:348
Iterator beforeBegin()
Definition: densevector.hh:368
DenseIterator< const DenseVector, const value_type > ConstIterator
ConstIterator class for sequential access.
Definition: densevector.hh:380
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:576
Traits::derived_type derived_type
type of derived vector class
Definition: densevector.hh:245
derived_type operator+(const DenseVector< Other > &b) const
Binary vector addition.
Definition: densevector.hh:440
size_type size() const
size method
Definition: densevector.hh:337
const value_type & front() const
return reference to first element
Definition: densevector.hh:313
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:592
size_type dim() const
dimension of the vector space
Definition: densevector.hh:734
ConstIterator beforeEnd() const
Definition: densevector.hh:398
FieldTraits< value_type >::real_type one_norm() const
one norm (sum over absolute values of entries)
Definition: densevector.hh:623
Iterator end()
end iterator
Definition: densevector.hh:354
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:546
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densevector.hh:257
DenseIterator< DenseVector, value_type > Iterator
Iterator class for sequential access.
Definition: densevector.hh:343
const value_type & back() const
return reference to last element
Definition: densevector.hh:325
DenseVector(const DenseVector &)=default
Iterator beforeEnd()
Definition: densevector.hh:361
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:642
derived_type & operator+=(const DenseVector< Other > &x)
vector space addition
Definition: densevector.hh:420
bool operator!=(const DenseVector< Other > &x) const
Binary vector incomparison.
Definition: densevector.hh:568
value_type & operator[](size_type i)
random access
Definition: densevector.hh:296
ConstIterator begin() const
begin ConstIterator
Definition: densevector.hh:385
derived_type & operator=(const value_type &k)
Assignment operator for scalar.
Definition: densevector.hh:267
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densevector.hh:678
constexpr DenseVector()=default
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (maximum of absolute values of entries)
Definition: densevector.hh:662
FieldTraits< value_type >::field_type field_type
export the type representing the field
Definition: densevector.hh:251
Traits::value_type block_type
export the type representing the components
Definition: densevector.hh:254
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:610
@ blocklevel
The number of block levels we contain.
Definition: densevector.hh:262
DenseVector & operator=(const DenseVector &)=default
Assignment operator for other DenseVector of same type.
derived_type operator-() const
Vector negation.
Definition: densevector.hh:455
derived_type & operator-=(const DenseVector< Other > &x)
vector space subtraction
Definition: densevector.hh:430
Iterator find(size_type i)
return iterator to given element or end()
Definition: densevector.hh:374
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition: densevector.hh:728
bool empty() const
checks whether the container is empty
Definition: densevector.hh:331
FieldTraits< value_type >::real_type one_norm_real() const
simplified one norm (uses Manhattan norm for complex values)
Definition: densevector.hh:633
value_type & front()
return reference to first element
Definition: densevector.hh:307
FieldTraits< typename DenseMatVecTraits< V >::value_type >::field_type field_type
Definition: densevector.hh:25
FieldTraits< typename DenseMatVecTraits< V >::value_type >::real_type real_type
Definition: densevector.hh:26
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:129
void increment()
Definition: densevector.hh:180
SizeType index() const
return index
Definition: densevector.hh:211
bool equals(const MutableIterator &other) const
Definition: densevector.hh:165
DenseIterator(const MutableIterator &other)
Definition: densevector.hh:156
bool equals(const ConstIterator &other) const
Definition: densevector.hh:171
R elementAt(DifferenceType i) const
Definition: densevector.hh:190
DifferenceType distanceTo(DenseIterator< const typename std::remove_const< C >::type, const typename std::remove_const< T >::type > other) const
Definition: densevector.hh:198
void decrement()
Definition: densevector.hh:185
DenseIterator(const ConstIterator &other)
Definition: densevector.hh:160
DifferenceType distanceTo(DenseIterator< typename std::remove_const< C >::type, typename std::remove_const< T >::type > other) const
Definition: densevector.hh:204
DenseIterator(C &cont, SizeType pos)
Definition: densevector.hh:152
std::ptrdiff_t DifferenceType
The type of the difference between two positions.
Definition: densevector.hh:140
R dereference() const
Definition: densevector.hh:176
void advance(DifferenceType n)
Definition: densevector.hh:194
C::size_type SizeType
The type to index the underlying container.
Definition: densevector.hh:145
Definition: ftraits.hh:24
T field_type
export the type representing the field
Definition: ftraits.hh:26
T real_type
export the type representing the real type of the field
Definition: ftraits.hh:28
get the 'mutable' version of a reference to a const object
Definition: genericiterator.hh:114
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:432
Definition: matvectraits.hh:29
Compute type of the result of an arithmetic operation involving two different number types.
Definition: promotiontraits.hh:25