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
11namespace Dune
12{
13
14 namespace GenericGeometry
15 {
16
17 // FieldHelper
18 // -----------
19
20 template< class Field >
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 >
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