00001
00002
00003 #ifndef DUNE_DIAGONAL_MATRIX_HH
00004 #define DUNE_DIAGONAL_MATRIX_HH
00005
00010 #include <algorithm>
00011 #include <cassert>
00012 #include <cmath>
00013 #include <complex>
00014 #include <cstddef>
00015 #include <initializer_list>
00016 #include <iostream>
00017 #include <memory>
00018
00019 #include <dune/common/boundschecking.hh>
00020 #include <dune/common/densematrix.hh>
00021 #include <dune/common/exceptions.hh>
00022 #include <dune/common/fmatrix.hh>
00023 #include <dune/common/fvector.hh>
00024 #include <dune/common/genericiterator.hh>
00025 #include <dune/common/typetraits.hh>
00026 #include <dune/common/unused.hh>
00027
00028
00029 namespace Dune {
00030
00031 template< class K, int n > class DiagonalRowVectorConst;
00032 template< class K, int n > class DiagonalRowVector;
00033 template< class DiagonalMatrixType > class DiagonalMatrixWrapper;
00034 template< class C, class T, class R> class ContainerWrapperIterator;
00035
00050 template<class K, int n>
00051 class DiagonalMatrix
00052 {
00053 typedef DiagonalMatrixWrapper< DiagonalMatrix<K,n> > WrapperType;
00054
00055 public:
00056
00057
00059 typedef K value_type;
00060 typedef value_type field_type;
00061
00063 typedef K block_type;
00064
00066 typedef std::size_t size_type;
00067
00069 enum {
00071 blocklevel = 1
00072 };
00073
00075 typedef DiagonalRowVector<K,n> row_type;
00076 typedef row_type reference;
00077 typedef row_type row_reference;
00078 typedef DiagonalRowVectorConst<K,n> const_row_type;
00079 typedef const_row_type const_reference;
00080 typedef const_row_type const_row_reference;
00081
00083 enum {
00085 rows = n,
00087 cols = n
00088 };
00089
00090
00091
00092 size_type size () const
00093 {
00094 return rows;
00095 }
00096
00097
00098
00100 DiagonalMatrix () {}
00101
00103 DiagonalMatrix (const K& k)
00104 : diag_(k)
00105 {}
00106
00108 DiagonalMatrix (const FieldVector<K,n>& diag)
00109 : diag_(diag)
00110 {}
00111
00120 DiagonalMatrix (std::initializer_list<K> const &l)
00121 {
00122 std::copy_n(l.begin(), std::min(static_cast<std::size_t>(rows),
00123 l.size()),
00124 diag_.begin());
00125 }
00126
00128 DiagonalMatrix& operator= (const K& k)
00129 {
00130 diag_ = k;
00131 return *this;
00132 }
00133
00135 bool identical(const DiagonalMatrix<K,n>& other) const
00136 {
00137 return (this==&other);
00138 }
00139
00140
00142 typedef ContainerWrapperIterator<const WrapperType, reference, reference> Iterator;
00144 typedef Iterator iterator;
00146 typedef Iterator RowIterator;
00148 typedef typename row_type::Iterator ColIterator;
00149
00151 Iterator begin ()
00152 {
00153 return Iterator(WrapperType(this),0);
00154 }
00155
00157 Iterator end ()
00158 {
00159 return Iterator(WrapperType(this),n);
00160 }
00161
00164 Iterator beforeEnd ()
00165 {
00166 return Iterator(WrapperType(this),n-1);
00167 }
00168
00171 Iterator beforeBegin ()
00172 {
00173 return Iterator(WrapperType(this),-1);
00174 }
00175
00176
00178 typedef ContainerWrapperIterator<const WrapperType, const_reference, const_reference> ConstIterator;
00180 typedef ConstIterator const_iterator;
00182 typedef ConstIterator ConstRowIterator;
00184 typedef typename const_row_type::ConstIterator ConstColIterator;
00185
00187 ConstIterator begin () const
00188 {
00189 return ConstIterator(WrapperType(this),0);
00190 }
00191
00193 ConstIterator end () const
00194 {
00195 return ConstIterator(WrapperType(this),n);
00196 }
00197
00200 ConstIterator beforeEnd() const
00201 {
00202 return ConstIterator(WrapperType(this),n-1);
00203 }
00204
00207 ConstIterator beforeBegin () const
00208 {
00209 return ConstIterator(WrapperType(this),-1);
00210 }
00211
00212
00213
00214
00215
00217 DiagonalMatrix& operator+= (const DiagonalMatrix& y)
00218 {
00219 diag_ += y.diag_;
00220 return *this;
00221 }
00222
00224 DiagonalMatrix& operator-= (const DiagonalMatrix& y)
00225 {
00226 diag_ -= y.diag_;
00227 return *this;
00228 }
00229
00231 DiagonalMatrix& operator+= (const K& k)
00232 {
00233 diag_ += k;
00234 return *this;
00235 }
00236
00238 DiagonalMatrix& operator-= (const K& k)
00239 {
00240 diag_ -= k;
00241 return *this;
00242 }
00243
00245 DiagonalMatrix& operator*= (const K& k)
00246 {
00247 diag_ *= k;
00248 return *this;
00249 }
00250
00252 DiagonalMatrix& operator/= (const K& k)
00253 {
00254 diag_ /= k;
00255 return *this;
00256 }
00257
00258
00259
00261 bool operator==(const DiagonalMatrix& other) const
00262 {
00263 return diag_==other.diagonal();
00264 }
00265
00267 bool operator!=(const DiagonalMatrix& other) const
00268 {
00269 return diag_!=other.diagonal();
00270 }
00271
00272
00273
00274
00276 template<class X, class Y>
00277 void mv (const X& x, Y& y) const
00278 {
00279 #ifdef DUNE_FMatrix_WITH_CHECKING
00280 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00281 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00282 #endif
00283 for (size_type i=0; i<n; ++i)
00284 y[i] = diag_[i] * x[i];
00285 }
00286
00288 template<class X, class Y>
00289 void mtv (const X& x, Y& y) const
00290 {
00291 mv(x, y);
00292 }
00293
00295 template<class X, class Y>
00296 void umv (const X& x, Y& y) const
00297 {
00298 #ifdef DUNE_FMatrix_WITH_CHECKING
00299 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00300 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00301 #endif
00302 for (size_type i=0; i<n; ++i)
00303 y[i] += diag_[i] * x[i];
00304 }
00305
00307 template<class X, class Y>
00308 void umtv (const X& x, Y& y) const
00309 {
00310 #ifdef DUNE_FMatrix_WITH_CHECKING
00311 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00312 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00313 #endif
00314 for (size_type i=0; i<n; ++i)
00315 y[i] += diag_[i] * x[i];
00316 }
00317
00319 template<class X, class Y>
00320 void umhv (const X& x, Y& y) const
00321 {
00322 #ifdef DUNE_FMatrix_WITH_CHECKING
00323 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00324 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00325 #endif
00326 for (size_type i=0; i<n; i++)
00327 y[i] += conjugateComplex(diag_[i])*x[i];
00328 }
00329
00331 template<class X, class Y>
00332 void mmv (const X& x, Y& y) const
00333 {
00334 #ifdef DUNE_FMatrix_WITH_CHECKING
00335 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00336 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00337 #endif
00338 for (size_type i=0; i<n; ++i)
00339 y[i] -= diag_[i] * x[i];
00340 }
00341
00343 template<class X, class Y>
00344 void mmtv (const X& x, Y& y) const
00345 {
00346 #ifdef DUNE_FMatrix_WITH_CHECKING
00347 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00348 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00349 #endif
00350 for (size_type i=0; i<n; ++i)
00351 y[i] -= diag_[i] * x[i];
00352 }
00353
00355 template<class X, class Y>
00356 void mmhv (const X& x, Y& y) const
00357 {
00358 #ifdef DUNE_FMatrix_WITH_CHECKING
00359 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00360 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00361 #endif
00362 for (size_type i=0; i<n; i++)
00363 y[i] -= conjugateComplex(diag_[i])*x[i];
00364 }
00365
00367 template<class X, class Y>
00368 void usmv (const typename FieldTraits<Y>::field_type & alpha,
00369 const X& x, Y& y) const
00370 {
00371 #ifdef DUNE_FMatrix_WITH_CHECKING
00372 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00373 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00374 #endif
00375 for (size_type i=0; i<n; i++)
00376 y[i] += alpha * diag_[i] * x[i];
00377 }
00378
00380 template<class X, class Y>
00381 void usmtv (const typename FieldTraits<Y>::field_type & alpha,
00382 const X& x, Y& y) const
00383 {
00384 #ifdef DUNE_FMatrix_WITH_CHECKING
00385 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00386 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00387 #endif
00388 for (size_type i=0; i<n; i++)
00389 y[i] += alpha * diag_[i] * x[i];
00390 }
00391
00393 template<class X, class Y>
00394 void usmhv (const typename FieldTraits<Y>::field_type & alpha,
00395 const X& x, Y& y) const
00396 {
00397 #ifdef DUNE_FMatrix_WITH_CHECKING
00398 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00399 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00400 #endif
00401 for (size_type i=0; i<n; i++)
00402 y[i] += alpha * conjugateComplex(diag_[i]) * x[i];
00403 }
00404
00405
00406
00408 double frobenius_norm () const
00409 {
00410 return diag_.two_norm();
00411 }
00412
00414 double frobenius_norm2 () const
00415 {
00416 return diag_.two_norm2();
00417 }
00418
00420 double infinity_norm () const
00421 {
00422 return diag_.infinity_norm();
00423 }
00424
00426 double infinity_norm_real () const
00427 {
00428 return diag_.infinity_norm_real();
00429 }
00430
00431
00432
00433
00434
00436 template<class V>
00437 void solve (V& x, const V& b) const
00438 {
00439 for (int i=0; i<n; i++)
00440 x[i] = b[i]/diag_[i];
00441 }
00442
00444 void invert()
00445 {
00446 for (int i=0; i<n; i++)
00447 diag_[i] = 1/diag_[i];
00448 }
00449
00451 K determinant () const
00452 {
00453 K det = diag_[0];
00454 for (int i=1; i<n; i++)
00455 det *= diag_[i];
00456 return det;
00457 }
00458
00459
00460
00461
00462
00464 size_type N () const
00465 {
00466 return n;
00467 }
00468
00470 size_type M () const
00471 {
00472 return n;
00473 }
00474
00475
00476
00477
00478
00480 bool exists (size_type i, size_type j) const
00481 {
00482 DUNE_ASSERT_BOUNDS(i >= 0 && i < n);
00483 DUNE_ASSERT_BOUNDS(j >= 0 && j < n);
00484 return i==j;
00485 }
00486
00487
00488
00490 friend std::ostream& operator<< (std::ostream& s, const DiagonalMatrix<K,n>& a)
00491 {
00492 for (size_type i=0; i<n; i++) {
00493 for (size_type j=0; j<n; j++)
00494 s << ((i==j) ? a.diag_[i] : 0) << " ";
00495 s << std::endl;
00496 }
00497 return s;
00498 }
00499
00501 reference operator[](size_type i)
00502 {
00503 return reference(const_cast<K*>(&diag_[i]), i);
00504 }
00505
00507 const_reference operator[](size_type i) const
00508 {
00509 return const_reference(const_cast<K*>(&diag_[i]), i);
00510 }
00511
00513 const K& diagonal(size_type i) const
00514 {
00515 return diag_[i];
00516 }
00517
00519 K& diagonal(size_type i)
00520 {
00521 return diag_[i];
00522 }
00523
00525 const FieldVector<K,n>& diagonal() const
00526 {
00527 return diag_;
00528 }
00529
00531 FieldVector<K,n>& diagonal()
00532 {
00533 return diag_;
00534 }
00535
00536 private:
00537
00538
00539 FieldVector<K,n> diag_;
00540 };
00541
00542 #ifndef DOXYGEN // hide specialization
00543
00545 template< class K >
00546 class DiagonalMatrix<K, 1> : public FieldMatrix<K, 1, 1>
00547 {
00548 typedef FieldMatrix<K,1,1> Base;
00549 public:
00551 typedef typename Base::size_type size_type;
00552
00554 enum {
00557 blocklevel = 1
00558 };
00559
00560 typedef typename Base::row_type row_type;
00561
00562 typedef typename Base::row_reference row_reference;
00563 typedef typename Base::const_row_reference const_row_reference;
00564
00566 enum {
00569 rows = 1,
00572 cols = 1
00573 };
00574
00575
00577 DiagonalMatrix()
00578 {}
00579
00581 DiagonalMatrix(const K& scalar)
00582 {
00583 (*this)[0][0] = scalar;
00584 }
00585
00587 const K& diagonal(size_type) const
00588 {
00589 return (*this)[0][0];
00590 }
00591
00593 K& diagonal(size_type)
00594 {
00595 return (*this)[0][0];
00596 }
00597
00599 const FieldVector<K,1>& diagonal() const
00600 {
00601 return (*this)[0];
00602 }
00603
00605 FieldVector<K,1>& diagonal()
00606 {
00607 return (*this)[0];
00608 }
00609
00610 };
00611 #endif
00612
00613
00614 template<class DiagonalMatrixType>
00615 class DiagonalMatrixWrapper
00616 {
00617 typedef typename DiagonalMatrixType::reference reference;
00618 typedef typename DiagonalMatrixType::const_reference const_reference;
00619 typedef typename DiagonalMatrixType::field_type K;
00620 typedef DiagonalRowVector<K, DiagonalMatrixType::rows> row_type;
00621 typedef std::size_t size_type;
00622 typedef DiagonalMatrixWrapper< DiagonalMatrixType> MyType;
00623
00624 friend class ContainerWrapperIterator<const MyType, reference, reference>;
00625 friend class ContainerWrapperIterator<const MyType, const_reference, const_reference>;
00626
00627 public:
00628
00629 DiagonalMatrixWrapper() :
00630 mat_(0)
00631 {}
00632
00633 DiagonalMatrixWrapper(const DiagonalMatrixType* mat) :
00634 mat_(const_cast<DiagonalMatrixType*>(mat))
00635 {}
00636
00637 size_type realIndex(int i) const
00638 {
00639 return i;
00640 }
00641
00642 row_type* pointer(int i) const
00643 {
00644 row_ = row_type(&(mat_->diagonal(i)), i);
00645 return &row_;
00646 }
00647
00648 bool identical(const DiagonalMatrixWrapper& other) const
00649 {
00650 return mat_==other.mat_;
00651 }
00652
00653 private:
00654
00655 mutable DiagonalMatrixType* mat_;
00656 mutable row_type row_;
00657 };
00658
00662 template< class K, int n >
00663 class DiagonalRowVectorConst
00664 {
00665 template<class DiagonalMatrixType>
00666 friend class DiagonalMatrixWrapper;
00667 friend class ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&>;
00668
00669 public:
00670
00671 enum { dimension = n };
00672
00673
00674
00675
00676
00678 typedef K field_type;
00679
00681 typedef K block_type;
00682
00684 typedef std::size_t size_type;
00685
00687 enum {
00689 blocklevel = 1
00690 };
00691
00693 enum {
00695 size = n
00696 };
00697
00699 DiagonalRowVectorConst() :
00700 p_(0),
00701 row_(0)
00702 {}
00703
00705 explicit DiagonalRowVectorConst (K* p, int col) :
00706 p_(p),
00707 row_(col)
00708 {}
00709
00710
00711
00713 const K& operator[] (size_type i) const
00714 {
00715 DUNE_UNUSED_PARAMETER(i);
00716 DUNE_ASSERT_BOUNDS(i == row_);
00717 return *p_;
00718 }
00719
00720
00721
00722 bool identical(const DiagonalRowVectorConst<K,n>& other) const
00723 {
00724 return ((p_ == other.p_)and (row_ == other.row_));
00725 }
00726
00728 typedef ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&> ConstIterator;
00730 typedef ConstIterator const_iterator;
00731
00733 ConstIterator begin () const
00734 {
00735 return ConstIterator(*this,0);
00736 }
00737
00739 ConstIterator end () const
00740 {
00741 return ConstIterator(*this,1);
00742 }
00743
00746 ConstIterator beforeEnd() const
00747 {
00748 return ConstIterator(*this,0);
00749 }
00750
00753 ConstIterator beforeBegin () const
00754 {
00755 return ConstIterator(*this,-1);
00756 }
00757
00759 bool operator== (const DiagonalRowVectorConst& y) const
00760 {
00761 return ((p_==y.p_)and (row_==y.row_));
00762 }
00763
00764
00765
00767 size_type N () const
00768 {
00769 return n;
00770 }
00771
00773 size_type dim () const
00774 {
00775 return n;
00776 }
00777
00779 size_type rowIndex() const
00780 {
00781 return row_;
00782 }
00783
00785 const K& diagonal() const
00786 {
00787 return *p_;
00788 }
00789
00790 protected:
00791
00792 size_type realIndex(int i) const
00793 {
00794 return rowIndex();
00795 }
00796
00797 K* pointer(size_type i) const
00798 {
00799 return const_cast<K*>(p_);
00800 }
00801
00802 DiagonalRowVectorConst* operator&()
00803 {
00804 return this;
00805 }
00806
00807
00808 K* p_;
00809 size_type row_;
00810 };
00811
00812 template< class K, int n >
00813 class DiagonalRowVector : public DiagonalRowVectorConst<K,n>
00814 {
00815 template<class DiagonalMatrixType>
00816 friend class DiagonalMatrixWrapper;
00817 friend class ContainerWrapperIterator<DiagonalRowVector<K,n>, K, K&>;
00818
00819 public:
00820
00821
00822
00823
00825 typedef K field_type;
00826
00828 typedef K block_type;
00829
00831 typedef std::size_t size_type;
00832
00834 DiagonalRowVector() : DiagonalRowVectorConst<K,n>()
00835 {}
00836
00838 explicit DiagonalRowVector (K* p, int col) : DiagonalRowVectorConst<K,n>(p, col)
00839 {}
00840
00841
00843 DiagonalRowVector& operator= (const K& k)
00844 {
00845 *p_ = k;
00846 return *this;
00847 }
00848
00849
00850
00852 K& operator[] (size_type i)
00853 {
00854 DUNE_UNUSED_PARAMETER(i);
00855 DUNE_ASSERT_BOUNDS(i == row_);
00856 return *p_;
00857 }
00858
00860 typedef ContainerWrapperIterator<DiagonalRowVector<K,n>, K, K&> Iterator;
00862 typedef Iterator iterator;
00863
00865 Iterator begin ()
00866 {
00867 return Iterator(*this, 0);
00868 }
00869
00871 Iterator end ()
00872 {
00873 return Iterator(*this, 1);
00874 }
00875
00878 Iterator beforeEnd ()
00879 {
00880 return Iterator(*this, 0);
00881 }
00882
00885 Iterator beforeBegin ()
00886 {
00887 return Iterator(*this, -1);
00888 }
00889
00891 typedef ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&> ConstIterator;
00893 typedef ConstIterator const_iterator;
00894
00895 using DiagonalRowVectorConst<K,n>::identical;
00896 using DiagonalRowVectorConst<K,n>::operator[];
00897 using DiagonalRowVectorConst<K,n>::operator==;
00898 using DiagonalRowVectorConst<K,n>::begin;
00899 using DiagonalRowVectorConst<K,n>::end;
00900 using DiagonalRowVectorConst<K,n>::beforeEnd;
00901 using DiagonalRowVectorConst<K,n>::beforeBegin;
00902 using DiagonalRowVectorConst<K,n>::N;
00903 using DiagonalRowVectorConst<K,n>::dim;
00904 using DiagonalRowVectorConst<K,n>::rowIndex;
00905 using DiagonalRowVectorConst<K,n>::diagonal;
00906
00907 protected:
00908
00909 DiagonalRowVector* operator&()
00910 {
00911 return this;
00912 }
00913
00914 private:
00915
00916 using DiagonalRowVectorConst<K,n>::p_;
00917 using DiagonalRowVectorConst<K,n>::row_;
00918 };
00919
00920
00921
00922 template<class K, int n>
00923 struct const_reference< DiagonalRowVector<K,n> >
00924 {
00925 typedef DiagonalRowVectorConst<K,n> type;
00926 };
00927
00928 template<class K, int n>
00929 struct const_reference< DiagonalRowVectorConst<K,n> >
00930 {
00931 typedef DiagonalRowVectorConst<K,n> type;
00932 };
00933
00934 template<class K, int n>
00935 struct mutable_reference< DiagonalRowVector<K,n> >
00936 {
00937 typedef DiagonalRowVector<K,n> type;
00938 };
00939
00940 template<class K, int n>
00941 struct mutable_reference< DiagonalRowVectorConst<K,n> >
00942 {
00943 typedef DiagonalRowVector<K,n> type;
00944 };
00945
00946
00947
00970 template<class CW, class T, class R>
00971 class ContainerWrapperIterator : public BidirectionalIteratorFacade<ContainerWrapperIterator<CW,T,R>,T, R, int>
00972 {
00973 typedef typename std::remove_const<CW>::type NonConstCW;
00974
00975 friend class ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type>;
00976 friend class ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type>;
00977
00978 typedef ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type> MyType;
00979 typedef ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type> MyConstType;
00980
00981 public:
00982
00983
00984 ContainerWrapperIterator() :
00985 containerWrapper_(),
00986 position_(0)
00987 {}
00988
00989 ContainerWrapperIterator(CW containerWrapper, int position) :
00990 containerWrapper_(containerWrapper),
00991 position_(position)
00992 {}
00993
00994 template<class OtherContainerWrapperIteratorType>
00995 ContainerWrapperIterator(OtherContainerWrapperIteratorType& other) :
00996 containerWrapper_(other.containerWrapper_),
00997 position_(other.position_)
00998 {}
00999
01000 ContainerWrapperIterator(const MyType& other) :
01001 containerWrapper_(other.containerWrapper_),
01002 position_(other.position_)
01003 {}
01004
01005 ContainerWrapperIterator(const MyConstType& other) :
01006 containerWrapper_(other.containerWrapper_),
01007 position_(other.position_)
01008 {}
01009
01010 template<class OtherContainerWrapperIteratorType>
01011 ContainerWrapperIterator& operator=(OtherContainerWrapperIteratorType& other)
01012 {
01013 containerWrapper_ = other.containerWrapper_;
01014 position_ = other.position_;
01015 return *this;
01016 }
01017
01018
01019
01020 T* operator->() const
01021 {
01022 return containerWrapper_.pointer(position_);
01023 }
01024
01025
01026 bool equals(const MyType& other) const
01027 {
01028 return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
01029 }
01030
01031 bool equals(const MyConstType& other) const
01032 {
01033 return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
01034 }
01035
01036 R dereference() const
01037 {
01038 return *containerWrapper_.pointer(position_);
01039 }
01040
01041 void increment()
01042 {
01043 ++position_;
01044 }
01045
01046
01047 void decrement()
01048 {
01049 --position_;
01050 }
01051
01052
01053 R elementAt(int i) const
01054 {
01055 return *containerWrapper_.pointer(position_+i);
01056 }
01057
01058 void advance(int n)
01059 {
01060 position_=position_+n;
01061 }
01062
01063 template<class OtherContainerWrapperIteratorType>
01064 std::ptrdiff_t distanceTo(OtherContainerWrapperIteratorType& other) const
01065 {
01066 assert(containerWrapper_.identical(other));
01067 return other.position_ - position_;
01068 }
01069
01070 std::ptrdiff_t index() const
01071 {
01072 return containerWrapper_.realIndex(position_);
01073 }
01074
01075 private:
01076 NonConstCW containerWrapper_;
01077 size_t position_;
01078 };
01079
01080 template <class DenseMatrix, class field, int N>
01081 struct DenseMatrixAssigner<DenseMatrix, DiagonalMatrix<field, N>> {
01082 static void apply(DenseMatrix& denseMatrix,
01083 DiagonalMatrix<field, N> const& rhs) {
01084 DUNE_ASSERT_BOUNDS(denseMatrix.M() == N);
01085 DUNE_ASSERT_BOUNDS(denseMatrix.N() == N);
01086 denseMatrix = field(0);
01087 for (int i = 0; i < N; ++i)
01088 denseMatrix[i][i] = rhs.diagonal()[i];
01089 }
01090 };
01091
01092 }
01093 #endif