dune-common  2.2.0
fmatrix.hh
Go to the documentation of this file.
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=8 sw=2 sts=2:
00003 // $Id: fmatrix.hh 6785 2012-05-31 22:07:47Z sander $
00004 #ifndef DUNE_FMATRIX_HH
00005 #define DUNE_FMATRIX_HH
00006 
00007 #include <cmath>
00008 #include <cstddef>
00009 #include <iostream>
00010 
00011 #include <dune/common/misc.hh>
00012 #include <dune/common/exceptions.hh>
00013 #include <dune/common/fvector.hh>
00014 #include <dune/common/densematrix.hh>
00015 #include <dune/common/precision.hh>
00016 #include <dune/common/static_assert.hh>
00017 
00018 namespace Dune
00019 {
00020    
00032   template< class K, int ROWS, int COLS > class FieldMatrix;
00033 
00034   template< class K, int ROWS, int COLS >
00035   struct DenseMatVecTraits< FieldMatrix<K,ROWS,COLS> >
00036   {
00037     typedef FieldMatrix<K,ROWS,COLS> derived_type;
00038 
00039     // each row is implemented by a field vector
00040     typedef FieldVector<K,COLS> row_type;
00041 
00042     typedef row_type &row_reference;
00043     typedef const row_type &const_row_reference;
00044 
00045     typedef Dune::array<row_type,ROWS> container_type;
00046     typedef K value_type;
00047     typedef typename container_type::size_type size_type;
00048   };
00049   
00050   template< class K, int ROWS, int COLS >
00051   struct FieldTraits< FieldMatrix<K,ROWS,COLS> >
00052   {
00053     typedef typename FieldTraits<K>::field_type field_type;
00054     typedef typename FieldTraits<K>::real_type real_type;
00055   };
00056 
00065   template<class K, int ROWS, int COLS>
00066   class FieldMatrix : public DenseMatrix< FieldMatrix<K,ROWS,COLS> >
00067   {
00068     Dune::array< FieldVector<K,COLS>, ROWS > _data;
00069     typedef DenseMatrix< FieldMatrix<K,ROWS,COLS> > Base;
00070   public:
00071 
00073     enum {
00075       rows = ROWS, 
00077       cols = COLS
00078     };
00079 
00080     typedef typename Base::size_type size_type;
00081     typedef typename Base::row_type row_type;
00082    
00083     typedef typename Base::row_reference row_reference;
00084     typedef typename Base::const_row_reference const_row_reference;
00085 
00086     //===== constructors
00089     FieldMatrix () {}
00090 
00093     explicit FieldMatrix (const K& k)
00094     {
00095       for (size_type i=0; i<rows; i++) _data[i] = k;
00096     }
00097 
00098     template<typename T>
00099     explicit FieldMatrix (const T& t)
00100     {
00101       DenseMatrixAssigner<Conversion<T,K>::exists>::assign(*this, t);
00102     }
00103     
00104     //===== assignment
00105     using Base::operator=;
00106 
00107     // To be removed!
00108 #if 0
00109 
00110     FieldMatrix& leftmultiply (const FieldMatrix<K,rows,rows>& M)
00111     {
00112       FieldMatrix<K,rows,cols> C(*this);
00113       
00114       for (size_type i=0; i<rows; i++)
00115         for (size_type j=0; j<cols; j++) {
00116           (*this)[i][j] = 0;
00117           for (size_type k=0; k<rows; k++)
00118             (*this)[i][j] += M[i][k]*C[k][j];
00119         }
00120       
00121       return *this;
00122     }
00123 #endif
00124  
00126     template<int l>
00127     FieldMatrix<K,l,cols> leftmultiplyany (const FieldMatrix<K,l,rows>& M) const
00128     {
00129       FieldMatrix<K,l,cols> C;
00130       
00131       for (size_type i=0; i<l; i++) {
00132         for (size_type j=0; j<cols; j++) {
00133           C[i][j] = 0;
00134           for (size_type k=0; k<rows; k++)
00135             C[i][j] += M[i][k]*(*this)[k][j];
00136         }
00137       }
00138       return C;
00139     }
00140 
00142     FieldMatrix& rightmultiply (const FieldMatrix<K,cols,cols>& M)
00143     {
00144       FieldMatrix<K,rows,cols> C(*this);
00145       
00146       for (size_type i=0; i<rows; i++)
00147         for (size_type j=0; j<cols; j++) {
00148           (*this)[i][j] = 0;
00149           for (size_type k=0; k<cols; k++)
00150             (*this)[i][j] += C[i][k]*M[k][j];
00151         }
00152       return *this;
00153     }
00154 
00156     template<int l>
00157     FieldMatrix<K,rows,l> rightmultiplyany (const FieldMatrix<K,cols,l>& M) const
00158     {
00159       FieldMatrix<K,rows,l> C;
00160       
00161       for (size_type i=0; i<rows; i++) {
00162         for (size_type j=0; j<l; j++) {
00163           C[i][j] = 0;
00164           for (size_type k=0; k<cols; k++)
00165             C[i][j] += (*this)[i][k]*M[k][j];
00166         }
00167       }
00168       return C;
00169     }
00170     
00171     // make this thing a matrix
00172     size_type mat_rows() const { return ROWS; }
00173     size_type mat_cols() const { return COLS; }
00174 
00175     row_reference mat_access ( size_type i )
00176     {
00177       assert(i < ROWS);
00178       return _data[i];
00179     }
00180 
00181     const_row_reference mat_access ( size_type i ) const
00182     {
00183       assert(i < ROWS);
00184       return _data[i];
00185     }
00186   };
00187 
00188 #ifndef DOXYGEN // hide specialization
00189 
00191   template<class K>
00192   class FieldMatrix<K,1,1> : public DenseMatrix< FieldMatrix<K,1,1> >
00193   {
00194     FieldVector<K,1> _data;
00195     typedef DenseMatrix< FieldMatrix<K,1,1> > Base;
00196   public:
00197     // standard constructor and everything is sufficient ...
00198     
00199     //===== type definitions and constants
00200     
00202     typedef typename Base::size_type size_type;
00203     
00205     enum {
00208       blocklevel = 1
00209     };
00210 
00211     typedef typename Base::row_type row_type;
00212    
00213     typedef typename Base::row_reference row_reference;
00214     typedef typename Base::const_row_reference const_row_reference;
00215 
00217     enum {
00220       rows = 1,
00223       cols = 1
00224     };
00225 
00226     //===== constructors
00229     FieldMatrix () {}
00230     
00233     FieldMatrix (const K& k)
00234     {
00235       _data[0] = k;
00236     }
00237     template<typename T>
00238     FieldMatrix(const T& t)
00239     {
00240       DenseMatrixAssigner<Conversion<T,K>::exists>::assign(*this, t);
00241     }
00242     
00243     //===== solve
00244 
00246     template<int l>
00247     FieldMatrix<K,l,1> leftmultiplyany (const FieldMatrix<K,l,1>& M) const
00248     {
00249       FieldMatrix<K,l,1> C;
00250       for (size_type j=0; j<l; j++)
00251         C[j][0] = M[j][0]*(*this)[0][0];
00252       return C;
00253     }
00254 
00256     FieldMatrix& rightmultiply (const FieldMatrix& M)
00257     {
00258       _data[0] *= M[0][0];
00259       return *this;
00260     }
00261 
00263     template<int l>
00264     FieldMatrix<K,1,l> rightmultiplyany (const FieldMatrix<K,1,l>& M) const
00265     {
00266       FieldMatrix<K,1,l> C;
00267       
00268       for (size_type j=0; j<l; j++)
00269         C[0][j] = M[0][j]*_data[0];
00270       return C;
00271     }
00272 
00273     // make this thing a matrix
00274     size_type mat_rows() const { return 1; }
00275     size_type mat_cols() const { return 1; }
00276 
00277     row_reference mat_access ( size_type i )
00278     {
00279       assert(i == 0);
00280       return _data;
00281     }
00282 
00283     const_row_reference mat_access ( size_type i ) const
00284     {
00285       assert(i == 0);
00286       return _data;
00287     }
00288 
00290     FieldMatrix& operator+= (const K& k)
00291     {
00292       _data[0] += k;
00293       return (*this);
00294     }
00295 
00297     FieldMatrix& operator-= (const K& k)
00298     {
00299       _data[0] -= k;
00300       return (*this);
00301     }
00302     
00304     FieldMatrix& operator*= (const K& k)
00305     {
00306       _data[0] *= k;
00307       return (*this);
00308     }
00309 
00311     FieldMatrix& operator/= (const K& k)
00312     {
00313       _data[0] /= k;
00314       return (*this);
00315     }
00316 
00317     //===== conversion operator
00318 
00319     operator K () const { return _data[0]; }
00320 
00321   };
00322 
00324   template<typename K>
00325   std::ostream& operator<< (std::ostream& s, const FieldMatrix<K,1,1>& a)
00326   {
00327     s << a[0][0];
00328     return s;
00329   }
00330 
00331 #endif // DOXYGEN
00332 
00333 namespace FMatrixHelp {
00334 
00336 template <typename K>
00337 static inline K invertMatrix (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
00338 {
00339   inverse[0][0] = 1.0/matrix[0][0];
00340   return matrix[0][0];
00341 }
00342 
00344 template <typename K>
00345 static inline K invertMatrix_retTransposed (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
00346 {
00347   return invertMatrix(matrix,inverse); 
00348 }
00349 
00350 
00352 template <typename K>
00353 static inline K invertMatrix (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
00354 {
00355   // code generated by maple 
00356   K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
00357   K det_1 = 1.0/det;
00358   inverse[0][0] =   matrix[1][1] * det_1;
00359   inverse[0][1] = - matrix[0][1] * det_1;
00360   inverse[1][0] = - matrix[1][0] * det_1;
00361   inverse[1][1] =   matrix[0][0] * det_1;
00362   return det;
00363 }
00364 
00367 template <typename K>
00368 static inline K invertMatrix_retTransposed (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
00369 {
00370   // code generated by maple 
00371   K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
00372   K det_1 = 1.0/det;
00373   inverse[0][0] =   matrix[1][1] * det_1;
00374   inverse[1][0] = - matrix[0][1] * det_1;
00375   inverse[0][1] = - matrix[1][0] * det_1;
00376   inverse[1][1] =   matrix[0][0] * det_1;
00377   return det;
00378 }
00379 
00381 template <typename K>
00382 static inline K invertMatrix (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
00383 {
00384   // code generated by maple 
00385   K t4  = matrix[0][0] * matrix[1][1];
00386   K t6  = matrix[0][0] * matrix[1][2];
00387   K t8  = matrix[0][1] * matrix[1][0];
00388   K t10 = matrix[0][2] * matrix[1][0];
00389   K t12 = matrix[0][1] * matrix[2][0];
00390   K t14 = matrix[0][2] * matrix[2][0];
00391 
00392   K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
00393            t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
00394   K t17 = 1.0/det;
00395 
00396   inverse[0][0] =  (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
00397   inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
00398   inverse[0][2] =  (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
00399   inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
00400   inverse[1][1] =  (matrix[0][0] * matrix[2][2] - t14) * t17;
00401   inverse[1][2] = -(t6-t10) * t17;
00402   inverse[2][0] =  (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
00403   inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
00404   inverse[2][2] =  (t4-t8) * t17;
00405 
00406   return det;
00407 }
00408 
00410 template <typename K>
00411 static inline K invertMatrix_retTransposed (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
00412 {
00413   // code generated by maple 
00414   K t4  = matrix[0][0] * matrix[1][1];
00415   K t6  = matrix[0][0] * matrix[1][2];
00416   K t8  = matrix[0][1] * matrix[1][0];
00417   K t10 = matrix[0][2] * matrix[1][0];
00418   K t12 = matrix[0][1] * matrix[2][0];
00419   K t14 = matrix[0][2] * matrix[2][0];
00420 
00421   K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
00422            t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
00423   K t17 = 1.0/det;
00424 
00425   inverse[0][0] =  (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
00426   inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
00427   inverse[2][0] =  (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
00428   inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
00429   inverse[1][1] =  (matrix[0][0] * matrix[2][2] - t14) * t17;
00430   inverse[2][1] = -(t6-t10) * t17;
00431   inverse[0][2] =  (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
00432   inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
00433   inverse[2][2] =  (t4-t8) * t17;
00434 
00435   return det;
00436 }
00437 
00439 template< class K, int m, int n, int p >
00440 static inline void multMatrix ( const FieldMatrix< K, m, n > &A,
00441                                 const FieldMatrix< K, n, p > &B,
00442                                 FieldMatrix< K, m, p > &ret )
00443 {
00444   typedef typename FieldMatrix< K, m, p > :: size_type size_type;
00445 
00446   for( size_type i = 0; i < m; ++i )
00447   {
00448     for( size_type j = 0; j < p; ++j )
00449     {
00450       ret[ i ][ j ] = K( 0 );
00451       for( size_type k = 0; k < n; ++k )
00452         ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
00453     }
00454   }
00455 }
00456 
00458 template <typename K, int rows, int cols>
00459 static inline void multTransposedMatrix(const FieldMatrix<K,rows,cols> &matrix, FieldMatrix<K,cols,cols>& ret)
00460 {
00461   typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
00462   
00463   for(size_type i=0; i<cols; i++)
00464     for(size_type j=0; j<cols; j++)
00465     {
00466       ret[i][j]=0.0;
00467       for(size_type k=0; k<rows; k++)
00468         ret[i][j]+=matrix[k][i]*matrix[k][j];
00469     }
00470 }
00471 
00472 #if 0
00473 
00474 template <typename K, int rows, int cols>
00475 static inline void multAssign(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x, FieldVector<K,rows> & ret) 
00476 {
00477   typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
00478   
00479   for(size_type i=0; i<rows; ++i)
00480   {
00481     ret[i] = 0.0;
00482     for(size_type j=0; j<cols; ++j)
00483     {
00484       ret[i] += matrix[i][j]*x[j];
00485     }
00486   }
00487 }
00488 #else
00489 using Dune::DenseMatrixHelp::multAssign;
00490 #endif
00491 
00493 template <typename K, int rows, int cols>
00494 static inline void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret) 
00495 {
00496   typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
00497   
00498   for(size_type i=0; i<cols; ++i)
00499   {
00500     ret[i] = 0.0;
00501     for(size_type j=0; j<rows; ++j)
00502       ret[i] += matrix[j][i]*x[j];
00503   }
00504 }
00505 
00507 template <typename K, int rows, int cols>
00508 static inline FieldVector<K,rows> mult(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x) 
00509 {
00510   FieldVector<K,rows> ret;
00511   multAssign(matrix,x,ret);
00512   return ret;
00513 }
00514 
00516 template <typename K, int rows, int cols>
00517 static inline FieldVector<K,cols> multTransposed(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x) 
00518 {
00519   FieldVector<K,cols> ret;
00520   multAssignTransposed( matrix, x, ret );
00521   return ret; 
00522 }
00523 
00524 } // end namespace FMatrixHelp 
00525 
00528 } // end namespace
00529 
00530 #include "fmatrixev.hh"
00531 #endif