dune-geometry  2.3.0
matrixhelper.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_GEOMETRY_GENERICGEOMETRY_MATRIXHELPER_HH
4 #define DUNE_GEOMETRY_GENERICGEOMETRY_MATRIXHELPER_HH
5 
6 #include <cmath>
7 
8 #include <dune/common/fvector.hh>
9 #include <dune/common/fmatrix.hh>
10 #include <dune/common/static_assert.hh>
11 
12 namespace Dune
13 {
14 
15  namespace GenericGeometry
16  {
17 
18  // FieldHelper
19  // -----------
20 
21  template< class Field >
22  struct FieldHelper
23  {
24  static Field abs ( const Field &x ) { return std::abs( x ); }
25  };
26 
27 
28 
29  // MatrixHelper
30  // ------------
31 
32  template< class Traits >
33  struct MatrixHelper
34  {
35  typedef typename Traits::ctype FieldType;
36 
37  static FieldType abs ( const FieldType &x )
38  {
39  //return std::abs( x );
41  }
42 
43  template< int m, int n >
44  static void
45  Ax ( const typename Traits :: template Matrix< m, n > :: type &A,
46  const typename Traits :: template Vector< n > :: type &x,
47  typename Traits :: template Vector< m > :: type &ret )
48  {
49  for( int i = 0; i < m; ++i )
50  {
51  ret[ i ] = FieldType( 0 );
52  for( int j = 0; j < n; ++j )
53  ret[ i ] += A[ i ][ j ] * x[ j ];
54  }
55  }
56 
57  template< int m, int n >
58  static void
59  ATx ( const typename Traits :: template Matrix< m, n > :: type &A,
60  const typename Traits :: template Vector< m > :: type &x,
61  typename Traits :: template Vector< n > :: type &ret )
62  {
63  for( int i = 0; i < n; ++i )
64  {
65  ret[ i ] = FieldType( 0 );
66  for( int j = 0; j < m; ++j )
67  ret[ i ] += A[ j ][ i ] * x[ j ];
68  }
69  }
70 
71  template< int m, int n, int p >
72  static void
73  AB ( const typename Traits :: template Matrix< m, n > :: type &A,
74  const typename Traits :: template Matrix< n, p > :: type &B,
75  typename Traits :: template Matrix< m, p > :: type &ret )
76  {
77  for( int i = 0; i < m; ++i )
78  {
79  for( int j = 0; j < p; ++j )
80  {
81  ret[ i ][ j ] = FieldType( 0 );
82  for( int k = 0; k < n; ++k )
83  ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
84  }
85  }
86  }
87 
88  template< int m, int n, int p >
89  static void
90  ATBT ( const typename Traits :: template Matrix< m, n > :: type &A,
91  const typename Traits :: template Matrix< p, m > :: type &B,
92  typename Traits :: template Matrix< n, p > :: type &ret )
93  {
94  for( int i = 0; i < n; ++i )
95  {
96  for( int j = 0; j < p; ++j )
97  {
98  ret[ i ][ j ] = FieldType( 0 );
99  for( int k = 0; k < m; ++k )
100  ret[ i ][ j ] += A[ k ][ i ] * B[ j ][ k ];
101  }
102  }
103  }
104 
105  template< int m, int n >
106  static void
107  ATA_L ( const typename Traits :: template Matrix< m, n > :: type &A,
108  typename Traits :: template Matrix< n, n > :: type &ret )
109  {
110  for( int i = 0; i < n; ++i )
111  {
112  for( int j = 0; j <= i; ++j )
113  {
114  ret[ i ][ j ] = FieldType( 0 );
115  for( int k = 0; k < m; ++k )
116  ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
117  }
118  }
119  }
120 
121  template< int m, int n >
122  static void
123  ATA ( const typename Traits :: template Matrix< m, n > :: type &A,
124  typename Traits :: template Matrix< n, n > :: type &ret )
125  {
126  for( int i = 0; i < n; ++i )
127  {
128  for( int j = 0; j <= i; ++j )
129  {
130  ret[ i ][ j ] = FieldType( 0 );
131  for( int k = 0; k < m; ++k )
132  ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
133  ret[ j ][ i ] = ret[ i ][ j ];
134  }
135 
136  ret[ i ][ i ] = FieldType( 0 );
137  for( int k = 0; k < m; ++k )
138  ret[ i ][ i ] += A[ k ][ i ] * A[ k ][ i ];
139  }
140  }
141 
142  template< int m, int n >
143  static void
144  AAT_L ( const typename Traits :: template Matrix< m, n > :: type &A,
145  typename Traits :: template Matrix< m, m > :: type &ret )
146  {
147  /*
148  if (m==2) {
149  ret[0][0] = A[0]*A[0];
150  ret[1][1] = A[1]*A[1];
151  ret[1][0] = A[0]*A[1];
152  }
153  else
154  */
155  for( int i = 0; i < m; ++i )
156  {
157  for( int j = 0; j <= i; ++j )
158  {
159  FieldType &retij = ret[ i ][ j ];
160  retij = A[ i ][ 0 ] * A[ j ][ 0 ];
161  for( int k = 1; k < n; ++k )
162  retij += A[ i ][ k ] * A[ j ][ k ];
163  }
164  }
165  }
166 
167  template< int m, int n >
168  static void
169  AAT ( const typename Traits :: template Matrix< m, n > :: type &A,
170  typename Traits :: template Matrix< m, m > :: type &ret )
171  {
172  for( int i = 0; i < m; ++i )
173  {
174  for( int j = 0; j < i; ++j )
175  {
176  ret[ i ][ j ] = FieldType( 0 );
177  for( int k = 0; k < n; ++k )
178  ret[ i ][ j ] += A[ i ][ k ] * A[ j ][ k ];
179  ret[ j ][ i ] = ret[ i ][ j ];
180  }
181  ret[ i ][ i ] = FieldType( 0 );
182  for( int k = 0; k < n; ++k )
183  ret[ i ][ i ] += A[ i ][ k ] * A[ i ][ k ];
184  }
185  }
186 
187  template< int n >
188  static void
189  Lx ( const typename Traits :: template Matrix< n, n > :: type &L,
190  const typename Traits :: template Vector< n > :: type &x,
191  typename Traits :: template Vector< n > :: type &ret )
192  {
193  for( int i = 0; i < n; ++i )
194  {
195  ret[ i ] = FieldType( 0 );
196  for( int j = 0; j <= i; ++j )
197  ret[ i ] += L[ i ][ j ] * x[ j ];
198  }
199  }
200 
201  template< int n >
202  static void
203  LTx ( const typename Traits :: template Matrix< n, n > :: type &L,
204  const typename Traits :: template Vector< n > :: type &x,
205  typename Traits :: template Vector< n > :: type &ret )
206  {
207  for( int i = 0; i < n; ++i )
208  {
209  ret[ i ] = FieldType( 0 );
210  for( int j = i; j < n; ++j )
211  ret[ i ] += L[ j ][ i ] * x[ j ];
212  }
213  }
214 
215  template< int n >
216  static void
217  LTL ( const typename Traits :: template Matrix< n, n > :: type &L,
218  typename Traits :: template Matrix< n, n > :: type &ret )
219  {
220  for( int i = 0; i < n; ++i )
221  {
222  for( int j = 0; j < i; ++j )
223  {
224  ret[ i ][ j ] = FieldType( 0 );
225  for( int k = i; k < n; ++k )
226  ret[ i ][ j ] += L[ k ][ i ] * L[ k ][ j ];
227  ret[ j ][ i ] = ret[ i ][ j ];
228  }
229  ret[ i ][ i ] = FieldType( 0 );
230  for( int k = i; k < n; ++k )
231  ret[ i ][ i ] += L[ k ][ i ] * L[ k ][ i ];
232  }
233  }
234 
235  template< int n >
236  static void
237  LLT ( const typename Traits :: template Matrix< n, n > :: type &L,
238  typename Traits :: template Matrix< n, n > :: type &ret )
239  {
240  for( int i = 0; i < n; ++i )
241  {
242  for( int j = 0; j < i; ++j )
243  {
244  ret[ i ][ j ] = FieldType( 0 );
245  for( int k = 0; k <= j; ++k )
246  ret[ i ][ j ] += L[ i ][ k ] * L[ j ][ k ];
247  ret[ j ][ i ] = ret[ i ][ j ];
248  }
249  ret[ i ][ i ] = FieldType( 0 );
250  for( int k = 0; k <= i; ++k )
251  ret[ i ][ i ] += L[ i ][ k ] * L[ i ][ k ];
252  }
253  }
254 
255  template< int n >
256  static void
257  cholesky_L ( const typename Traits :: template Matrix< n, n > :: type &A,
258  typename Traits :: template Matrix< n, n > :: type &ret )
259  {
260  for( int i = 0; i < n; ++i )
261  {
262  FieldType &rii = ret[ i ][ i ];
263 
264  FieldType x = A[ i ][ i ];
265  for( int j = 0; j < i; ++j )
266  x -= ret[ i ][ j ] * ret[ i ][ j ];
267  assert( x > FieldType( 0 ) );
268  rii = sqrt( x );
269 
270  FieldType invrii = FieldType( 1 ) / rii;
271  for( int k = i+1; k < n; ++k )
272  {
273  FieldType x = A[ k ][ i ];
274  for( int j = 0; j < i; ++j )
275  x -= ret[ i ][ j ] * ret[ k ][ j ];
276  ret[ k ][ i ] = invrii * x;
277  }
278  }
279  }
280 
281  template< int n >
282  static FieldType
283  detL ( const typename Traits :: template Matrix< n, n > :: type &L )
284  {
285  FieldType det = FieldType( 1 );
286  for( int i = 0; i < n; ++i )
287  det *= L[ i ][ i ];
288  return det;
289  }
290 
291  template< int n >
292  static FieldType
293  invL ( typename Traits :: template Matrix< n, n > :: type &L )
294  {
295  FieldType det = FieldType( 1 );
296  for( int i = 0; i < n; ++i )
297  {
298  FieldType &lii = L[ i ][ i ];
299  det *= lii;
300  lii = FieldType( 1 ) / lii;
301  for( int j = 0; j < i; ++j )
302  {
303  FieldType &lij = L[ i ][ j ];
304  FieldType x = lij * L[ j ][ j ];
305  for( int k = j+1; k < i; ++k )
306  x += L[ i ][ k ] * L[ k ][ j ];
307  lij = (-lii) * x;
308  }
309  }
310  return det;
311  }
312 
313  // calculates x := L^{-1} x
314  template< int n >
315  static void
316  invLx ( typename Traits :: template Matrix< n, n > :: type &L,
317  typename Traits :: template Vector< n > :: type &x )
318  {
319  for( int i = 0; i < n; ++i )
320  {
321  for( int j = 0; j < i; ++j )
322  x[ i ] -= L[ i ][ j ] * x[ j ];
323  x[ i ] /= L[ i ][ i ];
324  }
325  }
326 
327  // calculates x := L^{-T} x
328  template< int n >
329  static void
330  invLTx ( typename Traits :: template Matrix< n, n > :: type &L,
331  typename Traits :: template Vector< n > :: type &x )
332  {
333  for( int i = n; i > 0; --i )
334  {
335  for( int j = i; j < n; ++j )
336  x[ i-1 ] -= L[ j ][ i-1 ] * x[ j ];
337  x[ i-1 ] /= L[ i-1 ][ i-1 ];
338  }
339  }
340 
341  template< int n >
342  static FieldType
343  spdDetA ( const typename Traits :: template Matrix< n, n > :: type &A )
344  {
345  // return A[0][0]*A[1][1]-A[1][0]*A[1][0];
346  typename Traits :: template Matrix< n, n > :: type L;
347  cholesky_L< n >( A, L );
348  return detL< n >( L );
349  }
350 
351  template< int n >
352  static FieldType
353  spdInvA ( typename Traits :: template Matrix< n, n > :: type &A )
354  {
355  typename Traits :: template Matrix< n, n > :: type L;
356  cholesky_L< n >( A, L );
357  const FieldType det = invL< n >( L );
358  LTL< n >( L, A );
359  return det;
360  }
361 
362  // calculate x := A^{-1} x
363  template< int n >
364  static void
365  spdInvAx ( typename Traits :: template Matrix< n, n > :: type &A,
366  typename Traits :: template Vector< n > :: type &x )
367  {
368  typename Traits :: template Matrix< n, n > :: type L;
369  cholesky_L< n >( A, L );
370  invLx< n >( L, x );
371  invLTx< n >( L, x );
372  }
373 
374  template< int m, int n >
375  static FieldType
376  detATA ( const typename Traits :: template Matrix< m, n > :: type &A )
377  {
378  if( m >= n )
379  {
380  typename Traits :: template Matrix< n, n > :: type ata;
381  ATA_L< m, n >( A, ata );
382  return spdDetA< n >( ata );
383  }
384  else
385  return FieldType( 0 );
386  }
387 
393  template< int m, int n >
394  static FieldType
395  sqrtDetAAT ( const typename Traits::template Matrix< m, n >::type &A )
396  {
397  // These special cases are here not only for speed reasons:
398  // The general implementation aborts if the matrix is almost singular,
399  // and the special implementation provide a stable way to handle that case.
400  if( (n == 2) && (m == 2) )
401  {
402  // Special implementation for 2x2 matrices: faster and more stable
403  return abs( A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ] );
404  }
405  else if( (n == 3) && (m == 3) )
406  {
407  // Special implementation for 3x3 matrices
408  const FieldType v0 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 1 ][ 1 ] * A[ 0 ][ 2 ];
409  const FieldType v1 = A[ 0 ][ 2 ] * A[ 1 ][ 0 ] - A[ 1 ][ 2 ] * A[ 0 ][ 0 ];
410  const FieldType v2 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 1 ][ 0 ] * A[ 0 ][ 1 ];
411  return abs( v0 * A[ 2 ][ 0 ] + v1 * A[ 2 ][ 1 ] + v2 * A[ 2 ][ 2 ] );
412  }
413  else if ( (n == 3) && (m == 2) )
414  {
415  // Special implementation for 2x3 matrices
416  const FieldType v0 = A[ 0 ][ 0 ] * A[ 1 ][ 1 ] - A[ 0 ][ 1 ] * A[ 1 ][ 0 ];
417  const FieldType v1 = A[ 0 ][ 0 ] * A[ 1 ][ 2 ] - A[ 1 ][ 0 ] * A[ 0 ][ 2 ];
418  const FieldType v2 = A[ 0 ][ 1 ] * A[ 1 ][ 2 ] - A[ 0 ][ 2 ] * A[ 1 ][ 1 ];
419  return sqrt( v0*v0 + v1*v1 + v2*v2);
420  }
421  else if( n >= m )
422  {
423  // General case
424  typename Traits::template Matrix< m, m >::type aat;
425  AAT_L< m, n >( A, aat );
426  return spdDetA< m >( aat );
427  }
428  else
429  return FieldType( 0 );
430  }
431 
432  // A^{-1}_L = (A^T A)^{-1} A^T
433  // => A^{-1}_L A = I
434  template< int m, int n >
435  static FieldType
436  leftInvA ( const typename Traits :: template Matrix< m, n > :: type &A,
437  typename Traits :: template Matrix< n, m > :: type &ret )
438  {
439  dune_static_assert( (m >= n), "Matrix has no left inverse." );
440  typename Traits :: template Matrix< n, n > :: type ata;
441  ATA_L< m, n >( A, ata );
442  const FieldType det = spdInvA< n >( ata );
443  ATBT< n, n, m >( ata, A, ret );
444  return det;
445  }
446 
447  template< int m, int n >
448  static void
449  leftInvAx ( const typename Traits :: template Matrix< m, n > :: type &A,
450  const typename Traits :: template Vector< m > :: type &x,
451  typename Traits :: template Vector< n > :: type &y )
452  {
453  dune_static_assert( (m >= n), "Matrix has no left inverse." );
454  typename Traits :: template Matrix< n, n > :: type ata;
455  ATx< m, n >( A, x, y );
456  ATA_L< m, n >( A, ata );
457  spdInvAx< n >( ata, y );
458  }
459 
461  template< int m, int n >
462  static FieldType
463  rightInvA ( const typename Traits :: template Matrix< m, n > :: type &A,
464  typename Traits :: template Matrix< n, m > :: type &ret )
465  {
466  dune_static_assert( (n >= m), "Matrix has no right inverse." );
467  if( (n == 2) && (m == 2) )
468  {
469  const FieldType det = (A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ]);
470  const FieldType detInv = FieldType( 1 ) / det;
471  ret[ 0 ][ 0 ] = A[ 1 ][ 1 ] * detInv;
472  ret[ 1 ][ 1 ] = A[ 0 ][ 0 ] * detInv;
473  ret[ 1 ][ 0 ] = -A[ 1 ][ 0 ] * detInv;
474  ret[ 0 ][ 1 ] = -A[ 0 ][ 1 ] * detInv;
475  return abs( det );
476  }
477  else
478  {
479  typename Traits :: template Matrix< m , m > :: type aat;
480  AAT_L< m, n >( A, aat );
481  const FieldType det = spdInvA< m >( aat );
482  ATBT< m, n, m >( A , aat , ret );
483  return det;
484  }
485  }
486 
487  template< int m, int n >
488  static void
489  xTRightInvA ( const typename Traits :: template Matrix< m, n > :: type &A,
490  const typename Traits :: template Vector< n > :: type &x,
491  typename Traits :: template Vector< m > :: type &y )
492  {
493  dune_static_assert( (n >= m), "Matrix has no right inverse." );
494  typename Traits :: template Matrix< m, m > :: type aat;
495  Ax< m, n >( A, x, y );
496  AAT_L< m, n >( A, aat );
497  spdInvAx< m >( aat, y );
498  }
499  };
500 
501  } // namespace GenericGeometry
502 
503 } // namespace Dune
504 
505 #endif // #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_MATRIXHELPER_HH
static void Ax(const typename Traits::template Matrix< m, n >::type &A, const typename Traits::template Vector< n >::type &x, typename Traits::template Vector< m >::type &ret)
Definition: matrixhelper.hh:45
static FieldType abs(const FieldType &x)
Definition: matrixhelper.hh:37
static FieldType detL(const typename Traits::template Matrix< n, n >::type &L)
Definition: matrixhelper.hh:283
Definition: matrixhelper.hh:33
static Field abs(const Field &x)
Definition: matrixhelper.hh:24
static FieldType spdInvA(typename Traits::template Matrix< n, n >::type &A)
Definition: matrixhelper.hh:353
static void ATx(const typename Traits::template Matrix< m, n >::type &A, const typename Traits::template Vector< m >::type &x, typename Traits::template Vector< n >::type &ret)
Definition: matrixhelper.hh:59
static FieldType leftInvA(const typename Traits::template Matrix< m, n >::type &A, typename Traits::template Matrix< n, m >::type &ret)
Definition: matrixhelper.hh:436
static void ATA(const typename Traits::template Matrix< m, n >::type &A, typename Traits::template Matrix< n, n >::type &ret)
Definition: matrixhelper.hh:123
static FieldType sqrtDetAAT(const typename Traits::template Matrix< m, n >::type &A)
Compute the square root of the determinant of A times A transposed.
Definition: matrixhelper.hh:395
static void AB(const typename Traits::template Matrix< m, n >::type &A, const typename Traits::template Matrix< n, p >::type &B, typename Traits::template Matrix< m, p >::type &ret)
Definition: matrixhelper.hh:73
static void AAT(const typename Traits::template Matrix< m, n >::type &A, typename Traits::template Matrix< m, m >::type &ret)
Definition: matrixhelper.hh:169
static FieldType detATA(const typename Traits::template Matrix< m, n >::type &A)
Definition: matrixhelper.hh:376
static void leftInvAx(const typename Traits::template Matrix< m, n >::type &A, const typename Traits::template Vector< m >::type &x, typename Traits::template Vector< n >::type &y)
Definition: matrixhelper.hh:449
static void ATBT(const typename Traits::template Matrix< m, n >::type &A, const typename Traits::template Matrix< p, m >::type &B, typename Traits::template Matrix< n, p >::type &ret)
Definition: matrixhelper.hh:90
Definition: matrixhelper.hh:22
static void AAT_L(const typename Traits::template Matrix< m, n >::type &A, typename Traits::template Matrix< m, m >::type &ret)
Definition: matrixhelper.hh:144
static void xTRightInvA(const typename Traits::template Matrix< m, n >::type &A, const typename Traits::template Vector< n >::type &x, typename Traits::template Vector< m >::type &y)
Definition: matrixhelper.hh:489
static void LLT(const typename Traits::template Matrix< n, n >::type &L, typename Traits::template Matrix< n, n >::type &ret)
Definition: matrixhelper.hh:237
static FieldType spdDetA(const typename Traits::template Matrix< n, n >::type &A)
Definition: matrixhelper.hh:343
static void ATA_L(const typename Traits::template Matrix< m, n >::type &A, typename Traits::template Matrix< n, n >::type &ret)
Definition: matrixhelper.hh:107
static void invLx(typename Traits::template Matrix< n, n >::type &L, typename Traits::template Vector< n >::type &x)
Definition: matrixhelper.hh:316
static void LTx(const typename Traits::template Matrix< n, n >::type &L, const typename Traits::template Vector< n >::type &x, typename Traits::template Vector< n >::type &ret)
Definition: matrixhelper.hh:203
static FieldType rightInvA(const typename Traits::template Matrix< m, n >::type &A, typename Traits::template Matrix< n, m >::type &ret)
Compute right pseudo-inverse of matrix A.
Definition: matrixhelper.hh:463
static void invLTx(typename Traits::template Matrix< n, n >::type &L, typename Traits::template Vector< n >::type &x)
Definition: matrixhelper.hh:330
static void cholesky_L(const typename Traits::template Matrix< n, n >::type &A, typename Traits::template Matrix< n, n >::type &ret)
Definition: matrixhelper.hh:257
static void LTL(const typename Traits::template Matrix< n, n >::type &L, typename Traits::template Matrix< n, n >::type &ret)
Definition: matrixhelper.hh:217
static void Lx(const typename Traits::template Matrix< n, n >::type &L, const typename Traits::template Vector< n >::type &x, typename Traits::template Vector< n >::type &ret)
Definition: matrixhelper.hh:189
Traits::ctype FieldType
Definition: matrixhelper.hh:35
static FieldType invL(typename Traits::template Matrix< n, n >::type &L)
Definition: matrixhelper.hh:293
static void spdInvAx(typename Traits::template Matrix< n, n >::type &A, typename Traits::template Vector< n >::type &x)
Definition: matrixhelper.hh:365