dune-istl  2.2.0
diagonalmatrix.hh
Go to the documentation of this file.
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
00002 // vi: set et ts=4 sw=4 sts=4:
00003 #ifndef DUNE_DIAGONAL_MATRIX_HH
00004 #define DUNE_DIAGONAL_MATRIX_HH
00005 
00010 #include<cmath>
00011 #include<cstddef>
00012 #include<complex>
00013 #include<iostream>
00014 #include<memory>
00015 #include <dune/common/exceptions.hh>
00016 #include <dune/common/fmatrix.hh>
00017 #include <dune/common/fvector.hh>
00018 #include <dune/common/typetraits.hh>
00019 #include <dune/common/genericiterator.hh>
00020 
00021 
00022 
00023 namespace Dune {
00024 
00025 template< class K, int n > class DiagonalRowVectorConst;
00026 template< class K, int n > class DiagonalRowVector;
00027 template< class DiagonalMatrixType > class DiagonalMatrixWrapper;
00028 template< class C, class T, class R> class ContainerWrapperIterator;
00029 
00044 template<class K, int n>
00045 class DiagonalMatrix
00046 {
00047     typedef DiagonalMatrixWrapper< DiagonalMatrix<K,n> > WrapperType;
00048 
00049     public:
00050     //===== type definitions and constants
00051 
00053     typedef K field_type;
00054 
00056     typedef K block_type;
00057 
00059     typedef std::size_t size_type;
00060 
00062     enum {
00064         blocklevel = 1
00065     };
00066 
00068     typedef DiagonalRowVector<K,n> row_type;
00069     typedef row_type reference;
00070     typedef DiagonalRowVectorConst<K,n> const_row_type;
00071     typedef const_row_type const_reference;
00072 
00074     enum {
00076         rows = n,
00078         cols = n
00079     };
00080 
00081 
00082 
00083     //===== constructors
00084 
00086     DiagonalMatrix () {}
00087 
00089     DiagonalMatrix (const K& k)
00090         : diag_(k)
00091     {}
00092 
00094     DiagonalMatrix (const FieldVector<K,n>& diag)
00095         : diag_(diag)
00096     {}
00097 
00098 
00100     DiagonalMatrix& operator= (const K& k)
00101     {
00102         diag_ = k;
00103         return *this;
00104     }
00105 
00107     bool identical(const DiagonalMatrix<K,n>& other) const
00108     {
00109         return (this==&other);
00110     }
00111 
00112     //===== iterator interface to rows of the matrix
00114     typedef ContainerWrapperIterator<const WrapperType, reference, reference> Iterator;
00116     typedef Iterator iterator;
00118     typedef Iterator RowIterator;
00120     typedef typename row_type::Iterator ColIterator;
00121 
00123     Iterator begin ()
00124     {
00125         return Iterator(WrapperType(this),0);
00126     }
00127 
00129     Iterator end ()
00130     {
00131         return Iterator(WrapperType(this),n);
00132     }
00133 
00136         Iterator beforeEnd ()
00137     {
00138         return Iterator(WrapperType(this),n-1);
00139     }
00140 
00143         Iterator beforeBegin ()
00144     {
00145         return Iterator(WrapperType(this),-1);
00146     }
00147 
00148 
00150     typedef ContainerWrapperIterator<const WrapperType, const_reference, const_reference> ConstIterator;
00152     typedef ConstIterator const_iterator;
00154     typedef ConstIterator ConstRowIterator;
00156     typedef typename const_row_type::ConstIterator ConstColIterator;
00157 
00159     ConstIterator begin () const
00160     {
00161         return ConstIterator(WrapperType(this),0);
00162     }
00163 
00165     ConstIterator end () const
00166     {
00167         return ConstIterator(WrapperType(this),n);
00168     }
00169 
00172         ConstIterator beforeEnd() const
00173     {
00174         return ConstIterator(WrapperType(this),n-1);
00175     }
00176 
00179         ConstIterator beforeBegin () const
00180     {
00181         return ConstIterator(WrapperType(this),-1);
00182     }
00183 
00184 
00185 
00186     //===== vector space arithmetic
00187 
00189     DiagonalMatrix& operator+= (const DiagonalMatrix& y)
00190     {
00191         diag_ += y.diag_;
00192         return *this;
00193     }
00194 
00196     DiagonalMatrix& operator-= (const DiagonalMatrix& y)
00197     {
00198         diag_ -= y.diag_;
00199         return *this;
00200     }
00201 
00203     DiagonalMatrix& operator+= (const K& k)
00204     {
00205         diag_ += k;
00206         return *this;
00207     }
00208 
00210     DiagonalMatrix& operator-= (const K& k)
00211     {
00212         diag_ -= k;
00213         return *this;
00214     }
00215 
00217     DiagonalMatrix& operator*= (const K& k)
00218     {
00219         diag_ *= k;
00220         return *this;
00221     }
00222 
00224     DiagonalMatrix& operator/= (const K& k)
00225     {
00226         diag_ /= k;
00227         return *this;
00228     }
00229 
00230     //===== comparison ops
00231 
00233     bool operator==(const DiagonalMatrix& other) const
00234     {
00235         return diag_==other.diagonal();
00236     }
00237 
00239     bool operator!=(const DiagonalMatrix& other) const
00240     {
00241         return diag_!=other.diagonal();
00242     }
00243 
00244 
00245     //===== linear maps
00246 
00248     template<class X, class Y>
00249     void mv (const X& x, Y& y) const
00250     {
00251 #ifdef DUNE_FMatrix_WITH_CHECKING
00252         if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00253         if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00254 #endif
00255         for (size_type i=0; i<n; ++i)
00256             y[i] = diag_[i] * x[i];
00257     }
00258 
00260     template<class X, class Y>
00261     void mtv (const X& x, Y& y) const
00262     {
00263         mv(x, y);
00264     }
00265 
00267     template<class X, class Y>
00268     void umv (const X& x, Y& y) const
00269     {
00270 #ifdef DUNE_FMatrix_WITH_CHECKING
00271         if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00272         if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00273 #endif
00274         for (size_type i=0; i<n; ++i)
00275             y[i] += diag_[i] * x[i];
00276     }
00277 
00279     template<class X, class Y>
00280     void umtv (const X& x, Y& y) const
00281     {
00282 #ifdef DUNE_FMatrix_WITH_CHECKING
00283         if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00284         if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00285 #endif
00286         for (size_type i=0; i<n; ++i)
00287             y[i] += diag_[i] * x[i];
00288     }
00289 
00291     template<class X, class Y>
00292     void umhv (const X& x, Y& y) const
00293     {
00294 #ifdef DUNE_FMatrix_WITH_CHECKING
00295         if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00296         if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00297 #endif
00298         for (size_type i=0; i<n; i++)
00299             y[i] += conjugateComplex(diag_[i])*x[i];
00300     }
00301 
00303     template<class X, class Y>
00304     void mmv (const X& x, Y& y) const
00305     {
00306 #ifdef DUNE_FMatrix_WITH_CHECKING
00307         if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00308         if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00309 #endif
00310         for (size_type i=0; i<n; ++i)
00311             y[i] -= diag_[i] * x[i];
00312     }
00313 
00315     template<class X, class Y>
00316     void mmtv (const X& x, Y& y) const
00317     {
00318 #ifdef DUNE_FMatrix_WITH_CHECKING
00319         if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00320         if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00321 #endif
00322         for (size_type i=0; i<n; ++i)
00323             y[i] -= diag_[i] * x[i];
00324     }
00325 
00327     template<class X, class Y>
00328     void mmhv (const X& x, Y& y) const
00329     {
00330 #ifdef DUNE_FMatrix_WITH_CHECKING
00331         if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00332         if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00333 #endif
00334         for (size_type i=0; i<n; i++)
00335             y[i] -= conjugateComplex(diag_[i])*x[i];
00336     }
00337 
00339     template<class X, class Y>
00340     void usmv (const K& alpha, const X& x, Y& y) const
00341     {
00342 #ifdef DUNE_FMatrix_WITH_CHECKING
00343         if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00344         if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00345 #endif
00346         for (size_type i=0; i<n; i++)
00347             y[i] += alpha * diag_[i] * x[i];
00348     }
00349 
00351     template<class X, class Y>
00352     void usmtv (const K& alpha, const X& x, Y& y) const
00353     {
00354 #ifdef DUNE_FMatrix_WITH_CHECKING
00355         if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00356         if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00357 #endif
00358         for (size_type i=0; i<n; i++)
00359             y[i] += alpha * diag_[i] * x[i];
00360     }
00361 
00363     template<class X, class Y>
00364     void usmhv (const K& alpha, const X& x, Y& y) const
00365     {
00366 #ifdef DUNE_FMatrix_WITH_CHECKING
00367         if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
00368         if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
00369 #endif
00370         for (size_type i=0; i<n; i++)
00371             y[i] += alpha * conjugateComplex(diag_[i]) * x[i];
00372     }
00373 
00374     //===== norms
00375 
00377     double frobenius_norm () const
00378     {
00379         return diag_.two_norm();
00380     }
00381 
00383     double frobenius_norm2 () const
00384     {
00385         return diag_.two_norm2();
00386     }
00387 
00389     double infinity_norm () const
00390     {
00391         return diag_.infinity_norm();
00392     }
00393 
00395     double infinity_norm_real () const
00396     {
00397         return diag_.infinity_norm_real();
00398     }
00399 
00400 
00401 
00402     //===== solve
00403 
00405     template<class V>
00406     void solve (V& x, const V& b) const
00407     {
00408         for (int i=0; i<n; i++)
00409             x[i] = b[i]/diag_[i];
00410     }
00411 
00413     void invert()
00414     {
00415         for (int i=0; i<n; i++)
00416             diag_[i] = 1/diag_[i];
00417     }
00418 
00420     K determinant () const
00421     {
00422         K det = diag_[0];
00423         for (int i=1; i<n; i++)
00424             det *= diag_[i];
00425         return det;
00426     }
00427 
00428 
00429 
00430     //===== sizes
00431 
00433     size_type N () const
00434     {
00435         return n;
00436     }
00437 
00439     size_type M () const
00440     {
00441         return n;
00442     }
00443 
00444 
00445 
00446     //===== query
00447 
00449     bool exists (size_type i, size_type j) const
00450     {
00451 #ifdef DUNE_FMatrix_WITH_CHECKING
00452         if (i<0 || i>=n) DUNE_THROW(FMatrixError,"row index out of range");
00453         if (j<0 || j>=n) DUNE_THROW(FMatrixError,"column index out of range");
00454 #endif
00455         return i==j;
00456     }
00457 
00458 
00459 
00461     friend std::ostream& operator<< (std::ostream& s, const DiagonalMatrix<K,n>& a)
00462     {
00463         for (size_type i=0; i<n; i++) {
00464             for (size_type j=0; j<n; j++)
00465                 s << ((i==j) ? a.diag_[i] : 0) << " ";
00466             s << std::endl;
00467         }
00468         return s;
00469     }
00470 
00472     reference operator[](size_type i)
00473     {
00474         return reference(const_cast<K*>(&diag_[i]), i);
00475     }
00476 
00478     const_reference operator[](size_type i) const
00479     {
00480         return const_reference(const_cast<K*>(&diag_[i]), i);
00481     }
00482 
00484     const K& diagonal(size_type i) const
00485     {
00486         return diag_[i];
00487     }
00488 
00490     K& diagonal(size_type i)
00491     {
00492         return diag_[i];
00493     }
00494 
00496     const FieldVector<K,n>& diagonal() const
00497     {
00498         return diag_;
00499     }
00500 
00502     FieldVector<K,n>& diagonal()
00503     {
00504         return diag_;
00505     }
00506 
00507     private:
00508 
00509     // the data, a FieldVector storing the diagonal
00510     FieldVector<K,n> diag_;
00511 };
00512 
00513 #ifndef DOXYGEN // hide specialization
00514 
00516 template< class K >
00517 class DiagonalMatrix<K, 1> : public FieldMatrix<K, 1, 1>
00518 {
00519         typedef FieldMatrix<K,1,1> Base;
00520     public:
00522         typedef typename Base::size_type size_type;
00523 
00525         enum {
00528             blocklevel = 1
00529         };
00530 
00531         typedef typename Base::row_type row_type;
00532 
00533         typedef typename Base::row_reference row_reference;
00534         typedef typename Base::const_row_reference const_row_reference;
00535 
00537         enum {
00540             rows = 1,
00543             cols = 1
00544         };
00545 
00546 
00548         DiagonalMatrix()
00549         {}
00550 
00552         DiagonalMatrix(const K& scalar)
00553         {
00554             (*this)[0][0] = scalar;
00555         }
00556 
00558         template<typename T>
00559         DiagonalMatrix(const T& t)
00560         {
00561             DenseMatrixAssigner<Conversion<T,K>::exists>::assign(*this, t);
00562         }
00563 
00565         const K& diagonal(size_type) const
00566         {
00567             return (*this)[0][0];
00568         }
00569 
00571         K& diagonal(size_type)
00572         {
00573             return (*this)[0][0];
00574         }
00575 
00577         const FieldVector<K,1>& diagonal() const
00578         {
00579             return (*this)[0];
00580         }
00581 
00583         FieldVector<K,1>& diagonal()
00584         {
00585             return (*this)[0];
00586         }
00587 
00588 };
00589 #endif
00590 
00591 
00592 template<class DiagonalMatrixType>
00593 class DiagonalMatrixWrapper
00594 {
00595         typedef typename DiagonalMatrixType::reference reference;
00596         typedef typename DiagonalMatrixType::const_reference const_reference;
00597         typedef typename DiagonalMatrixType::field_type K;
00598         typedef DiagonalRowVector<K, DiagonalMatrixType::rows> row_type;
00599         typedef std::size_t size_type;
00600         typedef DiagonalMatrixWrapper< DiagonalMatrixType> MyType;
00601 
00602         friend class ContainerWrapperIterator<const MyType, reference, reference>;
00603         friend class ContainerWrapperIterator<const MyType, const_reference, const_reference>;
00604 
00605     public:
00606 
00607         DiagonalMatrixWrapper() :
00608             mat_(0)
00609         {}
00610 
00611         DiagonalMatrixWrapper(const DiagonalMatrixType* mat) :
00612             mat_(const_cast<DiagonalMatrixType*>(mat))
00613         {}
00614 
00615         size_type realIndex(int i) const
00616         {
00617             return i;
00618         }
00619 
00620         row_type* pointer(int i) const
00621         {
00622             row_ = row_type(&(mat_->diagonal(i)), i);
00623             return &row_;
00624         }
00625 
00626         bool identical(const DiagonalMatrixWrapper& other) const
00627         {
00628             return mat_==other.mat_;
00629         }
00630 
00631     private:
00632 
00633         mutable DiagonalMatrixType* mat_;
00634         mutable row_type row_;
00635 };
00636 
00640 template< class K, int n >
00641 class DiagonalRowVectorConst
00642 {
00643     template<class DiagonalMatrixType>
00644     friend class DiagonalMatrixWrapper;
00645     friend class ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&>;
00646 
00647 public:
00648     // remember size of vector
00649     enum { dimension = n };
00650 
00651     // standard constructor and everything is sufficient ...
00652 
00653     //===== type definitions and constants
00654 
00656     typedef K field_type;
00657 
00659     typedef K block_type;
00660 
00662     typedef std::size_t size_type;
00663 
00665     enum {
00667         blocklevel = 1
00668     };
00669 
00671     enum {
00673         size = n
00674     };
00675 
00677     DiagonalRowVectorConst() :
00678         p_(0),
00679         row_(0)
00680     {}
00681 
00683     explicit DiagonalRowVectorConst (K* p, int col) :
00684         p_(p),
00685         row_(col)
00686     {}
00687 
00688     //===== access to components
00689 
00691     const K& operator[] (size_type i) const
00692     {
00693 #ifdef DUNE_FMatrix_WITH_CHECKING
00694         if (i!=row_)
00695             DUNE_THROW(FMatrixError,"index is contained in pattern");
00696 #endif
00697         return *p_;
00698     }
00699 
00700     // check if row is identical to other row (not only identical values)
00701     // since this is a proxy class we need to check equality of the stored pointer
00702     bool identical(const DiagonalRowVectorConst<K,n>& other) const
00703     {
00704         return ((p_ == other.p_) and (row_ == other.row_));
00705     }
00706 
00708     typedef ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&> ConstIterator;
00710     typedef ConstIterator const_iterator;
00711 
00713     ConstIterator begin () const
00714     {
00715         return ConstIterator(*this,0);
00716     }
00717 
00719     ConstIterator end () const
00720     {
00721         return ConstIterator(*this,1);
00722     }
00723 
00726         ConstIterator beforeEnd() const
00727     {
00728         return ConstIterator(*this,0);
00729     }
00730 
00733         ConstIterator beforeBegin () const
00734     {
00735         return ConstIterator(*this,-1);
00736     }
00737 
00739     bool operator== (const DiagonalRowVectorConst& y) const
00740     {
00741         return ((p_==y.p_) and (row_==y.row_));
00742     }
00743 
00744     //===== sizes
00745 
00747     size_type N () const
00748     {
00749         return n;
00750     }
00751 
00753     size_type dim () const
00754     {
00755         return n;
00756     }
00757 
00759     size_type rowIndex() const
00760     {
00761         return row_;
00762     }
00763 
00765     const K& diagonal() const
00766     {
00767         return *p_;
00768     }
00769 
00770 protected:
00771 
00772     size_type realIndex(int i) const
00773     {
00774         return rowIndex();
00775     }
00776 
00777     K* pointer(size_type i) const
00778     {
00779         return const_cast<K*>(p_);
00780     }
00781 
00782     DiagonalRowVectorConst* operator&()
00783     {
00784         return this;
00785     }
00786 
00787     // the data, very simply a pointer to the diagonal value and the row number
00788     K* p_;
00789     size_type row_;
00790 };
00791 
00792 template< class K, int n >
00793 class DiagonalRowVector : public DiagonalRowVectorConst<K,n>
00794 {
00795     template<class DiagonalMatrixType>
00796     friend class DiagonalMatrixWrapper;
00797     friend class ContainerWrapperIterator<DiagonalRowVector<K,n>, K, K&>;
00798 
00799 public:
00800     // standard constructor and everything is sufficient ...
00801 
00802     //===== type definitions and constants
00803 
00805     typedef K field_type;
00806 
00808     typedef K block_type;
00809 
00811     typedef std::size_t size_type;
00812 
00814     DiagonalRowVector() : DiagonalRowVectorConst<K,n>()
00815     {}
00816 
00818     explicit DiagonalRowVector (K* p, int col) : DiagonalRowVectorConst<K,n>(p, col)
00819     {}
00820 
00821     //===== assignment from scalar
00823     DiagonalRowVector& operator= (const K& k)
00824     {
00825         *p_ = k;
00826         return *this;
00827     }
00828 
00829     //===== access to components
00830 
00832     K& operator[] (size_type i)
00833     {
00834 #ifdef DUNE_FMatrix_WITH_CHECKING
00835         if (i!=row_)
00836             DUNE_THROW(FMatrixError,"index is contained in pattern");
00837 #endif
00838         return *p_;
00839     }
00840 
00842     typedef ContainerWrapperIterator<DiagonalRowVector<K,n>, K, K&> Iterator;
00844     typedef Iterator iterator;
00845 
00847     Iterator begin ()
00848     {
00849         return Iterator(*this, 0);
00850     }
00851 
00853     Iterator end ()
00854     {
00855         return Iterator(*this, 1);
00856     }
00857 
00860         Iterator beforeEnd ()
00861     {
00862         return Iterator(*this, 0);
00863     }
00864 
00867         Iterator beforeBegin ()
00868     {
00869         return Iterator(*this, -1);
00870     }
00871 
00873     typedef ContainerWrapperIterator<DiagonalRowVectorConst<K,n>, const K, const K&> ConstIterator;
00875     typedef ConstIterator const_iterator;
00876 
00877     using DiagonalRowVectorConst<K,n>::identical;
00878     using DiagonalRowVectorConst<K,n>::operator[];
00879     using DiagonalRowVectorConst<K,n>::operator==;
00880     using DiagonalRowVectorConst<K,n>::begin;
00881     using DiagonalRowVectorConst<K,n>::end;
00882     using DiagonalRowVectorConst<K,n>::beforeEnd;
00883     using DiagonalRowVectorConst<K,n>::beforeBegin;
00884     using DiagonalRowVectorConst<K,n>::N;
00885     using DiagonalRowVectorConst<K,n>::dim;
00886     using DiagonalRowVectorConst<K,n>::rowIndex;
00887     using DiagonalRowVectorConst<K,n>::diagonal;
00888 
00889 protected:
00890 
00891     DiagonalRowVector* operator&()
00892     {
00893         return this;
00894     }
00895 
00896 private:
00897 
00898     using DiagonalRowVectorConst<K,n>::p_;
00899     using DiagonalRowVectorConst<K,n>::row_;
00900 };
00901 
00902 
00903 // implement type traits
00904 template<class K, int n>
00905 struct const_reference< DiagonalRowVector<K,n> >
00906 {
00907     typedef DiagonalRowVectorConst<K,n> type;
00908 };
00909 
00910 template<class K, int n>
00911 struct const_reference< DiagonalRowVectorConst<K,n> >
00912 {
00913     typedef DiagonalRowVectorConst<K,n> type;
00914 };
00915 
00916 template<class K, int n>
00917 struct mutable_reference< DiagonalRowVector<K,n> >
00918 {
00919     typedef DiagonalRowVector<K,n> type;
00920 };
00921 
00922 template<class K, int n>
00923 struct mutable_reference< DiagonalRowVectorConst<K,n> >
00924 {
00925     typedef DiagonalRowVector<K,n> type;
00926 };
00927 
00928 
00929 
00952 template<class CW, class T, class R>
00953 class ContainerWrapperIterator : public BidirectionalIteratorFacade<ContainerWrapperIterator<CW,T,R>,T, R, int>
00954 {
00955         typedef typename remove_const<CW>::type NonConstCW;
00956 
00957         friend class ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type>;
00958         friend class ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type>;
00959 
00960         typedef ContainerWrapperIterator<CW, typename mutable_reference<T>::type, typename mutable_reference<R>::type> MyType;
00961         typedef ContainerWrapperIterator<CW, typename const_reference<T>::type, typename const_reference<R>::type> MyConstType;
00962 
00963     public:
00964 
00965         // Constructors needed by the facade iterators.
00966         ContainerWrapperIterator():
00967             containerWrapper_(),
00968             position_(0)
00969         {}
00970 
00971         ContainerWrapperIterator(CW containerWrapper, int position) :
00972             containerWrapper_(containerWrapper),
00973             position_(position)
00974         {}
00975 
00976         template<class OtherContainerWrapperIteratorType>
00977         ContainerWrapperIterator(OtherContainerWrapperIteratorType& other):
00978             containerWrapper_(other.containerWrapper_),
00979             position_(other.position_)
00980         {}
00981 
00982         ContainerWrapperIterator(const MyType& other):
00983             containerWrapper_(other.containerWrapper_),
00984             position_(other.position_)
00985         {}
00986 
00987         ContainerWrapperIterator(const MyConstType& other):
00988             containerWrapper_(other.containerWrapper_),
00989             position_(other.position_)
00990         {}
00991 
00992         template<class OtherContainerWrapperIteratorType>
00993         ContainerWrapperIterator& operator=(OtherContainerWrapperIteratorType& other)
00994         {
00995             containerWrapper_ = other.containerWrapper_;
00996             position_ = other.position_;
00997         }
00998 
00999         // This operator is needed since we can not get the address of the
01000         // temporary object returned by dereference
01001         T* operator->() const
01002         {
01003             return containerWrapper_.pointer(position_);
01004         }
01005 
01006         // Methods needed by the forward iterator
01007         bool equals(const MyType& other) const
01008         {
01009             return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
01010         }
01011 
01012         bool equals(const MyConstType& other) const
01013         {
01014             return position_ == other.position_ && containerWrapper_.identical(other.containerWrapper_);
01015         }
01016 
01017         R dereference() const
01018         {
01019             return *containerWrapper_.pointer(position_);
01020         }
01021 
01022         void increment()
01023         {
01024             ++position_;
01025         }
01026 
01027         // Additional function needed by BidirectionalIterator
01028         void decrement()
01029         {
01030             --position_;
01031         }
01032 
01033         // Additional function needed by RandomAccessIterator
01034         R elementAt(int i) const
01035         {
01036             return *containerWrapper_.pointer(position_+i);
01037         }
01038 
01039         void advance(int n)
01040         {
01041             position_=position_+n;
01042         }
01043 
01044         template<class OtherContainerWrapperIteratorType>
01045         std::ptrdiff_t distanceTo(OtherContainerWrapperIteratorType& other) const
01046         {
01047             assert(containerWrapper_.identical(other));
01048             return other.position_ - position_;
01049         }
01050 
01051         std::ptrdiff_t index() const
01052         {
01053             return containerWrapper_.realIndex(position_);
01054         }
01055 
01056     private:
01057         NonConstCW containerWrapper_;
01058         size_t position_;
01059 };
01060 
01061 
01062 
01063 template<class K, int n>
01064 void istl_assign_to_fmatrix(FieldMatrix<K,n,n>& fm, const DiagonalMatrix<K,n>& s)
01065 {
01066     fm = K();
01067     for(int i=0; i<n; ++i)
01068         fm[i][i] = s.diagonal()[i];
01069 }
01070 /* @} */
01071 } // end namespace
01072 #endif