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