00001
00002
00003 #ifndef DUNE_FMATRIX_HH
00004 #define DUNE_FMATRIX_HH
00005
00006 #include <cmath>
00007 #include <cstddef>
00008 #include <iostream>
00009 #include <algorithm>
00010 #include <initializer_list>
00011
00012 #include <dune/common/boundschecking.hh>
00013 #include <dune/common/exceptions.hh>
00014 #include <dune/common/fvector.hh>
00015 #include <dune/common/densematrix.hh>
00016 #include <dune/common/precision.hh>
00017 #include <dune/common/typetraits.hh>
00018
00019 namespace Dune
00020 {
00021
00033 template< class K, int ROWS, int COLS > class FieldMatrix;
00034
00035 template< class K, int ROWS, int COLS >
00036 struct DenseMatVecTraits< FieldMatrix<K,ROWS,COLS> >
00037 {
00038 typedef FieldMatrix<K,ROWS,COLS> derived_type;
00039
00040
00041 typedef FieldVector<K,COLS> row_type;
00042
00043 typedef row_type &row_reference;
00044 typedef const row_type &const_row_reference;
00045
00046 typedef std::array<row_type,ROWS> container_type;
00047 typedef K value_type;
00048 typedef typename container_type::size_type size_type;
00049 };
00050
00051 template< class K, int ROWS, int COLS >
00052 struct FieldTraits< FieldMatrix<K,ROWS,COLS> >
00053 {
00054 typedef typename FieldTraits<K>::field_type field_type;
00055 typedef typename FieldTraits<K>::real_type real_type;
00056 };
00057
00066 template<class K, int ROWS, int COLS>
00067 class FieldMatrix : public DenseMatrix< FieldMatrix<K,ROWS,COLS> >
00068 {
00069 std::array< FieldVector<K,COLS>, ROWS > _data;
00070 typedef DenseMatrix< FieldMatrix<K,ROWS,COLS> > Base;
00071 public:
00072
00074 enum {
00076 rows = ROWS,
00078 cols = COLS
00079 };
00080
00081 typedef typename Base::size_type size_type;
00082 typedef typename Base::row_type row_type;
00083
00084 typedef typename Base::row_reference row_reference;
00085 typedef typename Base::const_row_reference const_row_reference;
00086
00087
00090 FieldMatrix () {}
00091
00094 FieldMatrix(std::initializer_list<Dune::FieldVector<K, cols> > const &l) {
00095 assert(l.size() == rows);
00096 std::copy_n(l.begin(), std::min(static_cast<std::size_t>(ROWS),
00097 l.size()),
00098 _data.begin());
00099 }
00100
00101 template <class T,
00102 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
00103 FieldMatrix(T const& rhs)
00104 {
00105 *this = rhs;
00106 }
00107
00108 using Base::operator=;
00109
00110
00111 template <typename T, int rows, int cols>
00112 FieldMatrix& operator=(FieldMatrix<T,rows,cols> const &rhs)
00113 {
00114 static_assert(rows == ROWS, "Size mismatch in matrix assignment (rows)");
00115 static_assert(cols == COLS, "Size mismatch in matrix assignment (columns)");
00116 _data = rhs._data;
00117 return *this;
00118 }
00119
00121 template<int l>
00122 FieldMatrix<K,l,cols> leftmultiplyany (const FieldMatrix<K,l,rows>& M) const
00123 {
00124 FieldMatrix<K,l,cols> C;
00125
00126 for (size_type i=0; i<l; i++) {
00127 for (size_type j=0; j<cols; j++) {
00128 C[i][j] = 0;
00129 for (size_type k=0; k<rows; k++)
00130 C[i][j] += M[i][k]*(*this)[k][j];
00131 }
00132 }
00133 return C;
00134 }
00135
00136 using Base::rightmultiply;
00137
00139 template <int r, int c>
00140 FieldMatrix& rightmultiply (const FieldMatrix<K,r,c>& M)
00141 {
00142 static_assert(r == c, "Cannot rightmultiply with non-square matrix");
00143 static_assert(r == cols, "Size mismatch");
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
00172 constexpr size_type mat_rows() const { return ROWS; }
00173 constexpr size_type mat_cols() const { return COLS; }
00174
00175 row_reference mat_access ( size_type i )
00176 {
00177 DUNE_ASSERT_BOUNDS(i < ROWS);
00178 return _data[i];
00179 }
00180
00181 const_row_reference mat_access ( size_type i ) const
00182 {
00183 DUNE_ASSERT_BOUNDS(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
00198
00199
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
00229 FieldMatrix () {}
00230
00233 FieldMatrix(std::initializer_list<Dune::FieldVector<K, 1>> const &l)
00234 {
00235 std::copy_n(l.begin(), std::min(static_cast< std::size_t >( 1 ), l.size()), &_data);
00236 }
00237
00238 template <class T,
00239 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
00240 FieldMatrix(T const& rhs)
00241 {
00242 *this = rhs;
00243 }
00244
00245 using Base::operator=;
00246
00247
00248
00250 template<int l>
00251 FieldMatrix<K,l,1> leftmultiplyany (const FieldMatrix<K,l,1>& M) const
00252 {
00253 FieldMatrix<K,l,1> C;
00254 for (size_type j=0; j<l; j++)
00255 C[j][0] = M[j][0]*(*this)[0][0];
00256 return C;
00257 }
00258
00260 FieldMatrix& rightmultiply (const FieldMatrix& M)
00261 {
00262 _data[0] *= M[0][0];
00263 return *this;
00264 }
00265
00267 template<int l>
00268 FieldMatrix<K,1,l> rightmultiplyany (const FieldMatrix<K,1,l>& M) const
00269 {
00270 FieldMatrix<K,1,l> C;
00271
00272 for (size_type j=0; j<l; j++)
00273 C[0][j] = M[0][j]*_data[0];
00274 return C;
00275 }
00276
00277
00278 constexpr size_type mat_rows() const { return 1; }
00279 constexpr size_type mat_cols() const { return 1; }
00280
00281 row_reference mat_access ( size_type i )
00282 {
00283 DUNE_UNUSED_PARAMETER(i);
00284 DUNE_ASSERT_BOUNDS(i == 0);
00285 return _data;
00286 }
00287
00288 const_row_reference mat_access ( size_type i ) const
00289 {
00290 DUNE_UNUSED_PARAMETER(i);
00291 DUNE_ASSERT_BOUNDS(i == 0);
00292 return _data;
00293 }
00294
00296 FieldMatrix& operator+= (const K& k)
00297 {
00298 _data[0] += k;
00299 return (*this);
00300 }
00301
00303 FieldMatrix& operator-= (const K& k)
00304 {
00305 _data[0] -= k;
00306 return (*this);
00307 }
00308
00310 FieldMatrix& operator*= (const K& k)
00311 {
00312 _data[0] *= k;
00313 return (*this);
00314 }
00315
00317 FieldMatrix& operator/= (const K& k)
00318 {
00319 _data[0] /= k;
00320 return (*this);
00321 }
00322
00323
00324
00325 operator const K& () const { return _data[0]; }
00326
00327 };
00328
00330 template<typename K>
00331 std::ostream& operator<< (std::ostream& s, const FieldMatrix<K,1,1>& a)
00332 {
00333 s << a[0][0];
00334 return s;
00335 }
00336
00337 #endif // DOXYGEN
00338
00339 namespace FMatrixHelp {
00340
00342 template <typename K>
00343 static inline K invertMatrix (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
00344 {
00345 inverse[0][0] = 1.0/matrix[0][0];
00346 return matrix[0][0];
00347 }
00348
00350 template <typename K>
00351 static inline K invertMatrix_retTransposed (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
00352 {
00353 return invertMatrix(matrix,inverse);
00354 }
00355
00356
00358 template <typename K>
00359 static inline K invertMatrix (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
00360 {
00361
00362 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
00363 K det_1 = 1.0/det;
00364 inverse[0][0] = matrix[1][1] * det_1;
00365 inverse[0][1] = - matrix[0][1] * det_1;
00366 inverse[1][0] = - matrix[1][0] * det_1;
00367 inverse[1][1] = matrix[0][0] * det_1;
00368 return det;
00369 }
00370
00373 template <typename K>
00374 static inline K invertMatrix_retTransposed (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
00375 {
00376
00377 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
00378 K det_1 = 1.0/det;
00379 inverse[0][0] = matrix[1][1] * det_1;
00380 inverse[1][0] = - matrix[0][1] * det_1;
00381 inverse[0][1] = - matrix[1][0] * det_1;
00382 inverse[1][1] = matrix[0][0] * det_1;
00383 return det;
00384 }
00385
00387 template <typename K>
00388 static inline K invertMatrix (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
00389 {
00390
00391 K t4 = matrix[0][0] * matrix[1][1];
00392 K t6 = matrix[0][0] * matrix[1][2];
00393 K t8 = matrix[0][1] * matrix[1][0];
00394 K t10 = matrix[0][2] * matrix[1][0];
00395 K t12 = matrix[0][1] * matrix[2][0];
00396 K t14 = matrix[0][2] * matrix[2][0];
00397
00398 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
00399 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
00400 K t17 = 1.0/det;
00401
00402 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
00403 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
00404 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
00405 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
00406 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
00407 inverse[1][2] = -(t6-t10) * t17;
00408 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
00409 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
00410 inverse[2][2] = (t4-t8) * t17;
00411
00412 return det;
00413 }
00414
00416 template <typename K>
00417 static inline K invertMatrix_retTransposed (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
00418 {
00419
00420 K t4 = matrix[0][0] * matrix[1][1];
00421 K t6 = matrix[0][0] * matrix[1][2];
00422 K t8 = matrix[0][1] * matrix[1][0];
00423 K t10 = matrix[0][2] * matrix[1][0];
00424 K t12 = matrix[0][1] * matrix[2][0];
00425 K t14 = matrix[0][2] * matrix[2][0];
00426
00427 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
00428 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
00429 K t17 = 1.0/det;
00430
00431 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
00432 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
00433 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
00434 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
00435 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
00436 inverse[2][1] = -(t6-t10) * t17;
00437 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
00438 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
00439 inverse[2][2] = (t4-t8) * t17;
00440
00441 return det;
00442 }
00443
00445 template< class K, int m, int n, int p >
00446 static inline void multMatrix ( const FieldMatrix< K, m, n > &A,
00447 const FieldMatrix< K, n, p > &B,
00448 FieldMatrix< K, m, p > &ret )
00449 {
00450 typedef typename FieldMatrix< K, m, p > :: size_type size_type;
00451
00452 for( size_type i = 0; i < m; ++i )
00453 {
00454 for( size_type j = 0; j < p; ++j )
00455 {
00456 ret[ i ][ j ] = K( 0 );
00457 for( size_type k = 0; k < n; ++k )
00458 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
00459 }
00460 }
00461 }
00462
00464 template <typename K, int rows, int cols>
00465 static inline void multTransposedMatrix(const FieldMatrix<K,rows,cols> &matrix, FieldMatrix<K,cols,cols>& ret)
00466 {
00467 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
00468
00469 for(size_type i=0; i<cols; i++)
00470 for(size_type j=0; j<cols; j++)
00471 {
00472 ret[i][j]=0.0;
00473 for(size_type k=0; k<rows; k++)
00474 ret[i][j]+=matrix[k][i]*matrix[k][j];
00475 }
00476 }
00477
00478 using Dune::DenseMatrixHelp::multAssign;
00479
00481 template <typename K, int rows, int cols>
00482 static inline void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret)
00483 {
00484 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
00485
00486 for(size_type i=0; i<cols; ++i)
00487 {
00488 ret[i] = 0.0;
00489 for(size_type j=0; j<rows; ++j)
00490 ret[i] += matrix[j][i]*x[j];
00491 }
00492 }
00493
00495 template <typename K, int rows, int cols>
00496 static inline FieldVector<K,rows> mult(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x)
00497 {
00498 FieldVector<K,rows> ret;
00499 multAssign(matrix,x,ret);
00500 return ret;
00501 }
00502
00504 template <typename K, int rows, int cols>
00505 static inline FieldVector<K,cols> multTransposed(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x)
00506 {
00507 FieldVector<K,cols> ret;
00508 multAssignTransposed( matrix, x, ret );
00509 return ret;
00510 }
00511
00512 }
00513
00516 }
00517
00518 #include "fmatrixev.hh"
00519 #endif