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 ) {
return std::abs( x ); }
31 template<
class Traits >
36 static FieldType
abs (
const FieldType &x )
42 template<
int m,
int n >
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 )
48 for(
int i = 0; i < m; ++i )
51 for(
int j = 0; j < n; ++j )
52 ret[ i ] += A[ i ][ j ] * x[ j ];
56 template<
int m,
int n >
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 )
62 for(
int i = 0; i < n; ++i )
65 for(
int j = 0; j < m; ++j )
66 ret[ i ] += A[ j ][ i ] * x[ j ];
70 template<
int m,
int n,
int p >
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 )
76 for(
int i = 0; i < m; ++i )
78 for(
int j = 0; j < p; ++j )
81 for(
int k = 0; k < n; ++k )
82 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
87 template<
int m,
int n,
int p >
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 )
93 for(
int i = 0; i < n; ++i )
95 for(
int j = 0; j < p; ++j )
98 for(
int k = 0; k < m; ++k )
99 ret[ i ][ j ] += A[ k ][ i ] * B[ j ][ k ];
104 template<
int m,
int n >
106 ATA_L (
const typename Traits :: template Matrix< m, n > :: type &A,
107 typename Traits :: template Matrix< n, n > :: type &ret )
109 for(
int i = 0; i < n; ++i )
111 for(
int j = 0; j <= i; ++j )
114 for(
int k = 0; k < m; ++k )
115 ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
120 template<
int m,
int n >
122 ATA (
const typename Traits :: template Matrix< m, n > :: type &A,
123 typename Traits :: template Matrix< n, n > :: type &ret )
125 for(
int i = 0; i < n; ++i )
127 for(
int j = 0; j <= i; ++j )
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 ];
136 for(
int k = 0; k < m; ++k )
137 ret[ i ][ i ] += A[ k ][ i ] * A[ k ][ i ];
141 template<
int m,
int n >
143 AAT_L (
const typename Traits :: template Matrix< m, n > :: type &A,
144 typename Traits :: template Matrix< m, m > :: type &ret )
154 for(
int i = 0; i < m; ++i )
156 for(
int j = 0; j <= i; ++j )
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 ];
166 template<
int m,
int n >
168 AAT (
const typename Traits :: template Matrix< m, n > :: type &A,
169 typename Traits :: template Matrix< m, m > :: type &ret )
171 for(
int i = 0; i < m; ++i )
173 for(
int j = 0; j < i; ++j )
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 ];
181 for(
int k = 0; k < n; ++k )
182 ret[ i ][ i ] += A[ i ][ k ] * A[ i ][ k ];
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 )
192 for(
int i = 0; i < n; ++i )
195 for(
int j = 0; j <= i; ++j )
196 ret[ i ] += L[ i ][ j ] * x[ j ];
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 )
206 for(
int i = 0; i < n; ++i )
209 for(
int j = i; j < n; ++j )
210 ret[ i ] += L[ j ][ i ] * x[ j ];
216 LTL (
const typename Traits :: template Matrix< n, n > :: type &L,
217 typename Traits :: template Matrix< n, n > :: type &ret )
219 for(
int i = 0; i < n; ++i )
221 for(
int j = 0; j < i; ++j )
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 ];
229 for(
int k = i; k < n; ++k )
230 ret[ i ][ i ] += L[ k ][ i ] * L[ k ][ i ];
236 LLT (
const typename Traits :: template Matrix< n, n > :: type &L,
237 typename Traits :: template Matrix< n, n > :: type &ret )
239 for(
int i = 0; i < n; ++i )
241 for(
int j = 0; j < i; ++j )
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 ];
249 for(
int k = 0; k <= i; ++k )
250 ret[ i ][ i ] += L[ i ][ k ] * L[ i ][ k ];
256 cholesky_L (
const typename Traits :: template Matrix< n, n > :: type &A,
257 typename Traits :: template Matrix< n, n > :: type &ret )
259 for(
int i = 0; i < n; ++i )
261 FieldType &rii = ret[ i ][ i ];
263 FieldType x = A[ i ][ i ];
264 for(
int j = 0; j < i; ++j )
265 x -= ret[ i ][ j ] * ret[ i ][ j ];
270 for(
int k = i+1; k < n; ++k )
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;
282 detL (
const typename Traits :: template Matrix< n, n > :: type &L )
285 for(
int i = 0; i < n; ++i )
292 invL (
typename Traits :: template Matrix< n, n > :: type &L )
295 for(
int i = 0; i < n; ++i )
297 FieldType &lii = L[ i ][ i ];
300 for(
int j = 0; j < i; ++j )
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 ];
315 invLx (
typename Traits :: template Matrix< n, n > :: type &L,
316 typename Traits :: template Vector< n > :: type &x )
318 for(
int i = 0; i < n; ++i )
320 for(
int j = 0; j < i; ++j )
321 x[ i ] -= L[ i ][ j ] * x[ j ];
322 x[ i ] /= L[ i ][ i ];
329 invLTx (
typename Traits :: template Matrix< n, n > :: type &L,
330 typename Traits :: template Vector< n > :: type &x )
332 for(
int i = n; i > 0; --i )
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 ];
342 spdDetA (
const typename Traits :: template Matrix< n, n > :: type &A )
345 typename Traits :: template Matrix< n, n > :: type L;
346 cholesky_L< n >( A, L );
347 return detL< n >( L );
352 spdInvA (
typename Traits :: template Matrix< n, n > :: type &A )
354 typename Traits :: template Matrix< n, n > :: type L;
355 cholesky_L< n >( A, L );
356 const FieldType det = invL< n >( L );
364 spdInvAx (
typename Traits :: template Matrix< n, n > :: type &A,
365 typename Traits :: template Vector< n > :: type &x )
367 typename Traits :: template Matrix< n, n > :: type L;
368 cholesky_L< n >( A, L );
373 template<
int m,
int n >
375 detATA (
const typename Traits :: template Matrix< m, n > :: type &A )
379 typename Traits :: template Matrix< n, n > :: type ata;
380 ATA_L< m, n >( A, ata );
381 return spdDetA< n >( ata );
392 template<
int m,
int n >
394 sqrtDetAAT (
const typename Traits::template Matrix< m, n >::type &A )
399 if( (n == 2) && (m == 2) )
402 return abs( A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ] );
404 else if( (n == 3) && (m == 3) )
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 ] );
412 else if ( (n == 3) && (m == 2) )
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);
423 typename Traits::template Matrix< m, m >::type aat;
424 AAT_L< m, n >( A, aat );
425 return spdDetA< m >( aat );
433 template<
int m,
int n >
435 leftInvA (
const typename Traits :: template Matrix< m, n > :: type &A,
436 typename Traits :: template Matrix< n, m > :: type &ret )
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 );
446 template<
int m,
int n >
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 )
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 );
460 template<
int m,
int n >
462 rightInvA (
const typename Traits :: template Matrix< m, n > :: type &A,
463 typename Traits :: template Matrix< n, m > :: type &ret )
465 static_assert((n >= m),
"Matrix has no right inverse.");
466 if( (n == 2) && (m == 2) )
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;
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 );
486 template<
int m,
int n >
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 )
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 );
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