dune-common  2.2.0
fmatrix.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 // $Id: fmatrix.hh 6785 2012-05-31 22:07:47Z sander $
4 #ifndef DUNE_FMATRIX_HH
5 #define DUNE_FMATRIX_HH
6 
7 #include <cmath>
8 #include <cstddef>
9 #include <iostream>
10 
11 #include <dune/common/misc.hh>
13 #include <dune/common/fvector.hh>
15 #include <dune/common/precision.hh>
17 
18 namespace 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 
43  typedef const row_type &const_row_reference;
44 
46  typedef K value_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  {
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  explicit FieldMatrix (const K& k)
94  {
95  for (size_type i=0; i<rows; i++) _data[i] = k;
96  }
97 
98  template<typename T>
99  explicit FieldMatrix (const T& t)
100  {
101  DenseMatrixAssigner<Conversion<T,K>::exists>::assign(*this, t);
102  }
103 
104  //===== assignment
105  using Base::operator=;
106 
107  // To be removed!
108 #if 0
109 
111  {
112  FieldMatrix<K,rows,cols> C(*this);
113 
114  for (size_type i=0; i<rows; i++)
115  for (size_type j=0; j<cols; j++) {
116  (*this)[i][j] = 0;
117  for (size_type k=0; k<rows; k++)
118  (*this)[i][j] += M[i][k]*C[k][j];
119  }
120 
121  return *this;
122  }
123 #endif
124 
126  template<int l>
128  {
130 
131  for (size_type i=0; i<l; i++) {
132  for (size_type j=0; j<cols; j++) {
133  C[i][j] = 0;
134  for (size_type k=0; k<rows; k++)
135  C[i][j] += M[i][k]*(*this)[k][j];
136  }
137  }
138  return C;
139  }
140 
143  {
144  FieldMatrix<K,rows,cols> C(*this);
145 
146  for (size_type i=0; i<rows; i++)
147  for (size_type j=0; j<cols; j++) {
148  (*this)[i][j] = 0;
149  for (size_type k=0; k<cols; k++)
150  (*this)[i][j] += C[i][k]*M[k][j];
151  }
152  return *this;
153  }
154 
156  template<int l>
158  {
160 
161  for (size_type i=0; i<rows; i++) {
162  for (size_type j=0; j<l; j++) {
163  C[i][j] = 0;
164  for (size_type k=0; k<cols; k++)
165  C[i][j] += (*this)[i][k]*M[k][j];
166  }
167  }
168  return C;
169  }
170 
171  // make this thing a matrix
172  size_type mat_rows() const { return ROWS; }
173  size_type mat_cols() const { return COLS; }
174 
176  {
177  assert(i < ROWS);
178  return _data[i];
179  }
180 
182  {
183  assert(i < ROWS);
184  return _data[i];
185  }
186  };
187 
188 #ifndef DOXYGEN // hide specialization
189 
191  template<class K>
192  class FieldMatrix<K,1,1> : public DenseMatrix< FieldMatrix<K,1,1> >
193  {
194  FieldVector<K,1> _data;
195  typedef DenseMatrix< FieldMatrix<K,1,1> > Base;
196  public:
197  // standard constructor and everything is sufficient ...
198 
199  //===== type definitions and constants
200 
202  typedef typename Base::size_type size_type;
203 
205  enum {
208  blocklevel = 1
209  };
210 
211  typedef typename Base::row_type row_type;
212 
213  typedef typename Base::row_reference row_reference;
215 
217  enum {
220  rows = 1,
223  cols = 1
224  };
225 
226  //===== constructors
229  FieldMatrix () {}
230 
233  FieldMatrix (const K& k)
234  {
235  _data[0] = k;
236  }
237  template<typename T>
238  FieldMatrix(const T& t)
239  {
240  DenseMatrixAssigner<Conversion<T,K>::exists>::assign(*this, t);
241  }
242 
243  //===== solve
244 
246  template<int l>
247  FieldMatrix<K,l,1> leftmultiplyany (const FieldMatrix<K,l,1>& M) const
248  {
249  FieldMatrix<K,l,1> C;
250  for (size_type j=0; j<l; j++)
251  C[j][0] = M[j][0]*(*this)[0][0];
252  return C;
253  }
254 
257  {
258  _data[0] *= M[0][0];
259  return *this;
260  }
261 
263  template<int l>
264  FieldMatrix<K,1,l> rightmultiplyany (const FieldMatrix<K,1,l>& M) const
265  {
266  FieldMatrix<K,1,l> C;
267 
268  for (size_type j=0; j<l; j++)
269  C[0][j] = M[0][j]*_data[0];
270  return C;
271  }
272 
273  // make this thing a matrix
274  size_type mat_rows() const { return 1; }
275  size_type mat_cols() const { return 1; }
276 
278  {
279  assert(i == 0);
280  return _data;
281  }
282 
284  {
285  assert(i == 0);
286  return _data;
287  }
288 
290  FieldMatrix& operator+= (const K& k)
291  {
292  _data[0] += k;
293  return (*this);
294  }
295 
297  FieldMatrix& operator-= (const K& k)
298  {
299  _data[0] -= k;
300  return (*this);
301  }
302 
304  FieldMatrix& operator*= (const K& k)
305  {
306  _data[0] *= k;
307  return (*this);
308  }
309 
311  FieldMatrix& operator/= (const K& k)
312  {
313  _data[0] /= k;
314  return (*this);
315  }
316 
317  //===== conversion operator
318 
319  operator K () const { return _data[0]; }
320 
321  };
322 
324  template<typename K>
325  std::ostream& operator<< (std::ostream& s, const FieldMatrix<K,1,1>& a)
326  {
327  s << a[0][0];
328  return s;
329  }
330 
331 #endif // DOXYGEN
332 
333 namespace FMatrixHelp {
334 
336 template <typename K>
337 static inline K invertMatrix (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
338 {
339  inverse[0][0] = 1.0/matrix[0][0];
340  return matrix[0][0];
341 }
342 
344 template <typename K>
345 static inline K invertMatrix_retTransposed (const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
346 {
347  return invertMatrix(matrix,inverse);
348 }
349 
350 
352 template <typename K>
353 static inline K invertMatrix (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
354 {
355  // code generated by maple
356  K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
357  K det_1 = 1.0/det;
358  inverse[0][0] = matrix[1][1] * det_1;
359  inverse[0][1] = - matrix[0][1] * det_1;
360  inverse[1][0] = - matrix[1][0] * det_1;
361  inverse[1][1] = matrix[0][0] * det_1;
362  return det;
363 }
364 
367 template <typename K>
368 static inline K invertMatrix_retTransposed (const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
369 {
370  // code generated by maple
371  K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
372  K det_1 = 1.0/det;
373  inverse[0][0] = matrix[1][1] * det_1;
374  inverse[1][0] = - matrix[0][1] * det_1;
375  inverse[0][1] = - matrix[1][0] * det_1;
376  inverse[1][1] = matrix[0][0] * det_1;
377  return det;
378 }
379 
381 template <typename K>
382 static inline K invertMatrix (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
383 {
384  // code generated by maple
385  K t4 = matrix[0][0] * matrix[1][1];
386  K t6 = matrix[0][0] * matrix[1][2];
387  K t8 = matrix[0][1] * matrix[1][0];
388  K t10 = matrix[0][2] * matrix[1][0];
389  K t12 = matrix[0][1] * matrix[2][0];
390  K t14 = matrix[0][2] * matrix[2][0];
391 
392  K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
393  t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
394  K t17 = 1.0/det;
395 
396  inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
397  inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
398  inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
399  inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
400  inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
401  inverse[1][2] = -(t6-t10) * t17;
402  inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
403  inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
404  inverse[2][2] = (t4-t8) * t17;
405 
406  return det;
407 }
408 
410 template <typename K>
411 static inline K invertMatrix_retTransposed (const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
412 {
413  // code generated by maple
414  K t4 = matrix[0][0] * matrix[1][1];
415  K t6 = matrix[0][0] * matrix[1][2];
416  K t8 = matrix[0][1] * matrix[1][0];
417  K t10 = matrix[0][2] * matrix[1][0];
418  K t12 = matrix[0][1] * matrix[2][0];
419  K t14 = matrix[0][2] * matrix[2][0];
420 
421  K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
422  t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
423  K t17 = 1.0/det;
424 
425  inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
426  inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
427  inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
428  inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
429  inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
430  inverse[2][1] = -(t6-t10) * t17;
431  inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
432  inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
433  inverse[2][2] = (t4-t8) * t17;
434 
435  return det;
436 }
437 
439 template< class K, int m, int n, int p >
440 static inline void multMatrix ( const FieldMatrix< K, m, n > &A,
441  const FieldMatrix< K, n, p > &B,
443 {
444  typedef typename FieldMatrix< K, m, p > :: size_type size_type;
445 
446  for( size_type i = 0; i < m; ++i )
447  {
448  for( size_type j = 0; j < p; ++j )
449  {
450  ret[ i ][ j ] = K( 0 );
451  for( size_type k = 0; k < n; ++k )
452  ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
453  }
454  }
455 }
456 
458 template <typename K, int rows, int cols>
460 {
461  typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
462 
463  for(size_type i=0; i<cols; i++)
464  for(size_type j=0; j<cols; j++)
465  {
466  ret[i][j]=0.0;
467  for(size_type k=0; k<rows; k++)
468  ret[i][j]+=matrix[k][i]*matrix[k][j];
469  }
470 }
471 
472 #if 0
473 
474 template <typename K, int rows, int cols>
475 static inline void multAssign(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x, FieldVector<K,rows> & ret)
476 {
477  typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
478 
479  for(size_type i=0; i<rows; ++i)
480  {
481  ret[i] = 0.0;
482  for(size_type j=0; j<cols; ++j)
483  {
484  ret[i] += matrix[i][j]*x[j];
485  }
486  }
487 }
488 #else
490 #endif
491 
493 template <typename K, int rows, int cols>
494 static inline void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret)
495 {
496  typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
497 
498  for(size_type i=0; i<cols; ++i)
499  {
500  ret[i] = 0.0;
501  for(size_type j=0; j<rows; ++j)
502  ret[i] += matrix[j][i]*x[j];
503  }
504 }
505 
507 template <typename K, int rows, int cols>
509 {
511  multAssign(matrix,x,ret);
512  return ret;
513 }
514 
516 template <typename K, int rows, int cols>
518 {
520  multAssignTransposed( matrix, x, ret );
521  return ret;
522 }
523 
524 } // end namespace FMatrixHelp
525 
528 } // end namespace
529 
530 #include "fmatrixev.hh"
531 #endif