dune-istl
2.2.0
|
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*- 00002 // vi: set et ts=8 sw=4 sts=4: 00003 #ifndef DUNE_SCALED_IDENTITY_MATRIX_HH 00004 #define DUNE_SCALED_IDENTITY_MATRIX_HH 00005 00012 #include<cmath> 00013 #include<cstddef> 00014 #include<complex> 00015 #include<iostream> 00016 #include <dune/common/exceptions.hh> 00017 #include <dune/common/fmatrix.hh> 00018 #include <dune/istl/diagonalmatrix.hh> 00019 00020 namespace Dune { 00021 00025 template<class K, int n> 00026 class ScaledIdentityMatrix 00027 { 00028 typedef DiagonalMatrixWrapper< ScaledIdentityMatrix<K,n> > WrapperType; 00029 00030 public: 00031 //===== type definitions and constants 00032 00034 typedef K field_type; 00035 00037 typedef K block_type; 00038 00040 typedef std::size_t size_type; 00041 00043 enum { 00045 blocklevel = 1 00046 }; 00047 00049 typedef DiagonalRowVector<K,n> row_type; 00050 typedef row_type reference; 00051 typedef DiagonalRowVectorConst<K,n> const_row_type; 00052 typedef const_row_type const_reference; 00053 00055 enum { 00057 rows = n, 00059 cols = n 00060 }; 00061 00062 //===== constructors 00065 ScaledIdentityMatrix () {} 00066 00069 ScaledIdentityMatrix (const K& k) 00070 : p_(k) 00071 {} 00072 00073 //===== assignment from scalar 00074 ScaledIdentityMatrix& operator= (const K& k) 00075 { 00076 p_ = k; 00077 return *this; 00078 } 00079 00080 // check if matrix is identical to other matrix (not only identical values) 00081 bool identical(const ScaledIdentityMatrix<K,n>& other) const 00082 { 00083 return (this==&other); 00084 } 00085 00086 //===== iterator interface to rows of the matrix 00088 typedef ContainerWrapperIterator<const WrapperType, reference, reference> Iterator; 00090 typedef Iterator iterator; 00092 typedef Iterator RowIterator; 00094 typedef typename row_type::Iterator ColIterator; 00095 00097 Iterator begin () 00098 { 00099 return Iterator(WrapperType(this),0); 00100 } 00101 00103 Iterator end () 00104 { 00105 return Iterator(WrapperType(this),n); 00106 } 00107 00110 Iterator beforeEnd () 00111 { 00112 return Iterator(WrapperType(this),n-1); 00113 } 00114 00117 Iterator beforeBegin () 00118 { 00119 return Iterator(WrapperType(this),-1); 00120 } 00121 00122 00124 typedef ContainerWrapperIterator<const WrapperType, const_reference, const_reference> ConstIterator; 00126 typedef ConstIterator const_iterator; 00128 typedef ConstIterator ConstRowIterator; 00130 typedef typename const_row_type::ConstIterator ConstColIterator; 00131 00133 ConstIterator begin () const 00134 { 00135 return ConstIterator(WrapperType(this),0); 00136 } 00137 00139 ConstIterator end () const 00140 { 00141 return ConstIterator(WrapperType(this),n); 00142 } 00143 00146 ConstIterator beforeEnd() const 00147 { 00148 return ConstIterator(WrapperType(this),n-1); 00149 } 00150 00153 ConstIterator beforeBegin () const 00154 { 00155 return ConstIterator(WrapperType(this),-1); 00156 } 00157 00158 //===== vector space arithmetic 00159 00161 ScaledIdentityMatrix& operator+= (const ScaledIdentityMatrix& y) 00162 { 00163 p_ += y.p_; 00164 return *this; 00165 } 00166 00168 ScaledIdentityMatrix& operator-= (const ScaledIdentityMatrix& y) 00169 { 00170 p_ -= y.p_; 00171 return *this; 00172 } 00173 00175 ScaledIdentityMatrix& operator+= (const K& k) 00176 { 00177 p_ += k; 00178 return *this; 00179 } 00180 00182 ScaledIdentityMatrix& operator-= (const K& k) 00183 { 00184 p_ -= k; 00185 return *this; 00186 } 00188 ScaledIdentityMatrix& operator*= (const K& k) 00189 { 00190 p_ *= k; 00191 return *this; 00192 } 00193 00195 ScaledIdentityMatrix& operator/= (const K& k) 00196 { 00197 p_ /= k; 00198 return *this; 00199 } 00200 00201 //===== comparison ops 00202 00204 bool operator==(const ScaledIdentityMatrix& other) const 00205 { 00206 return p_==other.scalar(); 00207 } 00208 00210 bool operator!=(const ScaledIdentityMatrix& other) const 00211 { 00212 return p_!=other.scalar(); 00213 } 00214 00215 //===== linear maps 00216 00218 template<class X, class Y> 00219 void mv (const X& x, Y& y) const 00220 { 00221 #ifdef DUNE_FMatrix_WITH_CHECKING 00222 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range"); 00223 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range"); 00224 #endif 00225 for (size_type i=0; i<n; ++i) 00226 y[i] = p_ * x[i]; 00227 } 00228 00230 template<class X, class Y> 00231 void mtv (const X& x, Y& y) const 00232 { 00233 mv(x, y); 00234 } 00235 00237 template<class X, class Y> 00238 void umv (const X& x, Y& y) const 00239 { 00240 #ifdef DUNE_FMatrix_WITH_CHECKING 00241 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range"); 00242 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range"); 00243 #endif 00244 for (size_type i=0; i<n; ++i) 00245 y[i] += p_ * x[i]; 00246 } 00247 00249 template<class X, class Y> 00250 void umtv (const X& x, Y& y) const 00251 { 00252 #ifdef DUNE_FMatrix_WITH_CHECKING 00253 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range"); 00254 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range"); 00255 #endif 00256 for (size_type i=0; i<n; ++i) 00257 y[i] += p_ * x[i]; 00258 } 00259 00261 template<class X, class Y> 00262 void umhv (const X& x, Y& y) const 00263 { 00264 #ifdef DUNE_FMatrix_WITH_CHECKING 00265 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range"); 00266 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range"); 00267 #endif 00268 for (size_type i=0; i<n; i++) 00269 y[i] += conjugateComplex(p_)*x[i]; 00270 } 00271 00273 template<class X, class Y> 00274 void mmv (const X& x, Y& y) const 00275 { 00276 #ifdef DUNE_FMatrix_WITH_CHECKING 00277 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range"); 00278 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range"); 00279 #endif 00280 for (size_type i=0; i<n; ++i) 00281 y[i] -= p_ * x[i]; 00282 } 00283 00285 template<class X, class Y> 00286 void mmtv (const X& x, Y& y) const 00287 { 00288 #ifdef DUNE_FMatrix_WITH_CHECKING 00289 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range"); 00290 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range"); 00291 #endif 00292 for (size_type i=0; i<n; ++i) 00293 y[i] -= p_ * x[i]; 00294 } 00295 00297 template<class X, class Y> 00298 void mmhv (const X& x, Y& y) const 00299 { 00300 #ifdef DUNE_FMatrix_WITH_CHECKING 00301 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range"); 00302 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range"); 00303 #endif 00304 for (size_type i=0; i<n; i++) 00305 y[i] -= conjugateComplex(p_)*x[i]; 00306 } 00307 00309 template<class X, class Y> 00310 void usmv (const K& alpha, const X& x, Y& y) const 00311 { 00312 #ifdef DUNE_FMatrix_WITH_CHECKING 00313 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range"); 00314 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range"); 00315 #endif 00316 for (size_type i=0; i<n; i++) 00317 y[i] += alpha * p_ * x[i]; 00318 } 00319 00321 template<class X, class Y> 00322 void usmtv (const K& alpha, const X& x, Y& y) const 00323 { 00324 #ifdef DUNE_FMatrix_WITH_CHECKING 00325 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range"); 00326 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range"); 00327 #endif 00328 for (size_type i=0; i<n; i++) 00329 y[i] += alpha * p_ * x[i]; 00330 } 00331 00333 template<class X, class Y> 00334 void usmhv (const K& alpha, const X& x, Y& y) const 00335 { 00336 #ifdef DUNE_FMatrix_WITH_CHECKING 00337 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range"); 00338 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range"); 00339 #endif 00340 for (size_type i=0; i<n; i++) 00341 y[i] += alpha * conjugateComplex(p_) * x[i]; 00342 } 00343 00344 //===== norms 00345 00347 double frobenius_norm () const 00348 { 00349 return fvmeta::sqrt(n*p_*p_); 00350 } 00351 00353 double frobenius_norm2 () const 00354 { 00355 return n*p_*p_; 00356 } 00357 00359 double infinity_norm () const 00360 { 00361 return std::abs(p_); 00362 } 00363 00365 double infinity_norm_real () const 00366 { 00367 return fvmeta::absreal(p_); 00368 } 00369 00370 //===== solve 00371 00374 template<class V> 00375 void solve (V& x, const V& b) const 00376 { 00377 for (int i=0; i<n; i++) 00378 x[i] = b[i]/p_; 00379 } 00380 00383 void invert() 00384 { 00385 p_ = 1/p_; 00386 } 00387 00389 K determinant () const { 00390 return std::pow(p_,n); 00391 } 00392 00393 //===== sizes 00394 00396 size_type N () const 00397 { 00398 return n; 00399 } 00400 00402 size_type M () const 00403 { 00404 return n; 00405 } 00406 00407 //===== query 00408 00410 bool exists (size_type i, size_type j) const 00411 { 00412 #ifdef DUNE_FMatrix_WITH_CHECKING 00413 if (i<0 || i>=n) DUNE_THROW(FMatrixError,"row index out of range"); 00414 if (j<0 || j>=n) DUNE_THROW(FMatrixError,"column index out of range"); 00415 #endif 00416 return i==j; 00417 } 00418 00419 //===== conversion operator 00420 00422 friend std::ostream& operator<< (std::ostream& s, const ScaledIdentityMatrix<K,n>& a) 00423 { 00424 for (size_type i=0; i<n; i++) { 00425 for (size_type j=0; j<n; j++) 00426 s << ((i==j) ? a.p_ : 0) << " "; 00427 s << std::endl; 00428 } 00429 return s; 00430 } 00431 00433 reference operator[](size_type i) 00434 { 00435 return reference(const_cast<K*>(&p_), i); 00436 } 00437 00439 const_reference operator[](size_type i) const 00440 { 00441 return const_reference(const_cast<K*>(&p_), i); 00442 } 00443 00445 const K& diagonal(size_type i) const 00446 { 00447 return p_; 00448 } 00449 00451 K& diagonal(size_type i) 00452 { 00453 return p_; 00454 } 00455 00458 const K& scalar() const 00459 { 00460 return p_; 00461 } 00462 00465 K& scalar() 00466 { 00467 return p_; 00468 } 00469 00470 private: 00471 // the data, very simply a single number 00472 K p_; 00473 00474 }; 00475 00476 template<class M, class K, int n> 00477 void istl_assign_to_fmatrix(DenseMatrix<M>& fm, const ScaledIdentityMatrix<K,n>& s) 00478 { 00479 fm = K(); 00480 for(int i=0; i<n; ++i) 00481 fm[i][i] = s.scalar(); 00482 } 00483 00484 } // end namespace 00485 00486 #endif