Dune Core Modules (2.5.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 <limits>
7#include <type_traits>
8
9#include "genericiterator.hh"
10#include "ftraits.hh"
11#include "matvectraits.hh"
12#include "promotiontraits.hh"
13#include "dotproduct.hh"
14#include "boundschecking.hh"
15
16namespace Dune {
17
18 // forward declaration of template
19 template<typename V> class DenseVector;
20
21 template<typename V>
22 struct FieldTraits< DenseVector<V> >
23 {
24 typedef typename FieldTraits< typename DenseMatVecTraits<V>::value_type >::field_type field_type;
25 typedef typename FieldTraits< typename DenseMatVecTraits<V>::value_type >::real_type real_type;
26 };
27
37 namespace fvmeta
38 {
43 template<class K>
44 inline typename FieldTraits<K>::real_type absreal (const K& k)
45 {
46 using std::abs;
47 return abs(k);
48 }
49
54 template<class K>
55 inline typename FieldTraits<K>::real_type absreal (const std::complex<K>& c)
56 {
57 using std::abs;
58 return abs(c.real()) + abs(c.imag());
59 }
60
65 template<class K>
66 inline typename FieldTraits<K>::real_type abs2 (const K& k)
67 {
68 return k*k;
69 }
70
75 template<class K>
76 inline typename FieldTraits<K>::real_type abs2 (const std::complex<K>& c)
77 {
78 return c.real()*c.real() + c.imag()*c.imag();
79 }
80
85 template<class K, bool isInteger = std::numeric_limits<K>::is_integer>
86 struct Sqrt
87 {
88 static inline typename FieldTraits<K>::real_type sqrt (const K& k)
89 {
90 using std::sqrt;
91 return sqrt(k);
92 }
93 };
94
99 template<class K>
100 struct Sqrt<K, true>
101 {
102 static inline typename FieldTraits<K>::real_type sqrt (const K& k)
103 {
104 using std::sqrt;
105 return typename FieldTraits<K>::real_type(sqrt(double(k)));
106 }
107 };
108
113 template<class K>
114 inline typename FieldTraits<K>::real_type sqrt (const K& k)
115 {
116 return Sqrt<K>::sqrt(k);
117 }
118
119 }
120
125 template<class C, class T, class R =T&>
127 public Dune::RandomAccessIteratorFacade<DenseIterator<C,T,R>,T, R, std::ptrdiff_t>
128 {
129 friend class DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type, typename mutable_reference<R>::type >;
130 friend class DenseIterator<const typename std::remove_const<C>::type, const typename std::remove_const<T>::type, typename const_reference<R>::type >;
131
132 typedef DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type, typename mutable_reference<R>::type > MutableIterator;
133 typedef DenseIterator<const typename std::remove_const<C>::type, const typename std::remove_const<T>::type, typename const_reference<R>::type > ConstIterator;
134 public:
135
139 typedef std::ptrdiff_t DifferenceType;
140
144 typedef typename C::size_type SizeType;
145
146 // Constructors needed by the base iterators.
148 : container_(0), position_()
149 {}
150
151 DenseIterator(C& cont, SizeType pos)
152 : container_(&cont), position_(pos)
153 {}
154
155 DenseIterator(const MutableIterator & other)
156 : container_(other.container_), position_(other.position_)
157 {}
158
159 DenseIterator(const ConstIterator & other)
160 : container_(other.container_), position_(other.position_)
161 {}
162
163 // Methods needed by the forward iterator
164 bool equals(const MutableIterator &other) const
165 {
166 return position_ == other.position_ && container_ == other.container_;
167 }
168
169
170 bool equals(const ConstIterator & other) const
171 {
172 return position_ == other.position_ && container_ == other.container_;
173 }
174
175 R dereference() const {
176 return container_->operator[](position_);
177 }
178
179 void increment(){
180 ++position_;
181 }
182
183 // Additional function needed by BidirectionalIterator
184 void decrement(){
185 --position_;
186 }
187
188 // Additional function needed by RandomAccessIterator
189 R elementAt(DifferenceType i) const {
190 return container_->operator[](position_+i);
191 }
192
193 void advance(DifferenceType n){
194 position_=position_+n;
195 }
196
197 DifferenceType distanceTo(DenseIterator<const typename std::remove_const<C>::type,const typename std::remove_const<T>::type> other) const
198 {
199 assert(other.container_==container_);
200 return static_cast< DifferenceType >( other.position_ ) - static_cast< DifferenceType >( position_ );
201 }
202
203 DifferenceType distanceTo(DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type> other) const
204 {
205 assert(other.container_==container_);
206 return static_cast< DifferenceType >( other.position_ ) - static_cast< DifferenceType >( position_ );
207 }
208
211 {
212 return this->position_;
213 }
214
215 private:
216 C *container_;
217 SizeType position_;
218 };
219
233 template<typename V>
235 {
236 typedef DenseMatVecTraits<V> Traits;
237 // typedef typename Traits::value_type K;
238
239 // Curiously recurring template pattern
240 V & asImp() { return static_cast<V&>(*this); }
241 const V & asImp() const { return static_cast<const V&>(*this); }
242
243 // prohibit copying
244 DenseVector ( const DenseVector & );
245
246 protected:
247 // construction allowed to derived classes only
248 constexpr DenseVector () {}
249
250 public:
251 //===== type definitions and constants
252
254 typedef typename Traits::derived_type derived_type;
255
257 typedef typename Traits::value_type value_type;
258
260 typedef typename Traits::value_type field_type;
261
263 typedef typename Traits::value_type block_type;
264
266 typedef typename Traits::size_type size_type;
267
269 enum {
271 blocklevel = 1
272 };
273
274 //===== assignment from scalar
277 {
278 for (size_type i=0; i<size(); i++)
279 asImp()[i] = k;
280 return asImp();
281 }
282
283 //===== access to components
284
287 {
288 return asImp()[i];
289 }
290
291 const value_type & operator[] (size_type i) const
292 {
293 return asImp()[i];
294 }
295
298 {
299 return asImp().size();
300 }
301
306
309 {
310 return Iterator(*this,0);
311 }
312
315 {
316 return Iterator(*this,size());
317 }
318
322 {
323 return Iterator(*this,size()-1);
324 }
325
329 {
330 return Iterator(*this,-1);
331 }
332
335 {
336 return Iterator(*this,std::min(i,size()));
337 }
338
343
346 {
347 return ConstIterator(*this,0);
348 }
349
352 {
353 return ConstIterator(*this,size());
354 }
355
359 {
360 return ConstIterator(*this,size()-1);
361 }
362
366 {
367 return ConstIterator(*this,-1);
368 }
369
372 {
373 return ConstIterator(*this,std::min(i,size()));
374 }
375
376 //===== vector space arithmetic
377
379 template <class Other>
381 {
382 DUNE_ASSERT_BOUNDS(y.size() == size());
383 for (size_type i=0; i<size(); i++)
384 (*this)[i] += y[i];
385 return asImp();
386 }
387
389 template <class Other>
391 {
392 DUNE_ASSERT_BOUNDS(y.size() == size());
393 for (size_type i=0; i<size(); i++)
394 (*this)[i] -= y[i];
395 return asImp();
396 }
397
399 template <class Other>
401 {
402 derived_type z = asImp();
403 return (z+=b);
404 }
405
407 template <class Other>
409 {
410 derived_type z = asImp();
411 return (z-=b);
412 }
413
415
423 template <typename ValueType>
424 typename std::enable_if<
425 std::is_convertible<ValueType, value_type>::value,
427 >::type&
428 operator+= (const ValueType& kk)
429 {
430 const value_type& k = kk;
431 for (size_type i=0; i<size(); i++)
432 (*this)[i] += k;
433 return asImp();
434 }
435
437
445 template <typename ValueType>
446 typename std::enable_if<
447 std::is_convertible<ValueType, value_type>::value,
449 >::type&
450 operator-= (const ValueType& kk)
451 {
452 const value_type& k = kk;
453 for (size_type i=0; i<size(); i++)
454 (*this)[i] -= k;
455 return asImp();
456 }
457
459
467 template <typename ValueType>
468 typename std::enable_if<
469 std::is_convertible<ValueType, value_type>::value,
471 >::type&
472 operator*= (const ValueType& kk)
473 {
474 const value_type& k = kk;
475 for (size_type i=0; i<size(); i++)
476 (*this)[i] *= k;
477 return asImp();
478 }
479
481
489 template <typename ValueType>
490 typename std::enable_if<
491 std::is_convertible<ValueType, value_type>::value,
493 >::type&
494 operator/= (const ValueType& kk)
495 {
496 const value_type& k = kk;
497 for (size_type i=0; i<size(); i++)
498 (*this)[i] /= k;
499 return asImp();
500 }
501
503 template <class Other>
504 bool operator== (const DenseVector<Other>& y) const
505 {
506 DUNE_ASSERT_BOUNDS(y.size() == size());
507 for (size_type i=0; i<size(); i++)
508 if ((*this)[i]!=y[i])
509 return false;
510
511 return true;
512 }
513
515 template <class Other>
516 bool operator!= (const DenseVector<Other>& y) const
517 {
518 return !operator==(y);
519 }
520
521
523 template <class Other>
525 {
526 DUNE_ASSERT_BOUNDS(y.size() == size());
527 for (size_type i=0; i<size(); i++)
528 (*this)[i] += a*y[i];
529 return asImp();
530 }
531
539 template<class Other>
541 typedef typename PromotionTraits<field_type, typename DenseVector<Other>::field_type>::PromotedType PromotedType;
542 PromotedType result(0);
543 assert(y.size() == size());
544 for (size_type i=0; i<size(); i++) {
545 result += PromotedType((*this)[i]*y[i]);
546 }
547 return result;
548 }
549
557 template<class Other>
559 typedef typename PromotionTraits<field_type, typename DenseVector<Other>::field_type>::PromotedType PromotedType;
560 PromotedType result(0);
561 assert(y.size() == size());
562 for (size_type i=0; i<size(); i++) {
563 result += Dune::dot((*this)[i],y[i]);
564 }
565 return result;
566 }
567
568 //===== norms
569
571 typename FieldTraits<value_type>::real_type one_norm() const {
572 using std::abs;
573 typename FieldTraits<value_type>::real_type result( 0 );
574 for (size_type i=0; i<size(); i++)
575 result += abs((*this)[i]);
576 return result;
577 }
578
579
581 typename FieldTraits<value_type>::real_type one_norm_real () const
582 {
583 typename FieldTraits<value_type>::real_type result( 0 );
584 for (size_type i=0; i<size(); i++)
585 result += fvmeta::absreal((*this)[i]);
586 return result;
587 }
588
590 typename FieldTraits<value_type>::real_type two_norm () const
591 {
592 typename FieldTraits<value_type>::real_type result( 0 );
593 for (size_type i=0; i<size(); i++)
594 result += fvmeta::abs2((*this)[i]);
595 return fvmeta::sqrt(result);
596 }
597
599 typename FieldTraits<value_type>::real_type two_norm2 () const
600 {
601 typename FieldTraits<value_type>::real_type result( 0 );
602 for (size_type i=0; i<size(); i++)
603 result += fvmeta::abs2((*this)[i]);
604 return result;
605 }
606
608 template <typename vt = value_type,
609 typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
610 typename FieldTraits<vt>::real_type infinity_norm() const {
611 using real_type = typename FieldTraits<vt>::real_type;
612 using std::abs;
613 using std::max;
614
615 real_type norm = 0;
616 for (auto const &x : *this) {
617 real_type const a = abs(x);
618 norm = max(a, norm);
619 }
620 return norm;
621 }
622
624 template <typename vt = value_type,
625 typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
626 typename FieldTraits<vt>::real_type infinity_norm_real() const {
627 using real_type = typename FieldTraits<vt>::real_type;
628 using std::max;
629
630 real_type norm = 0;
631 for (auto const &x : *this) {
632 real_type const a = fvmeta::absreal(x);
633 norm = max(a, norm);
634 }
635 return norm;
636 }
637
639 template <typename vt = value_type,
640 typename std::enable_if<has_nan<vt>::value, int>::type = 0>
641 typename FieldTraits<vt>::real_type infinity_norm() const {
642 using real_type = typename FieldTraits<vt>::real_type;
643 using std::abs;
644 using std::max;
645
646 real_type norm = 0;
647 real_type isNaN = 1;
648 for (auto const &x : *this) {
649 real_type const a = abs(x);
650 norm = max(a, norm);
651 isNaN += a;
652 }
653 isNaN /= isNaN;
654 return norm * isNaN;
655 }
656
658 template <typename vt = value_type,
659 typename std::enable_if<has_nan<vt>::value, int>::type = 0>
660 typename FieldTraits<vt>::real_type infinity_norm_real() const {
661 using real_type = typename FieldTraits<vt>::real_type;
662 using std::max;
663
664 real_type norm = 0;
665 real_type isNaN = 1;
666 for (auto const &x : *this) {
667 real_type const a = fvmeta::absreal(x);
668 norm = max(a, norm);
669 isNaN += a;
670 }
671 isNaN /= isNaN;
672 return norm * isNaN;
673 }
674
675 //===== sizes
676
678 size_type N () const
679 {
680 return size();
681 }
682
684 size_type dim () const
685 {
686 return size();
687 }
688
689 };
690
699 template<typename V>
700 std::ostream& operator<< (std::ostream& s, const DenseVector<V>& v)
701 {
702 for (typename DenseVector<V>::size_type i=0; i<v.size(); i++)
703 s << ((i>0) ? " " : "") << v[i];
704 return s;
705 }
706
709} // end namespace
710
711#endif // DUNE_DENSEVECTOR_HH
Macro for wrapping boundary checks.
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:128
SizeType index() const
return index
Definition: densevector.hh:210
std::ptrdiff_t DifferenceType
The type of the difference between two positions.
Definition: densevector.hh:139
C::size_type SizeType
The type to index the underlying container.
Definition: densevector.hh:144
Interface for a class of dense vectors over a given field.
Definition: densevector.hh:235
Traits::value_type value_type
export the type representing the field
Definition: densevector.hh:257
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:599
ConstIterator const_iterator
typedef for stl compliant access
Definition: densevector.hh:342
Iterator iterator
typedef for stl compliant access
Definition: densevector.hh:305
ConstIterator find(size_type i) const
return iterator to given element or end()
Definition: densevector.hh:371
ConstIterator end() const
end ConstIterator
Definition: densevector.hh:351
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition: densevector.hh:590
ConstIterator beforeBegin() const
Definition: densevector.hh:365
PromotionTraits< field_type, typenameDenseVector< Other >::field_type >::PromotedType dot(const DenseVector< Other > &y) const
vector dot product which corresponds to Petsc's VecDot
Definition: densevector.hh:558
Iterator begin()
begin iterator
Definition: densevector.hh:308
derived_type & axpy(const value_type &a, const DenseVector< Other > &y)
vector space axpy operation ( *this += a y )
Definition: densevector.hh:524
Iterator beforeBegin()
Definition: densevector.hh:328
DenseIterator< const DenseVector, const value_type > ConstIterator
ConstIterator class for sequential access.
Definition: densevector.hh:340
Traits::derived_type derived_type
type of derived vector class
Definition: densevector.hh:254
derived_type operator+(const DenseVector< Other > &b) const
Binary vector addition.
Definition: densevector.hh:400
@ blocklevel
The number of block levels we contain.
Definition: densevector.hh:271
std::enable_if< std::is_convertible< ValueType, value_type >::value, derived_type >::type & operator/=(const ValueType &kk)
vector space division by scalar
Definition: densevector.hh:494
size_type size() const
size method
Definition: densevector.hh:297
derived_type & operator-=(const DenseVector< Other > &y)
vector space subtraction
Definition: densevector.hh:390
size_type dim() const
dimension of the vector space
Definition: densevector.hh:684
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (maximum of absolute values of entries)
Definition: densevector.hh:610
ConstIterator beforeEnd() const
Definition: densevector.hh:358
derived_type & operator=(const value_type &k)
Assignment operator for scalar.
Definition: densevector.hh:276
Iterator end()
end iterator
Definition: densevector.hh:314
PromotionTraits< field_type, typenameDenseVector< Other >::field_type >::PromotedType operator*(const DenseVector< Other > &y) const
indefinite vector dot product which corresponds to Petsc's VecTDot
Definition: densevector.hh:540
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densevector.hh:266
derived_type operator-(const DenseVector< Other > &b) const
Binary vector subtraction.
Definition: densevector.hh:408
DenseIterator< DenseVector, value_type > Iterator
Iterator class for sequential access.
Definition: densevector.hh:303
bool operator==(const DenseVector< Other > &y) const
Binary vector comparison.
Definition: densevector.hh:504
Iterator beforeEnd()
Definition: densevector.hh:321
std::enable_if< std::is_convertible< ValueType, value_type >::value, derived_type >::type & operator*=(const ValueType &kk)
vector space multiplication with scalar
Definition: densevector.hh:472
derived_type & operator+=(const DenseVector< Other > &y)
vector space addition
Definition: densevector.hh:380
ConstIterator begin() const
begin ConstIterator
Definition: densevector.hh:345
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densevector.hh:626
Traits::value_type block_type
export the type representing the components
Definition: densevector.hh:263
value_type & operator[](size_type i)
random access
Definition: densevector.hh:286
Traits::value_type field_type
export the type representing the field
Definition: densevector.hh:260
bool operator!=(const DenseVector< Other > &y) const
Binary vector incomparison.
Definition: densevector.hh:516
FieldTraits< value_type >::real_type one_norm_real() const
simplified one norm (uses Manhattan norm for complex values)
Definition: densevector.hh:581
Iterator find(size_type i)
return iterator to given element or end()
Definition: densevector.hh:334
FieldTraits< value_type >::real_type one_norm() const
one norm (sum over absolute values of entries)
Definition: densevector.hh:571
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition: densevector.hh:678
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:433
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:43
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition: boundschecking.hh:28
Documentation of the traits classes you need to write for each implementation of DenseVector or Dense...
Dune namespace.
Definition: alignment.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 23, 23:29, 2024)