dune-common  2.2.1
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 #ifndef DUNE_DENSEVECTOR_HH
4 #define DUNE_DENSEVECTOR_HH
5 
6 #include <limits>
7 
8 #include "genericiterator.hh"
9 #include "ftraits.hh"
10 #include "matvectraits.hh"
11 
12 
13 namespace Dune {
14 
15  // forward declaration of template
16  template<typename V> class DenseVector;
17 
18  template<typename V>
19  struct FieldTraits< DenseVector<V> >
20  {
23  };
24 
34  namespace fvmeta
35  {
40  template<class K>
41  inline typename FieldTraits<K>::real_type absreal (const K& k)
42  {
43  return std::abs(k);
44  }
45 
50  template<class K>
51  inline typename FieldTraits<K>::real_type absreal (const std::complex<K>& c)
52  {
53  return std::abs(c.real()) + std::abs(c.imag());
54  }
55 
60  template<class K>
61  inline typename FieldTraits<K>::real_type abs2 (const K& k)
62  {
63  return k*k;
64  }
65 
70  template<class K>
71  inline typename FieldTraits<K>::real_type abs2 (const std::complex<K>& c)
72  {
73  return c.real()*c.real() + c.imag()*c.imag();
74  }
75 
80  template<class K, bool isInteger = std::numeric_limits<K>::is_integer>
81  struct Sqrt
82  {
83  static inline typename FieldTraits<K>::real_type sqrt (const K& k)
84  {
85  return std::sqrt(k);
86  }
87  };
88 
93  template<class K>
94  struct Sqrt<K, true>
95  {
96  static inline typename FieldTraits<K>::real_type sqrt (const K& k)
97  {
98  return typename FieldTraits<K>::real_type(std::sqrt(double(k)));
99  }
100  };
101 
106  template<class K>
107  inline typename FieldTraits<K>::real_type sqrt (const K& k)
108  {
109  return Sqrt<K>::sqrt(k);
110  }
111 
112  }
113 
118  template<class C, class T>
119  class DenseIterator :
120  public Dune::RandomAccessIteratorFacade<DenseIterator<C,T>,T, T&, std::ptrdiff_t>
121  {
122  friend class DenseIterator<typename remove_const<C>::type, typename remove_const<T>::type >;
123  friend class DenseIterator<const typename remove_const<C>::type, const typename remove_const<T>::type >;
124 
125  public:
126 
130  typedef std::ptrdiff_t DifferenceType;
131 
135  typedef typename C::size_type SizeType;
136 
137  // Constructors needed by the base iterators.
139  : container_(0), position_()
140  {}
141 
142  DenseIterator(C& cont, SizeType pos)
143  : container_(&cont), position_(pos)
144  {}
145 
147  : container_(other.container_), position_(other.position_)
148  {}
149 
150  // Methods needed by the forward iterator
151  bool equals(const DenseIterator<typename remove_const<C>::type,typename remove_const<T>::type>& other) const
152  {
153  return position_ == other.position_ && container_ == other.container_;
154  }
155 
156 
157  bool equals(const DenseIterator<const typename remove_const<C>::type,const typename remove_const<T>::type>& other) const
158  {
159  return position_ == other.position_ && container_ == other.container_;
160  }
161 
162  T& dereference() const{
163  return container_->operator[](position_);
164  }
165 
166  void increment(){
167  ++position_;
168  }
169 
170  // Additional function needed by BidirectionalIterator
171  void decrement(){
172  --position_;
173  }
174 
175  // Additional function needed by RandomAccessIterator
177  return container_->operator[](position_+i);
178  }
179 
181  position_=position_+n;
182  }
183 
185  {
186  assert(other.container_==container_);
187  return other.position_ - position_;
188  }
189 
191  {
192  assert(other.container_==container_);
193  return other.position_ - position_;
194  }
195 
197  SizeType index () const
198  {
199  return this->position_;
200  }
201 
202  private:
203  C *container_;
204  SizeType position_;
205  };
206 
220  template<typename V>
221  class DenseVector
222  {
223  typedef DenseMatVecTraits<V> Traits;
224  // typedef typename Traits::value_type K;
225 
226  // Curiously recurring template pattern
227  V & asImp() { return static_cast<V&>(*this); }
228  const V & asImp() const { return static_cast<const V&>(*this); }
229 
230  public:
231  //===== type definitions and constants
232 
235 
237  typedef typename Traits::value_type value_type;
238 
240  typedef typename Traits::value_type field_type;
241 
243  typedef typename Traits::value_type block_type;
244 
246  typedef typename Traits::size_type size_type;
247 
249  enum {
252  };
253 
254  //===== assignment from scalar
257  {
258  for (size_type i=0; i<size(); i++)
259  asImp()[i] = k;
260  return asImp();
261  }
262 
263  //===== access to components
264 
267  {
268  return asImp().vec_access(i);
269  }
270 
272  {
273  return asImp().vec_access(i);
274  }
275 
277  size_type size() const
278  {
279  return asImp().vec_size();
280  }
281 
286 
289  {
290  return Iterator(*this,0);
291  }
292 
295  {
296  return Iterator(*this,size());
297  }
298 
302  {
303  return Iterator(*this,size()-1);
304  }
305 
309  {
310  return Iterator(*this,-1);
311  }
312 
315  {
316  return Iterator(*this,std::min(i,size()));
317  }
318 
323 
326  {
327  return ConstIterator(*this,0);
328  }
329 
332  {
333  return ConstIterator(*this,size());
334  }
335 
339  {
340  return ConstIterator(*this,size()-1);
341  }
342 
346  {
347  return ConstIterator(*this,-1);
348  }
349 
352  {
353  return ConstIterator(*this,std::min(i,size()));
354  }
355 
356  //===== vector space arithmetic
357 
359  template <class Other>
361  {
362  assert(y.size() == size());
363  for (size_type i=0; i<size(); i++)
364  (*this)[i] += y[i];
365  return asImp();
366  }
367 
369  template <class Other>
371  {
372  assert(y.size() == size());
373  for (size_type i=0; i<size(); i++)
374  (*this)[i] -= y[i];
375  return asImp();
376  }
377 
379  template <class Other>
381  {
382  derived_type z = asImp();
383  return (z+=b);
384  }
385 
387  template <class Other>
389  {
390  derived_type z = asImp();
391  return (z-=b);
392  }
393 
396  {
397  for (size_type i=0; i<size(); i++)
398  (*this)[i] += k;
399  return asImp();
400  }
401 
404  {
405  for (size_type i=0; i<size(); i++)
406  (*this)[i] -= k;
407  return asImp();
408  }
409 
412  {
413  for (size_type i=0; i<size(); i++)
414  (*this)[i] *= k;
415  return asImp();
416  }
417 
420  {
421  for (size_type i=0; i<size(); i++)
422  (*this)[i] /= k;
423  return asImp();
424  }
425 
427  template <class Other>
428  bool operator== (const DenseVector<Other>& y) const
429  {
430  assert(y.size() == size());
431  for (size_type i=0; i<size(); i++)
432  if ((*this)[i]!=y[i])
433  return false;
434 
435  return true;
436  }
437 
439  template <class Other>
440  bool operator!= (const DenseVector<Other>& y) const
441  {
442  return !operator==(y);
443  }
444 
445 
447  template <class Other>
449  {
450  assert(y.size() == size());
451  for (size_type i=0; i<size(); i++)
452  (*this)[i] += a*y[i];
453  return asImp();
454  }
455 
456  //===== Euclidean scalar product
457 
459  template <class Other>
461  {
462  assert(y.size() == size());
463  value_type result( 0 );
464  for (size_type i=0; i<size(); i++)
465  result += (*this)[i]*y[i];
466  return result;
467  }
468 
469  //===== norms
470 
473  typename FieldTraits<value_type>::real_type result( 0 );
474  for (size_type i=0; i<size(); i++)
475  result += std::abs((*this)[i]);
476  return result;
477  }
478 
479 
482  {
483  typename FieldTraits<value_type>::real_type result( 0 );
484  for (size_type i=0; i<size(); i++)
485  result += fvmeta::absreal((*this)[i]);
486  return result;
487  }
488 
491  {
492  typename FieldTraits<value_type>::real_type result( 0 );
493  for (size_type i=0; i<size(); i++)
494  result += fvmeta::abs2((*this)[i]);
495  return fvmeta::sqrt(result);
496  }
497 
500  {
501  typename FieldTraits<value_type>::real_type result( 0 );
502  for (size_type i=0; i<size(); i++)
503  result += fvmeta::abs2((*this)[i]);
504  return result;
505  }
506 
509  {
510  if (size() == 0)
511  return 0.0;
512 
513  ConstIterator it = begin();
514  typename FieldTraits<value_type>::real_type max = std::abs(*it);
515  for (it = it + 1; it != end(); ++it)
516  max = std::max(max, std::abs(*it));
517 
518  return max;
519  }
520 
523  {
524  if (size() == 0)
525  return 0.0;
526 
527  ConstIterator it = begin();
528  typename FieldTraits<value_type>::real_type max = fvmeta::absreal(*it);
529  for (it = it + 1; it != end(); ++it)
530  max = std::max(max, fvmeta::absreal(*it));
531 
532  return max;
533  }
534 
535  //===== sizes
536 
538  size_type N () const
539  {
540  return size();
541  }
542 
544  size_type dim () const
545  {
546  return size();
547  }
548 
549  };
550 
559  template<typename V>
560  std::ostream& operator<< (std::ostream& s, const DenseVector<V>& v)
561  {
562  for (typename DenseVector<V>::size_type i=0; i<v.size(); i++)
563  s << ((i>0) ? " " : "") << v[i];
564  return s;
565  }
566 
569 } // end namespace
570 
571 #endif // DUNE_DENSEVECTOR_HH