dune-common 2.9.0
Loading...
Searching...
No Matches
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
19namespace Dune {
20
21 // forward declaration of template
22 template<typename V> class DenseVector;
23
24 template<typename V>
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
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
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:
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
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
Type traits to determine the type of reals (when working with complex numbers)
Compute type of the result of an arithmetic operation involving two different number types.
Provides the functions dot(a,b) := and dotT(a,b) := .
Macro for wrapping boundary checks.
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...
#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
STL namespace.
Dune namespace.
Definition alignedallocator.hh:13
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
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
ConstIterator const_iterator
typedef for stl compliant access
Definition densevector.hh:381
Iterator iterator
typedef for stl compliant access
Definition densevector.hh:344
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
value_type & front()
return reference to first element
Definition densevector.hh:306
FieldTraits< value_type >::real_type two_norm() const
two norm sqrt(sum over squared values of entries)
Definition densevector.hh:641
ConstIterator beforeBegin() const
Definition densevector.hh:404
bool operator==(const DenseVector< Other > &x) const
Binary vector comparison.
Definition densevector.hh:555
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
Traits::derived_type derived_type
type of derived vector class
Definition densevector.hh:247
const value_type & front() const
return reference to first element
Definition densevector.hh:312
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
size_type dim() const
dimension of the vector space
Definition densevector.hh:733
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (maximum of absolute values of entries)
Definition densevector.hh:661
ConstIterator beforeEnd() const
Definition densevector.hh:397
derived_type & axpy(const field_type &a, const DenseVector< Other > &x)
vector space axpy operation ( *this += a x )
Definition densevector.hh:575
derived_type & operator=(const value_type &k)
Assignment operator for scalar.
Definition densevector.hh:266
static constexpr int blocklevel
The number of block levels we contain. This is the leaf, that is, 1.
Definition densevector.hh:262
Iterator end()
end iterator
Definition densevector.hh:353
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
derived_type & operator-=(const DenseVector< Other > &x)
vector space subtraction
Definition densevector.hh:429
DenseVector(const DenseVector &)=default
Iterator beforeEnd()
Definition densevector.hh:360
derived_type & operator+=(const DenseVector< Other > &x)
vector space addition
Definition densevector.hh:419
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
bool operator!=(const DenseVector< Other > &x) const
Binary vector incomparison.
Definition densevector.hh:567
const value_type & back() const
return reference to last element
Definition densevector.hh:324
ConstIterator begin() const
begin ConstIterator
Definition densevector.hh:384
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:591
constexpr DenseVector()=default
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition densevector.hh:677
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:256
value_type & operator[](size_type i)
random access
Definition densevector.hh:295
FieldTraits< value_type >::field_type field_type
export the type representing the field
Definition densevector.hh:253
value_type & back()
return reference to last element
Definition densevector.hh:318
derived_type operator-() const
Vector negation.
Definition densevector.hh:454
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
FieldTraits< value_type >::real_type one_norm_real() const
simplified one norm (uses Manhattan norm for complex values)
Definition densevector.hh:632
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:609
Iterator find(size_type i)
return iterator to given element or end()
Definition densevector.hh:373
FieldTraits< value_type >::real_type one_norm() const
one norm (sum over absolute values of entries)
Definition densevector.hh:622
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< typenameDenseMatVecTraits< V >::value_type >::real_type real_type
Definition densevector.hh:28
FieldTraits< typenameDenseMatVecTraits< V >::value_type >::field_type field_type
Definition densevector.hh:27
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