Dune Core Modules (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
17namespace Dune {
18
19 // forward declaration of template
20 template<typename V> class DenseVector;
21
22 template<typename V>
23 struct FieldTraits< DenseVector<V> >
24 {
25 typedef typename FieldTraits< typename DenseMatVecTraits<V>::value_type >::field_type field_type;
26 typedef typename FieldTraits< typename DenseMatVecTraits<V>::value_type >::real_type real_type;
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
156 DenseIterator(const MutableIterator & other)
157 : container_(other.container_), position_(other.position_)
158 {}
159
160 DenseIterator(const ConstIterator & other)
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
190 R elementAt(DifferenceType i) const {
191 return container_->operator[](position_+i);
192 }
193
194 void advance(DifferenceType n){
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
212 {
213 return this->position_;
214 }
215
216 private:
217 C *container_;
218 SizeType position_;
219 };
220
225 template<typename V>
227 {
228 typedef DenseMatVecTraits<V> Traits;
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
251 typedef typename FieldTraits< value_type >::field_type field_type;
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:
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
301 const value_type & operator[] (size_type i) const
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
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
623 typename FieldTraits<value_type>::real_type one_norm() const {
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
633 typename FieldTraits<value_type>::real_type one_norm_real () const
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
642 typename FieldTraits<value_type>::real_type two_norm () const
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
651 typename FieldTraits<value_type>::real_type two_norm2 () const
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>
662 typename FieldTraits<vt>::real_type infinity_norm() const {
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>
678 typename FieldTraits<vt>::real_type infinity_norm_real() const {
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>
693 typename FieldTraits<vt>::real_type infinity_norm() const {
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>
711 typename FieldTraits<vt>::real_type infinity_norm_real() const {
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.
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:129
SizeType index() const
return index
Definition: densevector.hh:211
std::ptrdiff_t DifferenceType
The type of the difference between two positions.
Definition: densevector.hh:140
C::size_type SizeType
The type to index the underlying container.
Definition: densevector.hh:145
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
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
ConstIterator const_iterator
typedef for stl compliant access
Definition: densevector.hh:382
Iterator iterator
typedef for stl compliant access
Definition: densevector.hh:345
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
value_type & front()
return reference to first element
Definition: densevector.hh:307
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:642
ConstIterator beforeBegin() const
Definition: densevector.hh:405
bool operator==(const DenseVector< Other > &x) const
Binary vector comparison.
Definition: densevector.hh:556
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
Traits::derived_type derived_type
type of derived vector class
Definition: densevector.hh:245
const value_type & front() const
return reference to first element
Definition: densevector.hh:313
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
size_type dim() const
dimension of the vector space
Definition: densevector.hh:734
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (maximum of absolute values of entries)
Definition: densevector.hh:662
ConstIterator beforeEnd() const
Definition: densevector.hh:398
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition: densevector.hh:576
derived_type & operator=(const value_type &k)
Assignment operator for scalar.
Definition: densevector.hh:267
Iterator end()
end iterator
Definition: densevector.hh:354
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
derived_type & operator-=(const DenseVector< Other > &x)
vector space subtraction
Definition: densevector.hh:430
Iterator beforeEnd()
Definition: densevector.hh:361
derived_type & operator+=(const DenseVector< Other > &x)
vector space addition
Definition: densevector.hh:420
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
bool operator!=(const DenseVector< Other > &x) const
Binary vector incomparison.
Definition: densevector.hh:568
const value_type & back() const
return reference to last element
Definition: densevector.hh:325
ConstIterator begin() const
begin ConstIterator
Definition: densevector.hh:385
PromotionTraits< field_type, typenameDenseVector< Other >::field_type >::PromotedType operator*(const DenseVector< Other > &x) const
indefinite vector dot product which corresponds to Petsc's VecTDot
Definition: densevector.hh:592
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densevector.hh:678
DenseVector & operator=(const DenseVector &)=default
Assignment operator for other DenseVector of same type.
Traits::value_type block_type
export the type representing the components
Definition: densevector.hh:254
value_type & operator[](size_type i)
random access
Definition: densevector.hh:296
@ blocklevel
The number of block levels we contain.
Definition: densevector.hh:262
FieldTraits< value_type >::field_type field_type
export the type representing the field
Definition: densevector.hh:251
value_type & back()
return reference to last element
Definition: densevector.hh:319
derived_type operator-() const
Vector negation.
Definition: densevector.hh:455
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
FieldTraits< value_type >::real_type one_norm_real() const
simplified one norm (uses Manhattan norm for complex values)
Definition: densevector.hh:633
PromotionTraits< field_type, typenameDenseVector< Other >::field_type >::PromotedType dot(const DenseVector< Other > &x) const
vector dot product which corresponds to Petsc's VecDot
Definition: densevector.hh:610
Iterator find(size_type i)
return iterator to given element or end()
Definition: densevector.hh:374
FieldTraits< value_type >::real_type one_norm() const
one norm (sum over absolute values of entries)
Definition: densevector.hh:623
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
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:432
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.
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
auto min(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::min()
Definition: defaults.hh:87
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Documentation of the traits classes you need to write for each implementation of DenseVector or Dense...
Dune namespace.
Definition: alignedallocator.hh:11
STL namespace.
Compute type of the result of an arithmetic operation involving two different number types.
Compute type of the result of an arithmetic operation involving two different number types.
Definition: promotiontraits.hh:25
get the 'mutable' version of a reference to a const object
Definition: genericiterator.hh:114
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)