dune-istl  2.2.0
scaledidmatrix.hh
Go to the documentation of this file.
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