dune-istl
2.2.0
|
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