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