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