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>
14 namespace GenericGeometry
20 template<
class Field >
23 static Field
abs (
const Field &x ) {
34 template<
class Traits >
45 template<
int m,
int n >
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 )
51 for(
int i = 0; i < m; ++i )
54 for(
int j = 0; j < n; ++j )
55 ret[ i ] += A[ i ][ j ] * x[ j ];
59 template<
int m,
int n >
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 )
65 for(
int i = 0; i < n; ++i )
68 for(
int j = 0; j < m; ++j )
69 ret[ i ] += A[ j ][ i ] * x[ j ];
73 template<
int m,
int n,
int p >
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 )
79 for(
int i = 0; i < m; ++i )
81 for(
int j = 0; j < p; ++j )
84 for(
int k = 0; k < n; ++k )
85 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
90 template<
int m,
int n,
int p >
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 )
96 for(
int i = 0; i < n; ++i )
98 for(
int j = 0; j < p; ++j )
101 for(
int k = 0; k < m; ++k )
102 ret[ i ][ j ] += A[ k ][ i ] * B[ j ][ k ];
107 template<
int m,
int n >
109 ATA_L (
const typename Traits :: template Matrix< m, n > :: type &A,
110 typename Traits :: template Matrix< n, n > :: type &ret )
112 for(
int i = 0; i < n; ++i )
114 for(
int j = 0; j <= i; ++j )
117 for(
int k = 0; k < m; ++k )
118 ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
123 template<
int m,
int n >
125 ATA (
const typename Traits :: template Matrix< m, n > :: type &A,
126 typename Traits :: template Matrix< n, n > :: type &ret )
128 for(
int i = 0; i < n; ++i )
130 for(
int j = 0; j <= i; ++j )
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 ];
139 for(
int k = 0; k < m; ++k )
140 ret[ i ][ i ] += A[ k ][ i ] * A[ k ][ i ];
144 template<
int m,
int n >
146 AAT_L (
const typename Traits :: template Matrix< m, n > :: type &A,
147 typename Traits :: template Matrix< m, m > :: type &ret )
157 for(
int i = 0; i < m; ++i )
159 for(
int j = 0; j <= 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 ];
169 template<
int m,
int n >
171 AAT (
const typename Traits :: template Matrix< m, n > :: type &A,
172 typename Traits :: template Matrix< m, m > :: type &ret )
174 for(
int i = 0; i < m; ++i )
176 for(
int j = 0; j < i; ++j )
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 ];
184 for(
int k = 0; k < n; ++k )
185 ret[ i ][ i ] += A[ i ][ k ] * A[ i ][ k ];
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 )
195 for(
int i = 0; i < n; ++i )
198 for(
int j = 0; j <= i; ++j )
199 ret[ i ] += L[ i ][ j ] * x[ j ];
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 )
209 for(
int i = 0; i < n; ++i )
212 for(
int j = i; j < n; ++j )
213 ret[ i ] += L[ j ][ i ] * x[ j ];
219 LTL (
const typename Traits :: template Matrix< n, n > :: type &L,
220 typename Traits :: template Matrix< n, n > :: type &ret )
222 for(
int i = 0; i < n; ++i )
224 for(
int j = 0; j < i; ++j )
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 ];
232 for(
int k = i; k < n; ++k )
233 ret[ i ][ i ] += L[ k ][ i ] * L[ k ][ i ];
239 LLT (
const typename Traits :: template Matrix< n, n > :: type &L,
240 typename Traits :: template Matrix< n, n > :: type &ret )
242 for(
int i = 0; i < n; ++i )
244 for(
int j = 0; j < i; ++j )
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 ];
252 for(
int k = 0; k <= i; ++k )
253 ret[ i ][ i ] += L[ i ][ k ] * L[ i ][ k ];
259 cholesky_L (
const typename Traits :: template Matrix< n, n > :: type &A,
260 typename Traits :: template Matrix< n, n > :: type &ret )
262 for(
int i = 0; i < n; ++i )
267 for(
int j = 0; j < i; ++j )
268 xDiag -= ret[ i ][ j ] * ret[ i ][ j ];
273 for(
int k = i+1; k < n; ++k )
276 for(
int j = 0; j < i; ++j )
277 x -= ret[ i ][ j ] * ret[ k ][ j ];
278 ret[ k ][ i ] = invrii * x;
285 detL (
const typename Traits :: template Matrix< n, n > :: type &L )
288 for(
int i = 0; i < n; ++i )
295 invL (
typename Traits :: template Matrix< n, n > :: type &L )
298 for(
int i = 0; i < n; ++i )
303 for(
int j = 0; j < i; ++j )
307 for(
int k = j+1; k < i; ++k )
308 x += L[ i ][ k ] * L[ k ][ j ];
318 invLx (
typename Traits :: template Matrix< n, n > :: type &L,
319 typename Traits :: template Vector< n > :: type &x )
321 for(
int i = 0; i < n; ++i )
323 for(
int j = 0; j < i; ++j )
324 x[ i ] -= L[ i ][ j ] * x[ j ];
325 x[ i ] /= L[ i ][ i ];
332 invLTx (
typename Traits :: template Matrix< n, n > :: type &L,
333 typename Traits :: template Vector< n > :: type &x )
335 for(
int i = n; i > 0; --i )
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 ];
345 spdDetA (
const typename Traits :: template Matrix< n, n > :: type &A )
348 typename Traits :: template Matrix< n, n > :: type L;
349 cholesky_L< n >( A, L );
350 return detL< n >( L );
355 spdInvA (
typename Traits :: template Matrix< n, n > :: type &A )
357 typename Traits :: template Matrix< n, n > :: type L;
358 cholesky_L< n >( A, L );
367 spdInvAx (
typename Traits :: template Matrix< n, n > :: type &A,
368 typename Traits :: template Vector< n > :: type &x )
370 typename Traits :: template Matrix< n, n > :: type L;
371 cholesky_L< n >( A, L );
376 template<
int m,
int n >
378 detATA (
const typename Traits :: template Matrix< m, n > :: type &A )
382 typename Traits :: template Matrix< n, n > :: type ata;
383 ATA_L< m, n >( A, ata );
384 return spdDetA< n >( ata );
395 template<
int m,
int n >
397 sqrtDetAAT (
const typename Traits::template Matrix< m, n >::type &A )
402 if( (n == 2) && (m == 2) )
405 return abs( A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ] );
407 else if( (n == 3) && (m == 3) )
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 ] );
415 else if ( (n == 3) && (m == 2) )
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);
426 typename Traits::template Matrix< m, m >::type aat;
427 AAT_L< m, n >( A, aat );
428 return spdDetA< m >( aat );
436 template<
int m,
int n >
438 leftInvA (
const typename Traits :: template Matrix< m, n > :: type &A,
439 typename Traits :: template Matrix< n, m > :: type &ret )
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 );
449 template<
int m,
int n >
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 )
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 );
463 template<
int m,
int n >
465 rightInvA (
const typename Traits :: template Matrix< m, n > :: type &A,
466 typename Traits :: template Matrix< n, m > :: type &ret )
468 static_assert((n >= m),
"Matrix has no right inverse.");
469 if( (n == 2) && (m == 2) )
471 const FieldType det = (A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ]);
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;
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 );
489 template<
int m,
int n >
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 )
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 );
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