3 #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_MATRIXHELPER_HH
4 #define DUNE_GEOMETRY_GENERICGEOMETRY_MATRIXHELPER_HH
8 #include <dune/common/fvector.hh>
9 #include <dune/common/fmatrix.hh>
10 #include <dune/common/static_assert.hh>
15 namespace GenericGeometry
21 template<
class Field >
24 static Field
abs (
const Field &x ) {
return std::abs( x ); }
32 template<
class Traits >
43 template<
int m,
int n >
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 )
49 for(
int i = 0; i < m; ++i )
52 for(
int j = 0; j < n; ++j )
53 ret[ i ] += A[ i ][ j ] * x[ j ];
57 template<
int m,
int n >
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 )
63 for(
int i = 0; i < n; ++i )
66 for(
int j = 0; j < m; ++j )
67 ret[ i ] += A[ j ][ i ] * x[ j ];
71 template<
int m,
int n,
int p >
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 )
77 for(
int i = 0; i < m; ++i )
79 for(
int j = 0; j < p; ++j )
82 for(
int k = 0; k < n; ++k )
83 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
88 template<
int m,
int n,
int p >
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 )
94 for(
int i = 0; i < n; ++i )
96 for(
int j = 0; j < p; ++j )
99 for(
int k = 0; k < m; ++k )
100 ret[ i ][ j ] += A[ k ][ i ] * B[ j ][ k ];
105 template<
int m,
int n >
107 ATA_L (
const typename Traits :: template Matrix< m, n > :: type &A,
108 typename Traits :: template Matrix< n, n > :: type &ret )
110 for(
int i = 0; i < n; ++i )
112 for(
int j = 0; j <= i; ++j )
115 for(
int k = 0; k < m; ++k )
116 ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
121 template<
int m,
int n >
123 ATA (
const typename Traits :: template Matrix< m, n > :: type &A,
124 typename Traits :: template Matrix< n, n > :: type &ret )
126 for(
int i = 0; i < n; ++i )
128 for(
int j = 0; j <= i; ++j )
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 ];
137 for(
int k = 0; k < m; ++k )
138 ret[ i ][ i ] += A[ k ][ i ] * A[ k ][ i ];
142 template<
int m,
int n >
144 AAT_L (
const typename Traits :: template Matrix< m, n > :: type &A,
145 typename Traits :: template Matrix< m, m > :: type &ret )
155 for(
int i = 0; i < m; ++i )
157 for(
int j = 0; j <= 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 ];
167 template<
int m,
int n >
169 AAT (
const typename Traits :: template Matrix< m, n > :: type &A,
170 typename Traits :: template Matrix< m, m > :: type &ret )
172 for(
int i = 0; i < m; ++i )
174 for(
int j = 0; j < i; ++j )
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 ];
182 for(
int k = 0; k < n; ++k )
183 ret[ i ][ i ] += A[ i ][ k ] * A[ i ][ k ];
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 )
193 for(
int i = 0; i < n; ++i )
196 for(
int j = 0; j <= i; ++j )
197 ret[ i ] += L[ i ][ j ] * x[ j ];
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 )
207 for(
int i = 0; i < n; ++i )
210 for(
int j = i; j < n; ++j )
211 ret[ i ] += L[ j ][ i ] * x[ j ];
217 LTL (
const typename Traits :: template Matrix< n, n > :: type &L,
218 typename Traits :: template Matrix< n, n > :: type &ret )
220 for(
int i = 0; i < n; ++i )
222 for(
int j = 0; j < i; ++j )
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 ];
230 for(
int k = i; k < n; ++k )
231 ret[ i ][ i ] += L[ k ][ i ] * L[ k ][ i ];
237 LLT (
const typename Traits :: template Matrix< n, n > :: type &L,
238 typename Traits :: template Matrix< n, n > :: type &ret )
240 for(
int i = 0; i < n; ++i )
242 for(
int j = 0; j < i; ++j )
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 ];
250 for(
int k = 0; k <= i; ++k )
251 ret[ i ][ i ] += L[ i ][ k ] * L[ i ][ k ];
257 cholesky_L (
const typename Traits :: template Matrix< n, n > :: type &A,
258 typename Traits :: template Matrix< n, n > :: type &ret )
260 for(
int i = 0; i < n; ++i )
265 for(
int j = 0; j < i; ++j )
266 x -= ret[ i ][ j ] * ret[ i ][ j ];
271 for(
int k = i+1; k < n; ++k )
274 for(
int j = 0; j < i; ++j )
275 x -= ret[ i ][ j ] * ret[ k ][ j ];
276 ret[ k ][ i ] = invrii * x;
283 detL (
const typename Traits :: template Matrix< n, n > :: type &L )
286 for(
int i = 0; i < n; ++i )
293 invL (
typename Traits :: template Matrix< n, n > :: type &L )
296 for(
int i = 0; i < n; ++i )
301 for(
int j = 0; j < i; ++j )
305 for(
int k = j+1; k < i; ++k )
306 x += L[ i ][ k ] * L[ k ][ j ];
316 invLx (
typename Traits :: template Matrix< n, n > :: type &L,
317 typename Traits :: template Vector< n > :: type &x )
319 for(
int i = 0; i < n; ++i )
321 for(
int j = 0; j < i; ++j )
322 x[ i ] -= L[ i ][ j ] * x[ j ];
323 x[ i ] /= L[ i ][ i ];
330 invLTx (
typename Traits :: template Matrix< n, n > :: type &L,
331 typename Traits :: template Vector< n > :: type &x )
333 for(
int i = n; i > 0; --i )
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 ];
343 spdDetA (
const typename Traits :: template Matrix< n, n > :: type &A )
346 typename Traits :: template Matrix< n, n > :: type L;
347 cholesky_L< n >( A, L );
348 return detL< n >( L );
353 spdInvA (
typename Traits :: template Matrix< n, n > :: type &A )
355 typename Traits :: template Matrix< n, n > :: type L;
356 cholesky_L< n >( A, L );
365 spdInvAx (
typename Traits :: template Matrix< n, n > :: type &A,
366 typename Traits :: template Vector< n > :: type &x )
368 typename Traits :: template Matrix< n, n > :: type L;
369 cholesky_L< n >( A, L );
374 template<
int m,
int n >
376 detATA (
const typename Traits :: template Matrix< m, n > :: type &A )
380 typename Traits :: template Matrix< n, n > :: type ata;
381 ATA_L< m, n >( A, ata );
382 return spdDetA< n >( ata );
393 template<
int m,
int n >
395 sqrtDetAAT (
const typename Traits::template Matrix< m, n >::type &A )
400 if( (n == 2) && (m == 2) )
403 return abs( A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ] );
405 else if( (n == 3) && (m == 3) )
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 ] );
413 else if ( (n == 3) && (m == 2) )
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);
424 typename Traits::template Matrix< m, m >::type aat;
425 AAT_L< m, n >( A, aat );
426 return spdDetA< m >( aat );
434 template<
int m,
int n >
436 leftInvA (
const typename Traits :: template Matrix< m, n > :: type &A,
437 typename Traits :: template Matrix< n, m > :: type &ret )
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 );
447 template<
int m,
int n >
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 )
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 );
461 template<
int m,
int n >
463 rightInvA (
const typename Traits :: template Matrix< m, n > :: type &A,
464 typename Traits :: template Matrix< n, m > :: type &ret )
466 dune_static_assert( (n >= m),
"Matrix has no right inverse." );
467 if( (n == 2) && (m == 2) )
469 const FieldType det = (A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ]);
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;
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 );
487 template<
int m,
int n >
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 )
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 );
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