dune-common  2.2.0
densevector.hh
Go to the documentation of this file.
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