dune-common 3.0-git
fmatrix.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_FMATRIX_HH
4#define DUNE_FMATRIX_HH
5
6#include <cmath>
7#include <cstddef>
8#include <iostream>
9#include <algorithm>
10#include <initializer_list>
11
17
18namespace Dune
19{
20
32 template< class K, int ROWS, int COLS > class FieldMatrix;
33
34 template< class K, int ROWS, int COLS >
35 struct DenseMatVecTraits< FieldMatrix<K,ROWS,COLS> >
36 {
38
39 // each row is implemented by a field vector
41
44
45 typedef std::array<row_type,ROWS> container_type;
46 typedef K value_type;
47 typedef typename container_type::size_type size_type;
48 };
49
50 template< class K, int ROWS, int COLS >
51 struct FieldTraits< FieldMatrix<K,ROWS,COLS> >
52 {
55 };
56
65 template<class K, int ROWS, int COLS>
66 class FieldMatrix : public DenseMatrix< FieldMatrix<K,ROWS,COLS> >
67 {
68 std::array< FieldVector<K,COLS>, ROWS > _data;
70 public:
71
73 enum {
75 rows = ROWS,
77 cols = COLS
78 };
79
80 typedef typename Base::size_type size_type;
81 typedef typename Base::row_type row_type;
82
85
86 //===== constructors
90
93 template< class Other >
94 FieldMatrix ( const Other &other )
95 {
96 DenseMatrixAssigner< FieldMatrix< K, ROWS, COLS >, Other >::apply( *this, other );
97 }
98
101 FieldMatrix (std::initializer_list<std::initializer_list<K> > const &ll)
102 {
103 assert(ll.size() == rows); // Actually, this is not needed any more!
104 std::copy_n(ll.begin(), std::min(static_cast<std::size_t>(ROWS),
105 ll.size()),
106 _data.begin());
107 }
108
111 FieldMatrix(std::initializer_list<Dune::FieldVector<K, cols> > const &l) {
112 assert(l.size() == rows); // Actually, this is not needed any more!
113 std::copy_n(l.begin(), std::min(static_cast<std::size_t>(ROWS),
114 l.size()),
115 _data.begin());
116 }
117
118 //===== assignment
119 using Base::operator=;
120
122 template<int l>
124 {
126
127 for (size_type i=0; i<l; i++) {
128 for (size_type j=0; j<cols; j++) {
129 C[i][j] = 0;
130 for (size_type k=0; k<rows; k++)
131 C[i][j] += M[i][k]*(*this)[k][j];
132 }
133 }
134 return C;
135 }
136
138
140 template <int r, int c>
142 {
143 static_assert(r == c, "Cannot rightmultiply with non-square matrix");
144 static_assert(r == cols, "Size mismatch");
146
147 for (size_type i=0; i<rows; i++)
148 for (size_type j=0; j<cols; j++) {
149 (*this)[i][j] = 0;
150 for (size_type k=0; k<cols; k++)
151 (*this)[i][j] += C[i][k]*M[k][j];
152 }
153 return *this;
154 }
155
157 template<int l>
159 {
161
162 for (size_type i=0; i<rows; i++) {
163 for (size_type j=0; j<l; j++) {
164 C[i][j] = 0;
165 for (size_type k=0; k<cols; k++)
166 C[i][j] += (*this)[i][k]*M[k][j];
167 }
168 }
169 return C;
170 }
171
172 // make this thing a matrix
173 constexpr size_type mat_rows() const { return ROWS; }
174 constexpr size_type mat_cols() const { return COLS; }
175
177 {
178 DUNE_ASSERT_BOUNDS(i < ROWS);
179 return _data[i];
180 }
181
183 {
184 DUNE_ASSERT_BOUNDS(i < ROWS);
185 return _data[i];
186 }
187 };
188
189#ifndef DOXYGEN // hide specialization
192 template<class K>
193 class FieldMatrix<K,1,1> : public DenseMatrix< FieldMatrix<K,1,1> >
194 {
195 FieldVector<K,1> _data;
196 typedef DenseMatrix< FieldMatrix<K,1,1> > Base;
197 public:
198 // standard constructor and everything is sufficient ...
199
200 //===== type definitions and constants
201
203 typedef typename Base::size_type size_type;
204
206 enum {
209 blocklevel = 1
210 };
211
212 typedef typename Base::row_type row_type;
213
214 typedef typename Base::row_reference row_reference;
216
218 enum {
221 rows = 1,
224 cols = 1
225 };
226
227 //===== constructors
230 FieldMatrix () {}
231
234 FieldMatrix (const K& k)
235 {
236 _data[0] = k;
237 }
238
239 template< class Other >
240 FieldMatrix ( const Other &other )
241 {
242 DenseMatrixAssigner< FieldMatrix< K, 1, 1 >, Other >::apply( *this, other );
243 }
244
245 //===== solve
246
248 template<int l>
249 FieldMatrix<K,l,1> leftmultiplyany (const FieldMatrix<K,l,1>& M) const
250 {
251 FieldMatrix<K,l,1> C;
252 for (size_type j=0; j<l; j++)
253 C[j][0] = M[j][0]*(*this)[0][0];
254 return C;
255 }
256
259 {
260 _data[0] *= M[0][0];
261 return *this;
262 }
263
265 template<int l>
266 FieldMatrix<K,1,l> rightmultiplyany (const FieldMatrix<K,1,l>& M) const
267 {
268 FieldMatrix<K,1,l> C;
269
270 for (size_type j=0; j<l; j++)
271 C[0][j] = M[0][j]*_data[0];
272 return C;
273 }
274
275 // make this thing a matrix
276 constexpr size_type mat_rows() const { return 1; }
277 constexpr size_type mat_cols() const { return 1; }
278
280 {
282 DUNE_ASSERT_BOUNDS(i == 0);
283 return _data;
284 }
285
287 {
289 DUNE_ASSERT_BOUNDS(i == 0);
290 return _data;
291 }
292
294 FieldMatrix& operator+= (const K& k)
295 {
296 _data[0] += k;
297 return (*this);
298 }
299
301 FieldMatrix& operator-= (const K& k)
302 {
303 _data[0] -= k;
304 return (*this);
305 }
306
308 FieldMatrix& operator*= (const K& k)
309 {
310 _data[0] *= k;
311 return (*this);
312 }
313
315 FieldMatrix& operator/= (const K& k)
316 {
317 _data[0] /= k;
318 return (*this);
319 }
320
321 //===== conversion operator
322
323 operator const K& () const { return _data[0]; }
324
325 };
326
328 template<typename K>
329 std::ostream& operator<< (std::ostream& s, const FieldMatrix<K,1,1>& a)
330 {
331 s << a[0][0];
332 return s;
333 }
334
335#endif // DOXYGEN
336
337 namespace FMatrixHelp {
338
340 template <typename K>
341 static inline K invertMatrix (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
342 {
343 inverse[0][0] = 1.0/matrix[0][0];
344 return matrix[0][0];
345 }
346
348 template <typename K>
349 static inline K invertMatrix_retTransposed (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
350 {
351 return invertMatrix(matrix,inverse);
352 }
353
354
356 template <typename K>
357 static inline K invertMatrix (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
358 {
359 // code generated by maple
360 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
361 K det_1 = 1.0/det;
362 inverse[0][0] = matrix[1][1] * det_1;
363 inverse[0][1] = - matrix[0][1] * det_1;
364 inverse[1][0] = - matrix[1][0] * det_1;
365 inverse[1][1] = matrix[0][0] * det_1;
366 return det;
367 }
368
371 template <typename K>
372 static inline K invertMatrix_retTransposed (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
373 {
374 // code generated by maple
375 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
376 K det_1 = 1.0/det;
377 inverse[0][0] = matrix[1][1] * det_1;
378 inverse[1][0] = - matrix[0][1] * det_1;
379 inverse[0][1] = - matrix[1][0] * det_1;
380 inverse[1][1] = matrix[0][0] * det_1;
381 return det;
382 }
383
385 template <typename K>
386 static inline K invertMatrix (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
387 {
388 // code generated by maple
389 K t4 = matrix[0][0] * matrix[1][1];
390 K t6 = matrix[0][0] * matrix[1][2];
391 K t8 = matrix[0][1] * matrix[1][0];
392 K t10 = matrix[0][2] * matrix[1][0];
393 K t12 = matrix[0][1] * matrix[2][0];
394 K t14 = matrix[0][2] * matrix[2][0];
395
396 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
397 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
398 K t17 = 1.0/det;
399
400 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
401 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
402 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
403 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
404 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
405 inverse[1][2] = -(t6-t10) * t17;
406 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
407 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
408 inverse[2][2] = (t4-t8) * t17;
409
410 return det;
411 }
412
414 template <typename K>
415 static inline K invertMatrix_retTransposed (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
416 {
417 // code generated by maple
418 K t4 = matrix[0][0] * matrix[1][1];
419 K t6 = matrix[0][0] * matrix[1][2];
420 K t8 = matrix[0][1] * matrix[1][0];
421 K t10 = matrix[0][2] * matrix[1][0];
422 K t12 = matrix[0][1] * matrix[2][0];
423 K t14 = matrix[0][2] * matrix[2][0];
424
425 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
426 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
427 K t17 = 1.0/det;
428
429 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
430 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
431 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
432 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
433 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
434 inverse[2][1] = -(t6-t10) * t17;
435 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
436 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
437 inverse[2][2] = (t4-t8) * t17;
438
439 return det;
440 }
441
443 template< class K, int m, int n, int p >
444 static inline void multMatrix ( const FieldMatrix< K, m, n > &A,
445 const FieldMatrix< K, n, p > &B,
447 {
448 typedef typename FieldMatrix< K, m, p > :: size_type size_type;
449
450 for( size_type i = 0; i < m; ++i )
451 {
452 for( size_type j = 0; j < p; ++j )
453 {
454 ret[ i ][ j ] = K( 0 );
455 for( size_type k = 0; k < n; ++k )
456 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
457 }
458 }
459 }
460
462 template <typename K, int rows, int cols>
464 {
465 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
466
467 for(size_type i=0; i<cols; i++)
468 for(size_type j=0; j<cols; j++)
469 {
470 ret[i][j]=0.0;
471 for(size_type k=0; k<rows; k++)
472 ret[i][j]+=matrix[k][i]*matrix[k][j];
473 }
474 }
475
477
479 template <typename K, int rows, int cols>
481 {
482 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
483
484 for(size_type i=0; i<cols; ++i)
485 {
486 ret[i] = 0.0;
487 for(size_type j=0; j<rows; ++j)
488 ret[i] += matrix[j][i]*x[j];
489 }
490 }
491
493 template <typename K, int rows, int cols>
495 {
497 multAssign(matrix,x,ret);
498 return ret;
499 }
500
502 template <typename K, int rows, int cols>
504 {
506 multAssignTransposed( matrix, x, ret );
507 return ret;
508 }
509
510 } // end namespace FMatrixHelp
511
514} // end namespace
515
516#include "fmatrixev.hh"
517#endif
#define DUNE_ASSERT_BOUNDS(cond)
Definition boundschecking.hh:20
Implements a vector constructed from a given type representing a field and a compile-time given size.
Various precision settings for calculations with FieldMatrix and FieldVector.
Implements a matrix constructed from a given type representing a field and a compile-time given numbe...
A few common exception classes.
Eigenvalue computations for the FieldMatrix class.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition unused.hh:18
std::ostream & operator<<(std::ostream &s, const std::array< T, N > &e)
Output operator for array.
Definition array.hh:28
Dune namespace.
Definition alignment.hh:11
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition densematrix.hh:1147
static FieldVector< K, cols > multTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x)
calculates ret = matrix^T * x
Definition fmatrix.hh:503
static K invertMatrix_retTransposed(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition fmatrix.hh:349
static void multMatrix(const FieldMatrix< K, m, n > &A, const FieldMatrix< K, n, p > &B, FieldMatrix< K, m, p > &ret)
calculates ret = A * B
Definition fmatrix.hh:444
static K invertMatrix(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition fmatrix.hh:341
static FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition fmatrix.hh:494
static void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition fmatrix.hh:463
static void multAssignTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x, FieldVector< K, cols > &ret)
calculates ret = matrix^T * x
Definition fmatrix.hh:480
A dense n x m matrix.
Definition densematrix.hh:135
DenseMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition densematrix.hh:305
DenseMatrix & operator+=(const DenseMatrix< Other > &y)
vector space addition
Definition densematrix.hh:278
DenseMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition densematrix.hh:297
DenseMatrix & operator-=(const DenseMatrix< Other > &y)
vector space subtraction
Definition densematrix.hh:288
size_type M() const
number of columns
Definition densematrix.hh:658
FieldMatrix< K, ROWS, COLS > & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition densematrix.hh:600
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition densematrix.hh:161
Traits::size_type size_type
The type used for the index access and size operation.
Definition densematrix.hh:158
Traits::const_row_reference const_row_reference
The type used to represent a reference to a constant row (usually const row_type &)
Definition densematrix.hh:167
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition densematrix.hh:164
@ blocklevel
The number of block levels we contain. This is 1.
Definition densematrix.hh:172
A dense n x m matrix.
Definition fmatrix.hh:67
@ rows
The number of rows.
Definition fmatrix.hh:75
@ cols
The number of columns.
Definition fmatrix.hh:77
FieldMatrix(const Other &other)
Constructor initializing the whole matrix with a scalar.
Definition fmatrix.hh:94
Base::const_row_reference const_row_reference
Definition fmatrix.hh:84
FieldMatrix< K, rows, l > rightmultiplyany(const FieldMatrix< K, cols, l > &M) const
Multiplies M from the right to this matrix, this matrix is not modified.
Definition fmatrix.hh:158
Base::row_type row_type
Definition fmatrix.hh:81
Base::size_type size_type
Definition fmatrix.hh:80
FieldMatrix< K, l, cols > leftmultiplyany(const FieldMatrix< K, l, rows > &M) const
Multiplies M from the left to this matrix, this matrix is not modified.
Definition fmatrix.hh:123
FieldMatrix & rightmultiply(const FieldMatrix< K, r, c > &M)
Multiplies M from the right to this matrix.
Definition fmatrix.hh:141
FieldMatrix()
Default constructor.
Definition fmatrix.hh:89
Base::row_reference row_reference
Definition fmatrix.hh:83
constexpr size_type mat_cols() const
Definition fmatrix.hh:174
row_reference mat_access(size_type i)
Definition fmatrix.hh:176
FieldMatrix(std::initializer_list< std::initializer_list< K > > const &ll)
Constructor initializing the matrix from a list of lists of scalars.
Definition fmatrix.hh:101
FieldMatrix(std::initializer_list< Dune::FieldVector< K, cols > > const &l)
Constructor initializing the matrix from a list of vector.
Definition fmatrix.hh:111
constexpr size_type mat_rows() const
Definition fmatrix.hh:173
const_row_reference mat_access(size_type i) const
Definition fmatrix.hh:182
vector space out of a tensor product of fields.
Definition fvector.hh:93
you have to specialize this structure for any type that should be assignable to a DenseMatrix
Definition densematrix.hh:73
std::array< row_type, ROWS > container_type
Definition fmatrix.hh:45
FieldMatrix< K, ROWS, COLS > derived_type
Definition fmatrix.hh:37
const row_type & const_row_reference
Definition fmatrix.hh:43
FieldVector< K, COLS > row_type
Definition fmatrix.hh:40
container_type::size_type size_type
Definition fmatrix.hh:47
row_type & row_reference
Definition fmatrix.hh:42
FieldTraits< K >::field_type field_type
Definition fmatrix.hh:53
FieldTraits< K >::real_type real_type
Definition fmatrix.hh:54
Definition ftraits.hh:24
T field_type
export the type representing the field
Definition ftraits.hh:26
T real_type
export the type representing the real type of the field
Definition ftraits.hh:28
Definition matvectraits.hh:29