3#ifndef DUNE_GEOMETRY_GENERICGEOMETRY_MATRIXHELPER_HH
4#define DUNE_GEOMETRY_GENERICGEOMETRY_MATRIXHELPER_HH
14 namespace GenericGeometry
20 template<
class Field >
23 static Field abs (
const Field &x ) {
return std::abs( x ); }
31 template<
class Traits >
34 typedef typename Traits::ctype FieldType;
36 static FieldType abs (
const FieldType &x )
39 return FieldHelper< FieldType >::abs( 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 )
50 ret[ i ] = FieldType( 0 );
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 )
64 ret[ i ] = FieldType( 0 );
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 )
80 ret[ i ][ j ] = FieldType( 0 );
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 )
97 ret[ i ][ j ] = FieldType( 0 );
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 )
113 ret[ i ][ j ] = FieldType( 0 );
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 )
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 ];
135 ret[ i ][ i ] = FieldType( 0 );
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 )
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 ];
180 ret[ i ][ i ] = FieldType( 0 );
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 )
194 ret[ i ] = FieldType( 0 );
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 )
208 ret[ i ] = FieldType( 0 );
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 )
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 ];
228 ret[ i ][ i ] = FieldType( 0 );
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 )
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 ];
248 ret[ i ][ i ] = FieldType( 0 );
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 xDiag = A[ i ][ i ];
264 for(
int j = 0; j < i; ++j )
265 xDiag -= ret[ i ][ j ] * ret[ i ][ j ];
266 assert( xDiag > FieldType( 0 ) );
269 FieldType invrii = FieldType( 1 ) / rii;
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 )
284 FieldType det = FieldType( 1 );
285 for(
int i = 0; i < n; ++i )
292 invL (
typename Traits :: template Matrix< n, n > :: type &L )
294 FieldType det = FieldType( 1 );
295 for(
int i = 0; i < n; ++i )
297 FieldType &lii = L[ i ][ i ];
299 lii = FieldType( 1 ) / lii;
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 );
384 return FieldType( 0 );
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 );
428 return FieldType( 0 );
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 );
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune namespace.
Definition: alignment.hh:10