Dune Core Modules (2.3.1)

matrixhelper.hh
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#include <dune/common/static_assert.hh>
11
12namespace Dune
13{
14
15 namespace GenericGeometry
16 {
17
18 // FieldHelper
19 // -----------
20
21 template< class Field >
22 struct FieldHelper
23 {
24 static Field abs ( const Field &x ) { return std::abs( x ); }
25 };
26
27
28
29 // MatrixHelper
30 // ------------
31
32 template< class Traits >
33 struct MatrixHelper
34 {
35 typedef typename Traits::ctype FieldType;
36
37 static FieldType abs ( const FieldType &x )
38 {
39 //return std::abs( x );
40 return FieldHelper< FieldType >::abs( x );
41 }
42
43 template< int m, int n >
44 static void
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 )
48 {
49 for( int i = 0; i < m; ++i )
50 {
51 ret[ i ] = FieldType( 0 );
52 for( int j = 0; j < n; ++j )
53 ret[ i ] += A[ i ][ j ] * x[ j ];
54 }
55 }
56
57 template< int m, int n >
58 static void
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 )
62 {
63 for( int i = 0; i < n; ++i )
64 {
65 ret[ i ] = FieldType( 0 );
66 for( int j = 0; j < m; ++j )
67 ret[ i ] += A[ j ][ i ] * x[ j ];
68 }
69 }
70
71 template< int m, int n, int p >
72 static void
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 )
76 {
77 for( int i = 0; i < m; ++i )
78 {
79 for( int j = 0; j < p; ++j )
80 {
81 ret[ i ][ j ] = FieldType( 0 );
82 for( int k = 0; k < n; ++k )
83 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
84 }
85 }
86 }
87
88 template< int m, int n, int p >
89 static void
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 )
93 {
94 for( int i = 0; i < n; ++i )
95 {
96 for( int j = 0; j < p; ++j )
97 {
98 ret[ i ][ j ] = FieldType( 0 );
99 for( int k = 0; k < m; ++k )
100 ret[ i ][ j ] += A[ k ][ i ] * B[ j ][ k ];
101 }
102 }
103 }
104
105 template< int m, int n >
106 static void
107 ATA_L ( const typename Traits :: template Matrix< m, n > :: type &A,
108 typename Traits :: template Matrix< n, n > :: type &ret )
109 {
110 for( int i = 0; i < n; ++i )
111 {
112 for( int j = 0; j <= i; ++j )
113 {
114 ret[ i ][ j ] = FieldType( 0 );
115 for( int k = 0; k < m; ++k )
116 ret[ i ][ j ] += A[ k ][ i ] * A[ k ][ j ];
117 }
118 }
119 }
120
121 template< int m, int n >
122 static void
123 ATA ( const typename Traits :: template Matrix< m, n > :: type &A,
124 typename Traits :: template Matrix< n, n > :: type &ret )
125 {
126 for( int i = 0; i < n; ++i )
127 {
128 for( int j = 0; j <= i; ++j )
129 {
130 ret[ i ][ j ] = FieldType( 0 );
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 ];
134 }
135
136 ret[ i ][ i ] = FieldType( 0 );
137 for( int k = 0; k < m; ++k )
138 ret[ i ][ i ] += A[ k ][ i ] * A[ k ][ i ];
139 }
140 }
141
142 template< int m, int n >
143 static void
144 AAT_L ( const typename Traits :: template Matrix< m, n > :: type &A,
145 typename Traits :: template Matrix< m, m > :: type &ret )
146 {
147 /*
148 if (m==2) {
149 ret[0][0] = A[0]*A[0];
150 ret[1][1] = A[1]*A[1];
151 ret[1][0] = A[0]*A[1];
152 }
153 else
154 */
155 for( int i = 0; i < m; ++i )
156 {
157 for( int j = 0; j <= i; ++j )
158 {
159 FieldType &retij = ret[ 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 ];
163 }
164 }
165 }
166
167 template< int m, int n >
168 static void
169 AAT ( const typename Traits :: template Matrix< m, n > :: type &A,
170 typename Traits :: template Matrix< m, m > :: type &ret )
171 {
172 for( int i = 0; i < m; ++i )
173 {
174 for( int j = 0; j < i; ++j )
175 {
176 ret[ i ][ j ] = FieldType( 0 );
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 ];
180 }
181 ret[ i ][ i ] = FieldType( 0 );
182 for( int k = 0; k < n; ++k )
183 ret[ i ][ i ] += A[ i ][ k ] * A[ i ][ k ];
184 }
185 }
186
187 template< int n >
188 static void
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 )
192 {
193 for( int i = 0; i < n; ++i )
194 {
195 ret[ i ] = FieldType( 0 );
196 for( int j = 0; j <= i; ++j )
197 ret[ i ] += L[ i ][ j ] * x[ j ];
198 }
199 }
200
201 template< int n >
202 static void
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 )
206 {
207 for( int i = 0; i < n; ++i )
208 {
209 ret[ i ] = FieldType( 0 );
210 for( int j = i; j < n; ++j )
211 ret[ i ] += L[ j ][ i ] * x[ j ];
212 }
213 }
214
215 template< int n >
216 static void
217 LTL ( const typename Traits :: template Matrix< n, n > :: type &L,
218 typename Traits :: template Matrix< n, n > :: type &ret )
219 {
220 for( int i = 0; i < n; ++i )
221 {
222 for( int j = 0; j < i; ++j )
223 {
224 ret[ i ][ j ] = FieldType( 0 );
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 ];
228 }
229 ret[ i ][ i ] = FieldType( 0 );
230 for( int k = i; k < n; ++k )
231 ret[ i ][ i ] += L[ k ][ i ] * L[ k ][ i ];
232 }
233 }
234
235 template< int n >
236 static void
237 LLT ( const typename Traits :: template Matrix< n, n > :: type &L,
238 typename Traits :: template Matrix< n, n > :: type &ret )
239 {
240 for( int i = 0; i < n; ++i )
241 {
242 for( int j = 0; j < i; ++j )
243 {
244 ret[ i ][ j ] = FieldType( 0 );
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 ];
248 }
249 ret[ i ][ i ] = FieldType( 0 );
250 for( int k = 0; k <= i; ++k )
251 ret[ i ][ i ] += L[ i ][ k ] * L[ i ][ k ];
252 }
253 }
254
255 template< int n >
256 static void
257 cholesky_L ( const typename Traits :: template Matrix< n, n > :: type &A,
258 typename Traits :: template Matrix< n, n > :: type &ret )
259 {
260 for( int i = 0; i < n; ++i )
261 {
262 FieldType &rii = ret[ i ][ i ];
263
264 FieldType x = A[ i ][ i ];
265 for( int j = 0; j < i; ++j )
266 x -= ret[ i ][ j ] * ret[ i ][ j ];
267 assert( x > FieldType( 0 ) );
268 rii = sqrt( x );
269
270 FieldType invrii = FieldType( 1 ) / rii;
271 for( int k = i+1; k < n; ++k )
272 {
273 FieldType x = A[ k ][ i ];
274 for( int j = 0; j < i; ++j )
275 x -= ret[ i ][ j ] * ret[ k ][ j ];
276 ret[ k ][ i ] = invrii * x;
277 }
278 }
279 }
280
281 template< int n >
282 static FieldType
283 detL ( const typename Traits :: template Matrix< n, n > :: type &L )
284 {
285 FieldType det = FieldType( 1 );
286 for( int i = 0; i < n; ++i )
287 det *= L[ i ][ i ];
288 return det;
289 }
290
291 template< int n >
292 static FieldType
293 invL ( typename Traits :: template Matrix< n, n > :: type &L )
294 {
295 FieldType det = FieldType( 1 );
296 for( int i = 0; i < n; ++i )
297 {
298 FieldType &lii = L[ i ][ i ];
299 det *= lii;
300 lii = FieldType( 1 ) / lii;
301 for( int j = 0; j < i; ++j )
302 {
303 FieldType &lij = L[ i ][ j ];
304 FieldType x = lij * L[ j ][ j ];
305 for( int k = j+1; k < i; ++k )
306 x += L[ i ][ k ] * L[ k ][ j ];
307 lij = (-lii) * x;
308 }
309 }
310 return det;
311 }
312
313 // calculates x := L^{-1} x
314 template< int n >
315 static void
316 invLx ( typename Traits :: template Matrix< n, n > :: type &L,
317 typename Traits :: template Vector< n > :: type &x )
318 {
319 for( int i = 0; i < n; ++i )
320 {
321 for( int j = 0; j < i; ++j )
322 x[ i ] -= L[ i ][ j ] * x[ j ];
323 x[ i ] /= L[ i ][ i ];
324 }
325 }
326
327 // calculates x := L^{-T} x
328 template< int n >
329 static void
330 invLTx ( typename Traits :: template Matrix< n, n > :: type &L,
331 typename Traits :: template Vector< n > :: type &x )
332 {
333 for( int i = n; i > 0; --i )
334 {
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 ];
338 }
339 }
340
341 template< int n >
342 static FieldType
343 spdDetA ( const typename Traits :: template Matrix< n, n > :: type &A )
344 {
345 // return A[0][0]*A[1][1]-A[1][0]*A[1][0];
346 typename Traits :: template Matrix< n, n > :: type L;
347 cholesky_L< n >( A, L );
348 return detL< n >( L );
349 }
350
351 template< int n >
352 static FieldType
353 spdInvA ( typename Traits :: template Matrix< n, n > :: type &A )
354 {
355 typename Traits :: template Matrix< n, n > :: type L;
356 cholesky_L< n >( A, L );
357 const FieldType det = invL< n >( L );
358 LTL< n >( L, A );
359 return det;
360 }
361
362 // calculate x := A^{-1} x
363 template< int n >
364 static void
365 spdInvAx ( typename Traits :: template Matrix< n, n > :: type &A,
366 typename Traits :: template Vector< n > :: type &x )
367 {
368 typename Traits :: template Matrix< n, n > :: type L;
369 cholesky_L< n >( A, L );
370 invLx< n >( L, x );
371 invLTx< n >( L, x );
372 }
373
374 template< int m, int n >
375 static FieldType
376 detATA ( const typename Traits :: template Matrix< m, n > :: type &A )
377 {
378 if( m >= n )
379 {
380 typename Traits :: template Matrix< n, n > :: type ata;
381 ATA_L< m, n >( A, ata );
382 return spdDetA< n >( ata );
383 }
384 else
385 return FieldType( 0 );
386 }
387
393 template< int m, int n >
394 static FieldType
395 sqrtDetAAT ( const typename Traits::template Matrix< m, n >::type &A )
396 {
397 // These special cases are here not only for speed reasons:
398 // The general implementation aborts if the matrix is almost singular,
399 // and the special implementation provide a stable way to handle that case.
400 if( (n == 2) && (m == 2) )
401 {
402 // Special implementation for 2x2 matrices: faster and more stable
403 return abs( A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ] );
404 }
405 else if( (n == 3) && (m == 3) )
406 {
407 // Special implementation for 3x3 matrices
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 ] );
412 }
413 else if ( (n == 3) && (m == 2) )
414 {
415 // Special implementation for 2x3 matrices
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);
420 }
421 else if( n >= m )
422 {
423 // General case
424 typename Traits::template Matrix< m, m >::type aat;
425 AAT_L< m, n >( A, aat );
426 return spdDetA< m >( aat );
427 }
428 else
429 return FieldType( 0 );
430 }
431
432 // A^{-1}_L = (A^T A)^{-1} A^T
433 // => A^{-1}_L A = I
434 template< int m, int n >
435 static FieldType
436 leftInvA ( const typename Traits :: template Matrix< m, n > :: type &A,
437 typename Traits :: template Matrix< n, m > :: type &ret )
438 {
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 );
444 return det;
445 }
446
447 template< int m, int n >
448 static void
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 )
452 {
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 );
458 }
459
461 template< int m, int n >
462 static FieldType
463 rightInvA ( const typename Traits :: template Matrix< m, n > :: type &A,
464 typename Traits :: template Matrix< n, m > :: type &ret )
465 {
466 dune_static_assert( (n >= m), "Matrix has no right inverse." );
467 if( (n == 2) && (m == 2) )
468 {
469 const FieldType det = (A[ 0 ][ 0 ]*A[ 1 ][ 1 ] - A[ 1 ][ 0 ]*A[ 0 ][ 1 ]);
470 const FieldType detInv = FieldType( 1 ) / det;
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;
475 return abs( det );
476 }
477 else
478 {
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 );
483 return det;
484 }
485 }
486
487 template< int m, int n >
488 static void
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 )
492 {
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 );
498 }
499 };
500
501 } // namespace GenericGeometry
502
503} // namespace Dune
504
505#endif // #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_MATRIXHELPER_HH
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
dune_static_assert
#define dune_static_assert(COND, MSG)
Helper template so that compilation fails if condition is not true.
Definition: static_assert.hh:79
Dune
Dune namespace.
Definition: alignment.hh:14
static_assert.hh
Fallback implementation of the C++0x static_assert feature.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)