Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
CramerInversionSym.icc
Go to the documentation of this file.
1// @(#)root/smatrix:$Id$
2// Authors: L. Moneta 2005
3
4
5/**********************************************************************
6 * *
7 * Copyright (c) 2005 , LCG ROOT MathLib Team *
8 * *
9 * *
10 **********************************************************************/
11//
12// Cramer optimized inversion for symmetric matrices up to size 5x5.
13// Code from ROOT TMatrixDCramerInv which originates from CLHEP
14// (original author Mark Fischler)
15//
16// Modified by L. Moneta 22/03/07: specialize only until 5x5 (before was up to 6x6)
17// tests show that on 64 machines (like on slc4) it is faster the general method
18//
19
20#ifndef ROOT_Math_CramerInversionSym_icc
21#define ROOT_Math_CramerInversionSym_icc
22
23#ifndef ROOT_Math_Dinv
24#error "Do not use CramerInversionSym.icc directly. #include \"Math/Dinv.h\" instead."
25#endif // ROOT_Math_Dinv
26
27#include <cmath>
28
29
30namespace ROOT {
31
32 namespace Math {
33
34
35
36
37//==============================================================================
38/**
39 Inversion for a 3x3 symmetric matrix
40 */
41template <class T>
43
44 typedef T Scalar;
45
46 // check matrix sizes ??
47
48
49 const Scalar c00 = rhs[4] * rhs[8] - rhs[5] * rhs[5];
50 const Scalar c01 = rhs[5] * rhs[2] - rhs[1] * rhs[8];
51 const Scalar c02 = rhs[1] * rhs[5] - rhs[4] * rhs[2];
52 const Scalar c11 = rhs[8] * rhs[0] - rhs[2] * rhs[2];
53 const Scalar c12 = rhs[2] * rhs[1] - rhs[5] * rhs[0];
54 const Scalar c22 = rhs[0] * rhs[4] - rhs[1] * rhs[1];
55
56 const Scalar t0 = std::abs(rhs[0]);
57 const Scalar t1 = std::abs(rhs[1]);
58 const Scalar t2 = std::abs(rhs[2]);
59
60 Scalar det;
61 Scalar tmp;
62
63 if (t0 >= t1) {
64 if (t2 >= t0) {
65 tmp = rhs[2];
66 det = c12*c01-c11*c02;
67 } else {
68 tmp = rhs[0];
69 det = c11*c22-c12*c12;
70 }
71 } else if (t2 >= t1) {
72 tmp = rhs[2];
73 det = c12*c01-c11*c02;
74 } else {
75 tmp = rhs[1];
76 det = c02*c12-c01*c22;
77 }
78
79 if ( det == 0 || tmp == 0)
80 return false;
81
82 Scalar s = tmp/det;
83// if (determ)
84// *determ = 1./s;
85
86 rhs[0] = s*c00;
87 rhs[1] = s*c01;
88 rhs[2] = s*c02;
89 rhs[4] = s*c11;
90 rhs[5] = s*c12;
91 rhs[8] = s*c22;
92
93 return true;
94}
95
96
97//==============================================================================
98// Inversion for 4x4 matrices
99//==============================================================================
100
101// SFij are indices for a 4x4 symmetric matrix.
102
103#define SF00 0
104#define SF01 1
105#define SF02 2
106#define SF03 3
107
108#define SF10 1
109#define SF11 5
110#define SF12 6
111#define SF13 7
112
113#define SF20 2
114#define SF21 6
115#define SF22 10
116#define SF23 11
117
118#define SF30 3
119#define SF31 7
120#define SF32 11
121#define SF33 15
122
123
124/**
125 Inversion for a 4x4 symmetric matrix
126 */
127template <class T>
129
130 typedef T Scalar;
131
132
133 // Find all NECESSARY 2x2 dets: (14 of them)
134
149
150 // SFind all NECESSSFRY 3x3 dets: (10 of them)
151
172
173 // Find the 4x4 det:
174
177
178// if (determ)
179// *determ = det;
180
181 if ( det == 0 )
182 return false;
183
184 const Scalar oneOverDet = 1.0f / det;
185 const Scalar mn1OverDet = - oneOverDet;
186
191
195
198
200
201 return true;
202}
203
204
205//==============================================================================
206// Inversion for 5x5 matrices
207//==============================================================================
208
209// Mij are indices for a 5x5 matrix.
210
211#define SM00 0
212#define SM01 1
213#define SM02 2
214#define SM03 3
215#define SM04 4
216
217#define SM10 1
218#define SM11 6
219#define SM12 7
220#define SM13 8
221#define SM14 9
222
223#define SM20 2
224#define SM21 7
225#define SM22 12
226#define SM23 13
227#define SM24 14
228
229#define SM30 3
230#define SM31 8
231#define SM32 13
232#define SM33 18
233#define SM34 19
234
235#define SM40 4
236#define SM41 9
237#define SM42 14
238#define SM43 19
239#define SM44 24
240
241/**
242 Inversion for a 5x5 symmetric matrix
243 */
244template <class T>
246
247 typedef T Scalar;
248
249 // Find all NECESSARY 2x2 dets: (25 of them)
250
276
277 // Find all NECESSARY 3x3 dets: (30 of them)
278
309
310 // Find all NECESSARY 4x4 dets: (15 of them)
311
342
343 // Find the 5x5 det:
344
347// if (determ)
348// *determ = det;
349
350 if ( det == 0 )
351 return false;
352
353 const Scalar oneOverDet = 1.0f / det;
354 const Scalar mn1OverDet = - oneOverDet;
355
361
366
370
373
375
376
377 return true;
378}
379
380
381
382 } // namespace Math
383
384} // namespace ROOT
385
386
387// undef Smacros to avoid conflicts
388
389// undef SF 4x4
390#undef SF00
391#undef SF01
392#undef SF02
393#undef SF03
394
395#undef SF10
396#undef SF11
397#undef SF12
398#undef SF13
399
400#undef SF20
401#undef SF21
402#undef SF22
403#undef SF23
404
405#undef SF30
406#undef SF31
407#undef SF32
408#undef SF33
409
410// undef SM 5x5
411#undef SM00
412#undef SM01
413#undef SM02
414#undef SM03
415#undef SM04
416
417#undef SM10
418#undef SM11
419#undef SM12
420#undef SM13
421#undef SM14
422
423#undef SM20
424#undef SM21
425#undef SM22
426#undef SM23
427#undef SM24
428
429#undef SM30
430#undef SM31
431#undef SM32
432#undef SM33
433#undef SM34
434
435#undef SM40
436#undef SM41
437#undef SM42
438#undef SM43
439#undef SM44
440
441
442
443
444#endif
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
#define SM33
#define SM12
#define SM13
#define SM40
#define SF03
#define SM14
#define SF02
#define SF31
#define SF00
#define SF13
#define SF22
#define SF32
#define SM23
#define SF10
#define SM22
#define SF33
#define SM10
#define SM00
#define SF21
#define SM30
#define SF12
#define SM01
#define SF23
#define SM02
#define SF20
#define SM41
#define SM21
#define SM32
#define SM11
#define SM20
#define SM44
#define SM04
#define SF30
#define SF11
#define SM43
#define SM31
#define SM03
#define SM34
#define SM42
#define SF01
#define SM24
static bool Dinv(MatrixRep &rhs)
Definition Dinv.h:148
Namespace for new Math classes and functions.
Rotation3D::Scalar Scalar
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
auto * t1
Definition textangle.C:20