dune-common 3.0-git
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>
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
152 : container_(&cont), position_(pos)
153 {}
154
156 : container_(other.container_), position_(other.position_)
157 {}
158
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
190 return container_->operator[](position_+i);
191 }
192
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 {
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
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
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
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
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
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>
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>
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>
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>
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
Implements a generic iterator class for writing stl conformant iterators.
#define DUNE_ASSERT_BOUNDS(cond)
Definition boundschecking.hh:20
Compute type of the result of an arithmetic operation involving two different number types.
Provides the functions dot(a,b) := and dotT(a,b) := .
Documentation of the traits classes you need to write for each implementation of DenseVector or Dense...
Type traits to determine the type of reals (when working with complex numbers)
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
std::ostream & operator<<(std::ostream &s, const std::array< T, N > &e)
Output operator for array.
Definition array.hh:28
STL namespace.
Dune namespace.
Definition alignment.hh:11
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
constexpr DenseVector()
Definition densevector.hh:248
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
FieldTraits< typenameDenseMatVecTraits< V >::value_type >::real_type real_type
Definition densevector.hh:25
FieldTraits< typenameDenseMatVecTraits< V >::value_type >::field_type field_type
Definition densevector.hh:24
Generic iterator class for dense vector and matrix implementations.
Definition densevector.hh:128
void increment()
Definition densevector.hh:179
SizeType index() const
return index
Definition densevector.hh:210
bool equals(const MutableIterator &other) const
Definition densevector.hh:164
DenseIterator(const MutableIterator &other)
Definition densevector.hh:155
bool equals(const ConstIterator &other) const
Definition densevector.hh:170
R elementAt(DifferenceType i) const
Definition densevector.hh:189
DifferenceType distanceTo(DenseIterator< const typename std::remove_const< C >::type, const typename std::remove_const< T >::type > other) const
Definition densevector.hh:197
void decrement()
Definition densevector.hh:184
DenseIterator(const ConstIterator &other)
Definition densevector.hh:159
DifferenceType distanceTo(DenseIterator< typename std::remove_const< C >::type, typename std::remove_const< T >::type > other) const
Definition densevector.hh:203
DenseIterator(C &cont, SizeType pos)
Definition densevector.hh:151
std::ptrdiff_t DifferenceType
The type of the difference between two positions.
Definition densevector.hh:139
R dereference() const
Definition densevector.hh:175
void advance(DifferenceType n)
Definition densevector.hh:193
C::size_type SizeType
The type to index the underlying container.
Definition densevector.hh:144
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:433
Definition matvectraits.hh:29
Compute type of the result of an arithmetic operation involving two different number types.
Definition promotiontraits.hh:25
Definition typetraits.hh:74