dune-common
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=4 sw=2 sts=2: 00003 #ifndef DUNE_DENSEVECTOR_HH 00004 #define DUNE_DENSEVECTOR_HH 00005 00006 #include <limits> 00007 00008 #include "genericiterator.hh" 00009 #include "ftraits.hh" 00010 #include "matvectraits.hh" 00011 00012 00013 namespace Dune { 00014 00015 // forward declaration of template 00016 template<typename V> class DenseVector; 00017 00018 template<typename V> 00019 struct FieldTraits< DenseVector<V> > 00020 { 00021 typedef typename FieldTraits< typename DenseMatVecTraits<V>::value_type >::field_type field_type; 00022 typedef typename FieldTraits< typename DenseMatVecTraits<V>::value_type >::real_type real_type; 00023 }; 00024 00034 namespace fvmeta 00035 { 00040 template<class K> 00041 inline typename FieldTraits<K>::real_type absreal (const K& k) 00042 { 00043 return std::abs(k); 00044 } 00045 00050 template<class K> 00051 inline typename FieldTraits<K>::real_type absreal (const std::complex<K>& c) 00052 { 00053 return std::abs(c.real()) + std::abs(c.imag()); 00054 } 00055 00060 template<class K> 00061 inline typename FieldTraits<K>::real_type abs2 (const K& k) 00062 { 00063 return k*k; 00064 } 00065 00070 template<class K> 00071 inline typename FieldTraits<K>::real_type abs2 (const std::complex<K>& c) 00072 { 00073 return c.real()*c.real() + c.imag()*c.imag(); 00074 } 00075 00080 template<class K, bool isInteger = std::numeric_limits<K>::is_integer> 00081 struct Sqrt 00082 { 00083 static inline typename FieldTraits<K>::real_type sqrt (const K& k) 00084 { 00085 return std::sqrt(k); 00086 } 00087 }; 00088 00093 template<class K> 00094 struct Sqrt<K, true> 00095 { 00096 static inline typename FieldTraits<K>::real_type sqrt (const K& k) 00097 { 00098 return typename FieldTraits<K>::real_type(std::sqrt(double(k))); 00099 } 00100 }; 00101 00106 template<class K> 00107 inline typename FieldTraits<K>::real_type sqrt (const K& k) 00108 { 00109 return Sqrt<K>::sqrt(k); 00110 } 00111 00112 } 00113 00118 template<class C, class T> 00119 class DenseIterator : 00120 public Dune::RandomAccessIteratorFacade<DenseIterator<C,T>,T, T&, std::ptrdiff_t> 00121 { 00122 friend class DenseIterator<typename remove_const<C>::type, typename remove_const<T>::type >; 00123 friend class DenseIterator<const typename remove_const<C>::type, const typename remove_const<T>::type >; 00124 00125 public: 00126 00130 typedef std::ptrdiff_t DifferenceType; 00131 00135 typedef typename C::size_type SizeType; 00136 00137 // Constructors needed by the base iterators. 00138 DenseIterator() 00139 : container_(0), position_() 00140 {} 00141 00142 DenseIterator(C& cont, SizeType pos) 00143 : container_(&cont), position_(pos) 00144 {} 00145 00146 DenseIterator(const DenseIterator<typename remove_const<C>::type, typename remove_const<T>::type >& other) 00147 : container_(other.container_), position_(other.position_) 00148 {} 00149 00150 // Methods needed by the forward iterator 00151 bool equals(const DenseIterator<typename remove_const<C>::type,typename remove_const<T>::type>& other) const 00152 { 00153 return position_ == other.position_ && container_ == other.container_; 00154 } 00155 00156 00157 bool equals(const DenseIterator<const typename remove_const<C>::type,const typename remove_const<T>::type>& other) const 00158 { 00159 return position_ == other.position_ && container_ == other.container_; 00160 } 00161 00162 T& dereference() const{ 00163 return container_->operator[](position_); 00164 } 00165 00166 void increment(){ 00167 ++position_; 00168 } 00169 00170 // Additional function needed by BidirectionalIterator 00171 void decrement(){ 00172 --position_; 00173 } 00174 00175 // Additional function needed by RandomAccessIterator 00176 T& elementAt(DifferenceType i)const{ 00177 return container_->operator[](position_+i); 00178 } 00179 00180 void advance(DifferenceType n){ 00181 position_=position_+n; 00182 } 00183 00184 DifferenceType distanceTo(DenseIterator<const typename remove_const<C>::type,const typename remove_const<T>::type> other)const 00185 { 00186 assert(other.container_==container_); 00187 return other.position_ - position_; 00188 } 00189 00190 DifferenceType distanceTo(DenseIterator<typename remove_const<C>::type, typename remove_const<T>::type> other)const 00191 { 00192 assert(other.container_==container_); 00193 return other.position_ - position_; 00194 } 00195 00197 SizeType index () const 00198 { 00199 return this->position_; 00200 } 00201 00202 private: 00203 C *container_; 00204 SizeType position_; 00205 }; 00206 00220 template<typename V> 00221 class DenseVector 00222 { 00223 typedef DenseMatVecTraits<V> Traits; 00224 // typedef typename Traits::value_type K; 00225 00226 // Curiously recurring template pattern 00227 V & asImp() { return static_cast<V&>(*this); } 00228 const V & asImp() const { return static_cast<const V&>(*this); } 00229 00230 public: 00231 //===== type definitions and constants 00232 00234 typedef typename Traits::derived_type derived_type; 00235 00237 typedef typename Traits::value_type value_type; 00238 00240 typedef typename Traits::value_type field_type; 00241 00243 typedef typename Traits::value_type block_type; 00244 00246 typedef typename Traits::size_type size_type; 00247 00249 enum { 00251 blocklevel = 1 00252 }; 00253 00254 //===== assignment from scalar 00256 inline derived_type& operator= (const value_type& k) 00257 { 00258 for (size_type i=0; i<size(); i++) 00259 asImp()[i] = k; 00260 return asImp(); 00261 } 00262 00263 //===== access to components 00264 00266 value_type & operator[] (size_type i) 00267 { 00268 return asImp().vec_access(i); 00269 } 00270 00271 const value_type & operator[] (size_type i) const 00272 { 00273 return asImp().vec_access(i); 00274 } 00275 00277 size_type size() const 00278 { 00279 return asImp().vec_size(); 00280 } 00281 00283 typedef DenseIterator<DenseVector,value_type> Iterator; 00285 typedef Iterator iterator; 00286 00288 Iterator begin () 00289 { 00290 return Iterator(*this,0); 00291 } 00292 00294 Iterator end () 00295 { 00296 return Iterator(*this,size()); 00297 } 00298 00301 Iterator beforeEnd () 00302 { 00303 return Iterator(*this,size()-1); 00304 } 00305 00308 Iterator beforeBegin () 00309 { 00310 return Iterator(*this,-1); 00311 } 00312 00314 Iterator find (size_type i) 00315 { 00316 return Iterator(*this,std::min(i,size())); 00317 } 00318 00320 typedef DenseIterator<const DenseVector,const value_type> ConstIterator; 00322 typedef ConstIterator const_iterator; 00323 00325 ConstIterator begin () const 00326 { 00327 return ConstIterator(*this,0); 00328 } 00329 00331 ConstIterator end () const 00332 { 00333 return ConstIterator(*this,size()); 00334 } 00335 00338 ConstIterator beforeEnd () const 00339 { 00340 return ConstIterator(*this,size()-1); 00341 } 00342 00345 ConstIterator beforeBegin () const 00346 { 00347 return ConstIterator(*this,-1); 00348 } 00349 00351 ConstIterator find (size_type i) const 00352 { 00353 return ConstIterator(*this,std::min(i,size())); 00354 } 00355 00356 //===== vector space arithmetic 00357 00359 template <class Other> 00360 derived_type& operator+= (const DenseVector<Other>& y) 00361 { 00362 assert(y.size() == size()); 00363 for (size_type i=0; i<size(); i++) 00364 (*this)[i] += y[i]; 00365 return asImp(); 00366 } 00367 00369 template <class Other> 00370 derived_type& operator-= (const DenseVector<Other>& y) 00371 { 00372 assert(y.size() == size()); 00373 for (size_type i=0; i<size(); i++) 00374 (*this)[i] -= y[i]; 00375 return asImp(); 00376 } 00377 00379 template <class Other> 00380 derived_type operator+ (const DenseVector<Other>& b) const 00381 { 00382 derived_type z = asImp(); 00383 return (z+=b); 00384 } 00385 00387 template <class Other> 00388 derived_type operator- (const DenseVector<Other>& b) const 00389 { 00390 derived_type z = asImp(); 00391 return (z-=b); 00392 } 00393 00395 derived_type& operator+= (const value_type& k) 00396 { 00397 for (size_type i=0; i<size(); i++) 00398 (*this)[i] += k; 00399 return asImp(); 00400 } 00401 00403 derived_type& operator-= (const value_type& k) 00404 { 00405 for (size_type i=0; i<size(); i++) 00406 (*this)[i] -= k; 00407 return asImp(); 00408 } 00409 00411 derived_type& operator*= (const value_type& k) 00412 { 00413 for (size_type i=0; i<size(); i++) 00414 (*this)[i] *= k; 00415 return asImp(); 00416 } 00417 00419 derived_type& operator/= (const value_type& k) 00420 { 00421 for (size_type i=0; i<size(); i++) 00422 (*this)[i] /= k; 00423 return asImp(); 00424 } 00425 00427 template <class Other> 00428 bool operator== (const DenseVector<Other>& y) const 00429 { 00430 assert(y.size() == size()); 00431 for (size_type i=0; i<size(); i++) 00432 if ((*this)[i]!=y[i]) 00433 return false; 00434 00435 return true; 00436 } 00437 00439 template <class Other> 00440 bool operator!= (const DenseVector<Other>& y) const 00441 { 00442 return !operator==(y); 00443 } 00444 00445 00447 template <class Other> 00448 derived_type& axpy (const value_type& a, const DenseVector<Other>& y) 00449 { 00450 assert(y.size() == size()); 00451 for (size_type i=0; i<size(); i++) 00452 (*this)[i] += a*y[i]; 00453 return asImp(); 00454 } 00455 00456 //===== Euclidean scalar product 00457 00459 template <class Other> 00460 value_type operator* (const DenseVector<Other>& y) const 00461 { 00462 assert(y.size() == size()); 00463 value_type result( 0 ); 00464 for (size_type i=0; i<size(); i++) 00465 result += (*this)[i]*y[i]; 00466 return result; 00467 } 00468 00469 //===== norms 00470 00472 typename FieldTraits<value_type>::real_type one_norm() const { 00473 typename FieldTraits<value_type>::real_type result( 0 ); 00474 for (size_type i=0; i<size(); i++) 00475 result += std::abs((*this)[i]); 00476 return result; 00477 } 00478 00479 00481 typename FieldTraits<value_type>::real_type one_norm_real () const 00482 { 00483 typename FieldTraits<value_type>::real_type result( 0 ); 00484 for (size_type i=0; i<size(); i++) 00485 result += fvmeta::absreal((*this)[i]); 00486 return result; 00487 } 00488 00490 typename FieldTraits<value_type>::real_type two_norm () const 00491 { 00492 typename FieldTraits<value_type>::real_type result( 0 ); 00493 for (size_type i=0; i<size(); i++) 00494 result += fvmeta::abs2((*this)[i]); 00495 return fvmeta::sqrt(result); 00496 } 00497 00499 typename FieldTraits<value_type>::real_type two_norm2 () const 00500 { 00501 typename FieldTraits<value_type>::real_type result( 0 ); 00502 for (size_type i=0; i<size(); i++) 00503 result += fvmeta::abs2((*this)[i]); 00504 return result; 00505 } 00506 00508 typename FieldTraits<value_type>::real_type infinity_norm () const 00509 { 00510 typename FieldTraits<value_type>::real_type result( 0 ); 00511 for (size_type i=0; i<size(); i++) 00512 result = std::max(result, std::abs((*this)[i])); 00513 return result; 00514 } 00515 00517 typename FieldTraits<value_type>::real_type infinity_norm_real () const 00518 { 00519 typename FieldTraits<value_type>::real_type result( 0 ); 00520 for (size_type i=0; i<size(); i++) 00521 result = std::max(result, fvmeta::absreal((*this)[i])); 00522 return result; 00523 } 00524 00525 //===== sizes 00526 00528 size_type N () const 00529 { 00530 return size(); 00531 } 00532 00534 size_type dim () const 00535 { 00536 return size(); 00537 } 00538 00539 }; 00540 00549 template<typename V> 00550 std::ostream& operator<< (std::ostream& s, const DenseVector<V>& v) 00551 { 00552 for (typename DenseVector<V>::size_type i=0; i<v.size(); i++) 00553 s << ((i>0) ? " " : "") << v[i]; 00554 return s; 00555 } 00556 00559 } // end namespace 00560 00561 #endif // DUNE_DENSEVECTOR_HH