DUNE PDELab (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// SPDX-FileCopyrightInfo: Copyright © 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>
25 struct FieldTraits< DenseVector<V> >
26 {
27 typedef typename FieldTraits< typename DenseMatVecTraits<V>::value_type >::field_type field_type;
28 typedef typename FieldTraits< typename DenseMatVecTraits<V>::value_type >::real_type real_type;
29 };
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
154 DenseIterator(C& cont, SizeType pos)
155 : container_(&cont), position_(pos)
156 {}
157
158 DenseIterator(const MutableIterator & other)
159 : container_(other.container_), position_(other.position_)
160 {}
161
162 DenseIterator(const ConstIterator & other)
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
192 R elementAt(DifferenceType i) const {
193 return container_->operator[](position_+i);
194 }
195
196 void advance(DifferenceType n){
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 {
230 typedef DenseMatVecTraits<V> Traits;
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
253 typedef typename FieldTraits< value_type >::field_type field_type;
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
300 const value_type & operator[] (size_type i) const
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
622 typename FieldTraits<value_type>::real_type one_norm() const {
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
632 typename FieldTraits<value_type>::real_type one_norm_real () const
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
641 typename FieldTraits<value_type>::real_type two_norm () const
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
650 typename FieldTraits<value_type>::real_type two_norm2 () const
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>
661 typename FieldTraits<vt>::real_type infinity_norm() const {
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>
677 typename FieldTraits<vt>::real_type infinity_norm_real() const {
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>
692 typename FieldTraits<vt>::real_type infinity_norm() const {
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>
710 typename FieldTraits<vt>::real_type infinity_norm_real() const {
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
Macro for wrapping boundary checks.
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:131
SizeType index() const
return index
Definition: densevector.hh:213
std::ptrdiff_t DifferenceType
The type of the difference between two positions.
Definition: densevector.hh:142
C::size_type SizeType
The type to index the underlying container.
Definition: densevector.hh:147
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
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
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
Base class for stl conformant forward iterators.
Definition: iteratorfacades.hh:435
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.
#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
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto min
Function object that returns the smaller of the given values.
Definition: hybridutilities.hh:506
Documentation of the traits classes you need to write for each implementation of DenseVector or Dense...
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is NaN.
Definition: fvector.hh:627
Dune namespace.
Definition: alignedallocator.hh:13
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:27
get the 'mutable' version of a reference to a const object
Definition: genericiterator.hh:116
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)