00001
00002
00003 #ifndef DUNE_DENSEVECTOR_HH
00004 #define DUNE_DENSEVECTOR_HH
00005
00006 #include <limits>
00007 #include <type_traits>
00008
00009 #include "genericiterator.hh"
00010 #include "ftraits.hh"
00011 #include "matvectraits.hh"
00012 #include "promotiontraits.hh"
00013 #include "dotproduct.hh"
00014 #include "boundschecking.hh"
00015
00016 namespace Dune {
00017
00018
00019 template<typename V> class DenseVector;
00020
00021 template<typename V>
00022 struct FieldTraits< DenseVector<V> >
00023 {
00024 typedef typename FieldTraits< typename DenseMatVecTraits<V>::value_type >::field_type field_type;
00025 typedef typename FieldTraits< typename DenseMatVecTraits<V>::value_type >::real_type real_type;
00026 };
00027
00037 namespace fvmeta
00038 {
00043 template<class K>
00044 inline typename FieldTraits<K>::real_type absreal (const K& k)
00045 {
00046 using std::abs;
00047 return abs(k);
00048 }
00049
00054 template<class K>
00055 inline typename FieldTraits<K>::real_type absreal (const std::complex<K>& c)
00056 {
00057 using std::abs;
00058 return abs(c.real()) + abs(c.imag());
00059 }
00060
00065 template<class K>
00066 inline typename FieldTraits<K>::real_type abs2 (const K& k)
00067 {
00068 return k*k;
00069 }
00070
00075 template<class K>
00076 inline typename FieldTraits<K>::real_type abs2 (const std::complex<K>& c)
00077 {
00078 return c.real()*c.real() + c.imag()*c.imag();
00079 }
00080
00085 template<class K, bool isInteger = std::numeric_limits<K>::is_integer>
00086 struct Sqrt
00087 {
00088 static inline typename FieldTraits<K>::real_type sqrt (const K& k)
00089 {
00090 using std::sqrt;
00091 return sqrt(k);
00092 }
00093 };
00094
00099 template<class K>
00100 struct Sqrt<K, true>
00101 {
00102 static inline typename FieldTraits<K>::real_type sqrt (const K& k)
00103 {
00104 using std::sqrt;
00105 return typename FieldTraits<K>::real_type(sqrt(double(k)));
00106 }
00107 };
00108
00113 template<class K>
00114 inline typename FieldTraits<K>::real_type sqrt (const K& k)
00115 {
00116 return Sqrt<K>::sqrt(k);
00117 }
00118
00119 }
00120
00125 template<class C, class T, class R =T&>
00126 class DenseIterator :
00127 public Dune::RandomAccessIteratorFacade<DenseIterator<C,T,R>,T, R, std::ptrdiff_t>
00128 {
00129 friend class DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type, typename mutable_reference<R>::type >;
00130 friend class DenseIterator<const typename std::remove_const<C>::type, const typename std::remove_const<T>::type, typename const_reference<R>::type >;
00131
00132 typedef DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type, typename mutable_reference<R>::type > MutableIterator;
00133 typedef DenseIterator<const typename std::remove_const<C>::type, const typename std::remove_const<T>::type, typename const_reference<R>::type > ConstIterator;
00134 public:
00135
00139 typedef std::ptrdiff_t DifferenceType;
00140
00144 typedef typename C::size_type SizeType;
00145
00146
00147 DenseIterator()
00148 : container_(0), position_()
00149 {}
00150
00151 DenseIterator(C& cont, SizeType pos)
00152 : container_(&cont), position_(pos)
00153 {}
00154
00155 DenseIterator(const MutableIterator & other)
00156 : container_(other.container_), position_(other.position_)
00157 {}
00158
00159 DenseIterator(const ConstIterator & other)
00160 : container_(other.container_), position_(other.position_)
00161 {}
00162
00163
00164 bool equals(const MutableIterator &other) const
00165 {
00166 return position_ == other.position_ && container_ == other.container_;
00167 }
00168
00169
00170 bool equals(const ConstIterator & other) const
00171 {
00172 return position_ == other.position_ && container_ == other.container_;
00173 }
00174
00175 R dereference() const {
00176 return container_->operator[](position_);
00177 }
00178
00179 void increment(){
00180 ++position_;
00181 }
00182
00183
00184 void decrement(){
00185 --position_;
00186 }
00187
00188
00189 R elementAt(DifferenceType i) const {
00190 return container_->operator[](position_+i);
00191 }
00192
00193 void advance(DifferenceType n){
00194 position_=position_+n;
00195 }
00196
00197 DifferenceType distanceTo(DenseIterator<const typename std::remove_const<C>::type,const typename std::remove_const<T>::type> other) const
00198 {
00199 assert(other.container_==container_);
00200 return static_cast< DifferenceType >( other.position_ ) - static_cast< DifferenceType >( position_ );
00201 }
00202
00203 DifferenceType distanceTo(DenseIterator<typename std::remove_const<C>::type, typename std::remove_const<T>::type> other) const
00204 {
00205 assert(other.container_==container_);
00206 return static_cast< DifferenceType >( other.position_ ) - static_cast< DifferenceType >( position_ );
00207 }
00208
00210 SizeType index () const
00211 {
00212 return this->position_;
00213 }
00214
00215 private:
00216 C *container_;
00217 SizeType position_;
00218 };
00219
00233 template<typename V>
00234 class DenseVector
00235 {
00236 typedef DenseMatVecTraits<V> Traits;
00237
00238
00239
00240 V & asImp() { return static_cast<V&>(*this); }
00241 const V & asImp() const { return static_cast<const V&>(*this); }
00242
00243
00244 DenseVector ( const DenseVector & );
00245
00246 protected:
00247
00248 constexpr DenseVector () {}
00249
00250 public:
00251
00252
00254 typedef typename Traits::derived_type derived_type;
00255
00257 typedef typename Traits::value_type value_type;
00258
00260 typedef typename FieldTraits< value_type >::field_type field_type;
00261
00263 typedef typename Traits::value_type block_type;
00264
00266 typedef typename Traits::size_type size_type;
00267
00269 enum {
00271 blocklevel = 1
00272 };
00273
00274
00276 inline derived_type& operator= (const value_type& k)
00277 {
00278 for (size_type i=0; i<size(); i++)
00279 asImp()[i] = k;
00280 return asImp();
00281 }
00282
00283
00284
00286 value_type & operator[] (size_type i)
00287 {
00288 return asImp()[i];
00289 }
00290
00291 const value_type & operator[] (size_type i) const
00292 {
00293 return asImp()[i];
00294 }
00295
00297 size_type size() const
00298 {
00299 return asImp().size();
00300 }
00301
00303 typedef DenseIterator<DenseVector,value_type> Iterator;
00305 typedef Iterator iterator;
00306
00308 Iterator begin ()
00309 {
00310 return Iterator(*this,0);
00311 }
00312
00314 Iterator end ()
00315 {
00316 return Iterator(*this,size());
00317 }
00318
00321 Iterator beforeEnd ()
00322 {
00323 return Iterator(*this,size()-1);
00324 }
00325
00328 Iterator beforeBegin ()
00329 {
00330 return Iterator(*this,-1);
00331 }
00332
00334 Iterator find (size_type i)
00335 {
00336 return Iterator(*this,std::min(i,size()));
00337 }
00338
00340 typedef DenseIterator<const DenseVector,const value_type> ConstIterator;
00342 typedef ConstIterator const_iterator;
00343
00345 ConstIterator begin () const
00346 {
00347 return ConstIterator(*this,0);
00348 }
00349
00351 ConstIterator end () const
00352 {
00353 return ConstIterator(*this,size());
00354 }
00355
00358 ConstIterator beforeEnd () const
00359 {
00360 return ConstIterator(*this,size()-1);
00361 }
00362
00365 ConstIterator beforeBegin () const
00366 {
00367 return ConstIterator(*this,-1);
00368 }
00369
00371 ConstIterator find (size_type i) const
00372 {
00373 return ConstIterator(*this,std::min(i,size()));
00374 }
00375
00376
00377
00379 template <class Other>
00380 derived_type& operator+= (const DenseVector<Other>& y)
00381 {
00382 DUNE_ASSERT_BOUNDS(y.size() == size());
00383 for (size_type i=0; i<size(); i++)
00384 (*this)[i] += y[i];
00385 return asImp();
00386 }
00387
00389 template <class Other>
00390 derived_type& operator-= (const DenseVector<Other>& y)
00391 {
00392 DUNE_ASSERT_BOUNDS(y.size() == size());
00393 for (size_type i=0; i<size(); i++)
00394 (*this)[i] -= y[i];
00395 return asImp();
00396 }
00397
00399 template <class Other>
00400 derived_type operator+ (const DenseVector<Other>& b) const
00401 {
00402 derived_type z = asImp();
00403 return (z+=b);
00404 }
00405
00407 template <class Other>
00408 derived_type operator- (const DenseVector<Other>& b) const
00409 {
00410 derived_type z = asImp();
00411 return (z-=b);
00412 }
00413
00415
00423 template <typename ValueType>
00424 typename std::enable_if<
00425 std::is_convertible<ValueType, value_type>::value,
00426 derived_type
00427 >::type&
00428 operator+= (const ValueType& kk)
00429 {
00430 const value_type& k = kk;
00431 for (size_type i=0; i<size(); i++)
00432 (*this)[i] += k;
00433 return asImp();
00434 }
00435
00437
00445 template <typename ValueType>
00446 typename std::enable_if<
00447 std::is_convertible<ValueType, value_type>::value,
00448 derived_type
00449 >::type&
00450 operator-= (const ValueType& kk)
00451 {
00452 const value_type& k = kk;
00453 for (size_type i=0; i<size(); i++)
00454 (*this)[i] -= k;
00455 return asImp();
00456 }
00457
00459
00467 template <typename FieldType>
00468 typename std::enable_if<
00469 std::is_convertible<FieldType, field_type>::value,
00470 derived_type
00471 >::type&
00472 operator*= (const FieldType& kk)
00473 {
00474 const field_type& k = kk;
00475 for (size_type i=0; i<size(); i++)
00476 (*this)[i] *= k;
00477 return asImp();
00478 }
00479
00481
00489 template <typename FieldType>
00490 typename std::enable_if<
00491 std::is_convertible<FieldType, field_type>::value,
00492 derived_type
00493 >::type&
00494 operator/= (const FieldType& kk)
00495 {
00496 const field_type& k = kk;
00497 for (size_type i=0; i<size(); i++)
00498 (*this)[i] /= k;
00499 return asImp();
00500 }
00501
00503 template <class Other>
00504 bool operator== (const DenseVector<Other>& y) const
00505 {
00506 DUNE_ASSERT_BOUNDS(y.size() == size());
00507 for (size_type i=0; i<size(); i++)
00508 if ((*this)[i]!=y[i])
00509 return false;
00510
00511 return true;
00512 }
00513
00515 template <class Other>
00516 bool operator!= (const DenseVector<Other>& y) const
00517 {
00518 return !operator==(y);
00519 }
00520
00521
00523 template <class Other>
00524 derived_type& axpy (const field_type& a, const DenseVector<Other>& y)
00525 {
00526 DUNE_ASSERT_BOUNDS(y.size() == size());
00527 for (size_type i=0; i<size(); i++)
00528 (*this)[i] += a*y[i];
00529 return asImp();
00530 }
00531
00539 template<class Other>
00540 typename PromotionTraits<field_type,typename DenseVector<Other>::field_type>::PromotedType operator* (const DenseVector<Other>& y) const {
00541 typedef typename PromotionTraits<field_type, typename DenseVector<Other>::field_type>::PromotedType PromotedType;
00542 PromotedType result(0);
00543 assert(y.size() == size());
00544 for (size_type i=0; i<size(); i++) {
00545 result += PromotedType((*this)[i]*y[i]);
00546 }
00547 return result;
00548 }
00549
00557 template<class Other>
00558 typename PromotionTraits<field_type,typename DenseVector<Other>::field_type>::PromotedType dot(const DenseVector<Other>& y) const {
00559 typedef typename PromotionTraits<field_type, typename DenseVector<Other>::field_type>::PromotedType PromotedType;
00560 PromotedType result(0);
00561 assert(y.size() == size());
00562 for (size_type i=0; i<size(); i++) {
00563 result += Dune::dot((*this)[i],y[i]);
00564 }
00565 return result;
00566 }
00567
00568
00569
00571 typename FieldTraits<value_type>::real_type one_norm() const {
00572 using std::abs;
00573 typename FieldTraits<value_type>::real_type result( 0 );
00574 for (size_type i=0; i<size(); i++)
00575 result += abs((*this)[i]);
00576 return result;
00577 }
00578
00579
00581 typename FieldTraits<value_type>::real_type one_norm_real () const
00582 {
00583 typename FieldTraits<value_type>::real_type result( 0 );
00584 for (size_type i=0; i<size(); i++)
00585 result += fvmeta::absreal((*this)[i]);
00586 return result;
00587 }
00588
00590 typename FieldTraits<value_type>::real_type two_norm () const
00591 {
00592 typename FieldTraits<value_type>::real_type result( 0 );
00593 for (size_type i=0; i<size(); i++)
00594 result += fvmeta::abs2((*this)[i]);
00595 return fvmeta::sqrt(result);
00596 }
00597
00599 typename FieldTraits<value_type>::real_type two_norm2 () const
00600 {
00601 typename FieldTraits<value_type>::real_type result( 0 );
00602 for (size_type i=0; i<size(); i++)
00603 result += fvmeta::abs2((*this)[i]);
00604 return result;
00605 }
00606
00608 template <typename vt = value_type,
00609 typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
00610 typename FieldTraits<vt>::real_type infinity_norm() const {
00611 using real_type = typename FieldTraits<vt>::real_type;
00612 using std::abs;
00613 using std::max;
00614
00615 real_type norm = 0;
00616 for (auto const &x : *this) {
00617 real_type const a = abs(x);
00618 norm = max(a, norm);
00619 }
00620 return norm;
00621 }
00622
00624 template <typename vt = value_type,
00625 typename std::enable_if<!has_nan<vt>::value, int>::type = 0>
00626 typename FieldTraits<vt>::real_type infinity_norm_real() const {
00627 using real_type = typename FieldTraits<vt>::real_type;
00628 using std::max;
00629
00630 real_type norm = 0;
00631 for (auto const &x : *this) {
00632 real_type const a = fvmeta::absreal(x);
00633 norm = max(a, norm);
00634 }
00635 return norm;
00636 }
00637
00639 template <typename vt = value_type,
00640 typename std::enable_if<has_nan<vt>::value, int>::type = 0>
00641 typename FieldTraits<vt>::real_type infinity_norm() const {
00642 using real_type = typename FieldTraits<vt>::real_type;
00643 using std::abs;
00644 using std::max;
00645
00646 real_type norm = 0;
00647 real_type isNaN = 1;
00648 for (auto const &x : *this) {
00649 real_type const a = abs(x);
00650 norm = max(a, norm);
00651 isNaN += a;
00652 }
00653 isNaN /= isNaN;
00654 return norm * isNaN;
00655 }
00656
00658 template <typename vt = value_type,
00659 typename std::enable_if<has_nan<vt>::value, int>::type = 0>
00660 typename FieldTraits<vt>::real_type infinity_norm_real() const {
00661 using real_type = typename FieldTraits<vt>::real_type;
00662 using std::max;
00663
00664 real_type norm = 0;
00665 real_type isNaN = 1;
00666 for (auto const &x : *this) {
00667 real_type const a = fvmeta::absreal(x);
00668 norm = max(a, norm);
00669 isNaN += a;
00670 }
00671 isNaN /= isNaN;
00672 return norm * isNaN;
00673 }
00674
00675
00676
00678 size_type N () const
00679 {
00680 return size();
00681 }
00682
00684 size_type dim () const
00685 {
00686 return size();
00687 }
00688
00689 };
00690
00699 template<typename V>
00700 std::ostream& operator<< (std::ostream& s, const DenseVector<V>& v)
00701 {
00702 for (typename DenseVector<V>::size_type i=0; i<v.size(); i++)
00703 s << ((i>0) ? " " : "") << v[i];
00704 return s;
00705 }
00706
00709 }
00710
00711 #endif // DUNE_DENSEVECTOR_HH