Logo ROOT   6.12/07
Reference Guide
CramerInversion.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 optmized inversion for 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_CramerInversion_icc
21 #define ROOT_Math_CramerInversion_icc
22 
23 #ifndef ROOT_Math_Dinv
24 #error "Do not use CramerInversion.icc directly. #include \"Math/Dinv.h\" instead."
25 #endif // ROOT_Math_Dinv
26 
27 #include <cmath>
28 
29 
30 namespace ROOT {
31 
32  namespace Math {
33 
34 
35 
36 //==============================================================================
37 // Inversion for 3x3 matrices
38 //==============================================================================
39 
40 /**
41  Inversion for a 3x3 matrix
42  */
43 template <class MatrixRep>
44 bool FastInverter<3>::Dinv(MatrixRep & rhs) {
45 
46  typedef typename MatrixRep::value_type Scalar;
47 
48  // check matrix sizes ??
49 
50  // Scalar * pM = rhs.Array();
51 
52  const Scalar c00 = rhs[4] * rhs[8] - rhs[5] * rhs[7];
53  const Scalar c01 = rhs[5] * rhs[6] - rhs[3] * rhs[8];
54  const Scalar c02 = rhs[3] * rhs[7] - rhs[4] * rhs[6];
55  const Scalar c10 = rhs[7] * rhs[2] - rhs[8] * rhs[1];
56  const Scalar c11 = rhs[8] * rhs[0] - rhs[6] * rhs[2];
57  const Scalar c12 = rhs[6] * rhs[1] - rhs[7] * rhs[0];
58  const Scalar c20 = rhs[1] * rhs[5] - rhs[2] * rhs[4];
59  const Scalar c21 = rhs[2] * rhs[3] - rhs[0] * rhs[5];
60  const Scalar c22 = rhs[0] * rhs[4] - rhs[1] * rhs[3];
61 
62  const Scalar t0 = std::abs(rhs[0]);
63  const Scalar t1 = std::abs(rhs[3]);
64  const Scalar t2 = std::abs(rhs[6]);
65  Scalar det;
66  Scalar tmp;
67  if (t0 >= t1) {
68  if (t2 >= t0) {
69  tmp = rhs[6];
70  det = c12*c01-c11*c02;
71  } else {
72  tmp = rhs[0];
73  det = c11*c22-c12*c21;
74  }
75  } else if (t2 >= t1) {
76  tmp = rhs[6];
77  det = c12*c01-c11*c02;
78  } else {
79  tmp = rhs[3];
80  det = c02*c21-c01*c22;
81  }
82 
83  if ( det == 0 || tmp == 0) {
84  return false;
85  }
86 
87  const Scalar s = tmp/det;
88 
89 // if (determ)
90 // *determ = 1./s;
91 
92  rhs[0] = s*c00;
93  rhs[1] = s*c10;
94  rhs[2] = s*c20;
95  rhs[3] = s*c01;
96  rhs[4] = s*c11;
97  rhs[5] = s*c21;
98  rhs[6] = s*c02;
99  rhs[7] = s*c12;
100  rhs[8] = s*c22;
101 
102  return true;
103 }
104 
105 
106 //==============================================================================
107 // Inversion for 4x4 matrices
108 //==============================================================================
109 // Fij are indices for a 4x4 matrix.
110 
111 #define F00 0
112 #define F01 1
113 #define F02 2
114 #define F03 3
115 
116 #define F10 4
117 #define F11 5
118 #define F12 6
119 #define F13 7
120 
121 #define F20 8
122 #define F21 9
123 #define F22 10
124 #define F23 11
125 
126 #define F30 12
127 #define F31 13
128 #define F32 14
129 #define F33 15
130 
131 /**
132  Inversion for a 4x4 matrix
133  */
134 template <class MatrixRep>
135 bool FastInverter<4>::Dinv(MatrixRep & rhs) {
136 
137  typedef typename MatrixRep::value_type Scalar;
138 
139  // check matrix sizes ??
140 
141  // Scalar * pM = rhs.Array();
142 
143  // Find all NECESSARY 2x2 dets: (18 of them)
144 
145  const Scalar det2_12_01 = rhs[F10]*rhs[F21] - rhs[F11]*rhs[F20];
146  const Scalar det2_12_02 = rhs[F10]*rhs[F22] - rhs[F12]*rhs[F20];
147  const Scalar det2_12_03 = rhs[F10]*rhs[F23] - rhs[F13]*rhs[F20];
148  const Scalar det2_12_13 = rhs[F11]*rhs[F23] - rhs[F13]*rhs[F21];
149  const Scalar det2_12_23 = rhs[F12]*rhs[F23] - rhs[F13]*rhs[F22];
150  const Scalar det2_12_12 = rhs[F11]*rhs[F22] - rhs[F12]*rhs[F21];
151  const Scalar det2_13_01 = rhs[F10]*rhs[F31] - rhs[F11]*rhs[F30];
152  const Scalar det2_13_02 = rhs[F10]*rhs[F32] - rhs[F12]*rhs[F30];
153  const Scalar det2_13_03 = rhs[F10]*rhs[F33] - rhs[F13]*rhs[F30];
154  const Scalar det2_13_12 = rhs[F11]*rhs[F32] - rhs[F12]*rhs[F31];
155  const Scalar det2_13_13 = rhs[F11]*rhs[F33] - rhs[F13]*rhs[F31];
156  const Scalar det2_13_23 = rhs[F12]*rhs[F33] - rhs[F13]*rhs[F32];
157  const Scalar det2_23_01 = rhs[F20]*rhs[F31] - rhs[F21]*rhs[F30];
158  const Scalar det2_23_02 = rhs[F20]*rhs[F32] - rhs[F22]*rhs[F30];
159  const Scalar det2_23_03 = rhs[F20]*rhs[F33] - rhs[F23]*rhs[F30];
160  const Scalar det2_23_12 = rhs[F21]*rhs[F32] - rhs[F22]*rhs[F31];
161  const Scalar det2_23_13 = rhs[F21]*rhs[F33] - rhs[F23]*rhs[F31];
162  const Scalar det2_23_23 = rhs[F22]*rhs[F33] - rhs[F23]*rhs[F32];
163 
164  // Find all NECESSARY 3x3 dets: (16 of them)
165 
166  const Scalar det3_012_012 = rhs[F00]*det2_12_12 - rhs[F01]*det2_12_02
167  + rhs[F02]*det2_12_01;
168  const Scalar det3_012_013 = rhs[F00]*det2_12_13 - rhs[F01]*det2_12_03
169  + rhs[F03]*det2_12_01;
170  const Scalar det3_012_023 = rhs[F00]*det2_12_23 - rhs[F02]*det2_12_03
171  + rhs[F03]*det2_12_02;
172  const Scalar det3_012_123 = rhs[F01]*det2_12_23 - rhs[F02]*det2_12_13
173  + rhs[F03]*det2_12_12;
174  const Scalar det3_013_012 = rhs[F00]*det2_13_12 - rhs[F01]*det2_13_02
175  + rhs[F02]*det2_13_01;
176  const Scalar det3_013_013 = rhs[F00]*det2_13_13 - rhs[F01]*det2_13_03
177  + rhs[F03]*det2_13_01;
178  const Scalar det3_013_023 = rhs[F00]*det2_13_23 - rhs[F02]*det2_13_03
179  + rhs[F03]*det2_13_02;
180  const Scalar det3_013_123 = rhs[F01]*det2_13_23 - rhs[F02]*det2_13_13
181  + rhs[F03]*det2_13_12;
182  const Scalar det3_023_012 = rhs[F00]*det2_23_12 - rhs[F01]*det2_23_02
183  + rhs[F02]*det2_23_01;
184  const Scalar det3_023_013 = rhs[F00]*det2_23_13 - rhs[F01]*det2_23_03
185  + rhs[F03]*det2_23_01;
186  const Scalar det3_023_023 = rhs[F00]*det2_23_23 - rhs[F02]*det2_23_03
187  + rhs[F03]*det2_23_02;
188  const Scalar det3_023_123 = rhs[F01]*det2_23_23 - rhs[F02]*det2_23_13
189  + rhs[F03]*det2_23_12;
190  const Scalar det3_123_012 = rhs[F10]*det2_23_12 - rhs[F11]*det2_23_02
191  + rhs[F12]*det2_23_01;
192  const Scalar det3_123_013 = rhs[F10]*det2_23_13 - rhs[F11]*det2_23_03
193  + rhs[F13]*det2_23_01;
194  const Scalar det3_123_023 = rhs[F10]*det2_23_23 - rhs[F12]*det2_23_03
195  + rhs[F13]*det2_23_02;
196  const Scalar det3_123_123 = rhs[F11]*det2_23_23 - rhs[F12]*det2_23_13
197  + rhs[F13]*det2_23_12;
198 
199  // Find the 4x4 det:
200 
201  const Scalar det = rhs[F00]*det3_123_123 - rhs[F01]*det3_123_023
202  + rhs[F02]*det3_123_013 - rhs[F03]*det3_123_012;
203 
204 // if (determ)
205 // *determ = det;
206 
207  if ( det == 0 ) {
208  return false;
209  }
210 
211  // use 1.0f to remove warning C4244 on Windows when using float
212  const Scalar oneOverDet = 1.0f / det;
213  const Scalar mn1OverDet = - oneOverDet;
214 
215  rhs[F00] = det3_123_123 * oneOverDet;
216  rhs[F01] = det3_023_123 * mn1OverDet;
217  rhs[F02] = det3_013_123 * oneOverDet;
218  rhs[F03] = det3_012_123 * mn1OverDet;
219 
220  rhs[F10] = det3_123_023 * mn1OverDet;
221  rhs[F11] = det3_023_023 * oneOverDet;
222  rhs[F12] = det3_013_023 * mn1OverDet;
223  rhs[F13] = det3_012_023 * oneOverDet;
224 
225  rhs[F20] = det3_123_013 * oneOverDet;
226  rhs[F21] = det3_023_013 * mn1OverDet;
227  rhs[F22] = det3_013_013 * oneOverDet;
228  rhs[F23] = det3_012_013 * mn1OverDet;
229 
230  rhs[F30] = det3_123_012 * mn1OverDet;
231  rhs[F31] = det3_023_012 * oneOverDet;
232  rhs[F32] = det3_013_012 * mn1OverDet;
233  rhs[F33] = det3_012_012 * oneOverDet;
234 
235  return true;
236 }
237 
238 //==============================================================================
239 // Inversion for 5x5 matrices
240 //==============================================================================
241 // Mij are indices for a 5x5 matrix.
242 #define M00 0
243 #define M01 1
244 #define M02 2
245 #define M03 3
246 #define M04 4
247 
248 #define M10 5
249 #define M11 6
250 #define M12 7
251 #define M13 8
252 #define M14 9
253 
254 #define M20 10
255 #define M21 11
256 #define M22 12
257 #define M23 13
258 #define M24 14
259 
260 #define M30 15
261 #define M31 16
262 #define M32 17
263 #define M33 18
264 #define M34 19
265 
266 #define M40 20
267 #define M41 21
268 #define M42 22
269 #define M43 23
270 #define M44 24
271 
272 
273 /**
274  Inversion for a 5x5 matrix
275  */
276 template <class MatrixRep>
277 bool FastInverter<5>::Dinv(MatrixRep & rhs) {
278 
279  typedef typename MatrixRep::value_type Scalar;
280 
281  // check matrix sizes ??
282 
283  // Scalar * pM = rhs.Array();
284 
285 
286  // Find all NECESSARY 2x2 dets: (30 of them)
287 
288  const Scalar det2_23_01 = rhs[M20]*rhs[M31] - rhs[M21]*rhs[M30];
289  const Scalar det2_23_02 = rhs[M20]*rhs[M32] - rhs[M22]*rhs[M30];
290  const Scalar det2_23_03 = rhs[M20]*rhs[M33] - rhs[M23]*rhs[M30];
291  const Scalar det2_23_04 = rhs[M20]*rhs[M34] - rhs[M24]*rhs[M30];
292  const Scalar det2_23_12 = rhs[M21]*rhs[M32] - rhs[M22]*rhs[M31];
293  const Scalar det2_23_13 = rhs[M21]*rhs[M33] - rhs[M23]*rhs[M31];
294  const Scalar det2_23_14 = rhs[M21]*rhs[M34] - rhs[M24]*rhs[M31];
295  const Scalar det2_23_23 = rhs[M22]*rhs[M33] - rhs[M23]*rhs[M32];
296  const Scalar det2_23_24 = rhs[M22]*rhs[M34] - rhs[M24]*rhs[M32];
297  const Scalar det2_23_34 = rhs[M23]*rhs[M34] - rhs[M24]*rhs[M33];
298  const Scalar det2_24_01 = rhs[M20]*rhs[M41] - rhs[M21]*rhs[M40];
299  const Scalar det2_24_02 = rhs[M20]*rhs[M42] - rhs[M22]*rhs[M40];
300  const Scalar det2_24_03 = rhs[M20]*rhs[M43] - rhs[M23]*rhs[M40];
301  const Scalar det2_24_04 = rhs[M20]*rhs[M44] - rhs[M24]*rhs[M40];
302  const Scalar det2_24_12 = rhs[M21]*rhs[M42] - rhs[M22]*rhs[M41];
303  const Scalar det2_24_13 = rhs[M21]*rhs[M43] - rhs[M23]*rhs[M41];
304  const Scalar det2_24_14 = rhs[M21]*rhs[M44] - rhs[M24]*rhs[M41];
305  const Scalar det2_24_23 = rhs[M22]*rhs[M43] - rhs[M23]*rhs[M42];
306  const Scalar det2_24_24 = rhs[M22]*rhs[M44] - rhs[M24]*rhs[M42];
307  const Scalar det2_24_34 = rhs[M23]*rhs[M44] - rhs[M24]*rhs[M43];
308  const Scalar det2_34_01 = rhs[M30]*rhs[M41] - rhs[M31]*rhs[M40];
309  const Scalar det2_34_02 = rhs[M30]*rhs[M42] - rhs[M32]*rhs[M40];
310  const Scalar det2_34_03 = rhs[M30]*rhs[M43] - rhs[M33]*rhs[M40];
311  const Scalar det2_34_04 = rhs[M30]*rhs[M44] - rhs[M34]*rhs[M40];
312  const Scalar det2_34_12 = rhs[M31]*rhs[M42] - rhs[M32]*rhs[M41];
313  const Scalar det2_34_13 = rhs[M31]*rhs[M43] - rhs[M33]*rhs[M41];
314  const Scalar det2_34_14 = rhs[M31]*rhs[M44] - rhs[M34]*rhs[M41];
315  const Scalar det2_34_23 = rhs[M32]*rhs[M43] - rhs[M33]*rhs[M42];
316  const Scalar det2_34_24 = rhs[M32]*rhs[M44] - rhs[M34]*rhs[M42];
317  const Scalar det2_34_34 = rhs[M33]*rhs[M44] - rhs[M34]*rhs[M43];
318 
319  // Find all NECESSARY 3x3 dets: (40 of them)
320 
321  const Scalar det3_123_012 = rhs[M10]*det2_23_12 - rhs[M11]*det2_23_02 + rhs[M12]*det2_23_01;
322  const Scalar det3_123_013 = rhs[M10]*det2_23_13 - rhs[M11]*det2_23_03 + rhs[M13]*det2_23_01;
323  const Scalar det3_123_014 = rhs[M10]*det2_23_14 - rhs[M11]*det2_23_04 + rhs[M14]*det2_23_01;
324  const Scalar det3_123_023 = rhs[M10]*det2_23_23 - rhs[M12]*det2_23_03 + rhs[M13]*det2_23_02;
325  const Scalar det3_123_024 = rhs[M10]*det2_23_24 - rhs[M12]*det2_23_04 + rhs[M14]*det2_23_02;
326  const Scalar det3_123_034 = rhs[M10]*det2_23_34 - rhs[M13]*det2_23_04 + rhs[M14]*det2_23_03;
327  const Scalar det3_123_123 = rhs[M11]*det2_23_23 - rhs[M12]*det2_23_13 + rhs[M13]*det2_23_12;
328  const Scalar det3_123_124 = rhs[M11]*det2_23_24 - rhs[M12]*det2_23_14 + rhs[M14]*det2_23_12;
329  const Scalar det3_123_134 = rhs[M11]*det2_23_34 - rhs[M13]*det2_23_14 + rhs[M14]*det2_23_13;
330  const Scalar det3_123_234 = rhs[M12]*det2_23_34 - rhs[M13]*det2_23_24 + rhs[M14]*det2_23_23;
331  const Scalar det3_124_012 = rhs[M10]*det2_24_12 - rhs[M11]*det2_24_02 + rhs[M12]*det2_24_01;
332  const Scalar det3_124_013 = rhs[M10]*det2_24_13 - rhs[M11]*det2_24_03 + rhs[M13]*det2_24_01;
333  const Scalar det3_124_014 = rhs[M10]*det2_24_14 - rhs[M11]*det2_24_04 + rhs[M14]*det2_24_01;
334  const Scalar det3_124_023 = rhs[M10]*det2_24_23 - rhs[M12]*det2_24_03 + rhs[M13]*det2_24_02;
335  const Scalar det3_124_024 = rhs[M10]*det2_24_24 - rhs[M12]*det2_24_04 + rhs[M14]*det2_24_02;
336  const Scalar det3_124_034 = rhs[M10]*det2_24_34 - rhs[M13]*det2_24_04 + rhs[M14]*det2_24_03;
337  const Scalar det3_124_123 = rhs[M11]*det2_24_23 - rhs[M12]*det2_24_13 + rhs[M13]*det2_24_12;
338  const Scalar det3_124_124 = rhs[M11]*det2_24_24 - rhs[M12]*det2_24_14 + rhs[M14]*det2_24_12;
339  const Scalar det3_124_134 = rhs[M11]*det2_24_34 - rhs[M13]*det2_24_14 + rhs[M14]*det2_24_13;
340  const Scalar det3_124_234 = rhs[M12]*det2_24_34 - rhs[M13]*det2_24_24 + rhs[M14]*det2_24_23;
341  const Scalar det3_134_012 = rhs[M10]*det2_34_12 - rhs[M11]*det2_34_02 + rhs[M12]*det2_34_01;
342  const Scalar det3_134_013 = rhs[M10]*det2_34_13 - rhs[M11]*det2_34_03 + rhs[M13]*det2_34_01;
343  const Scalar det3_134_014 = rhs[M10]*det2_34_14 - rhs[M11]*det2_34_04 + rhs[M14]*det2_34_01;
344  const Scalar det3_134_023 = rhs[M10]*det2_34_23 - rhs[M12]*det2_34_03 + rhs[M13]*det2_34_02;
345  const Scalar det3_134_024 = rhs[M10]*det2_34_24 - rhs[M12]*det2_34_04 + rhs[M14]*det2_34_02;
346  const Scalar det3_134_034 = rhs[M10]*det2_34_34 - rhs[M13]*det2_34_04 + rhs[M14]*det2_34_03;
347  const Scalar det3_134_123 = rhs[M11]*det2_34_23 - rhs[M12]*det2_34_13 + rhs[M13]*det2_34_12;
348  const Scalar det3_134_124 = rhs[M11]*det2_34_24 - rhs[M12]*det2_34_14 + rhs[M14]*det2_34_12;
349  const Scalar det3_134_134 = rhs[M11]*det2_34_34 - rhs[M13]*det2_34_14 + rhs[M14]*det2_34_13;
350  const Scalar det3_134_234 = rhs[M12]*det2_34_34 - rhs[M13]*det2_34_24 + rhs[M14]*det2_34_23;
351  const Scalar det3_234_012 = rhs[M20]*det2_34_12 - rhs[M21]*det2_34_02 + rhs[M22]*det2_34_01;
352  const Scalar det3_234_013 = rhs[M20]*det2_34_13 - rhs[M21]*det2_34_03 + rhs[M23]*det2_34_01;
353  const Scalar det3_234_014 = rhs[M20]*det2_34_14 - rhs[M21]*det2_34_04 + rhs[M24]*det2_34_01;
354  const Scalar det3_234_023 = rhs[M20]*det2_34_23 - rhs[M22]*det2_34_03 + rhs[M23]*det2_34_02;
355  const Scalar det3_234_024 = rhs[M20]*det2_34_24 - rhs[M22]*det2_34_04 + rhs[M24]*det2_34_02;
356  const Scalar det3_234_034 = rhs[M20]*det2_34_34 - rhs[M23]*det2_34_04 + rhs[M24]*det2_34_03;
357  const Scalar det3_234_123 = rhs[M21]*det2_34_23 - rhs[M22]*det2_34_13 + rhs[M23]*det2_34_12;
358  const Scalar det3_234_124 = rhs[M21]*det2_34_24 - rhs[M22]*det2_34_14 + rhs[M24]*det2_34_12;
359  const Scalar det3_234_134 = rhs[M21]*det2_34_34 - rhs[M23]*det2_34_14 + rhs[M24]*det2_34_13;
360  const Scalar det3_234_234 = rhs[M22]*det2_34_34 - rhs[M23]*det2_34_24 + rhs[M24]*det2_34_23;
361 
362  // Find all NECESSARY 4x4 dets: (25 of them)
363 
364  const Scalar det4_0123_0123 = rhs[M00]*det3_123_123 - rhs[M01]*det3_123_023
365  + rhs[M02]*det3_123_013 - rhs[M03]*det3_123_012;
366  const Scalar det4_0123_0124 = rhs[M00]*det3_123_124 - rhs[M01]*det3_123_024
367  + rhs[M02]*det3_123_014 - rhs[M04]*det3_123_012;
368  const Scalar det4_0123_0134 = rhs[M00]*det3_123_134 - rhs[M01]*det3_123_034
369  + rhs[M03]*det3_123_014 - rhs[M04]*det3_123_013;
370  const Scalar det4_0123_0234 = rhs[M00]*det3_123_234 - rhs[M02]*det3_123_034
371  + rhs[M03]*det3_123_024 - rhs[M04]*det3_123_023;
372  const Scalar det4_0123_1234 = rhs[M01]*det3_123_234 - rhs[M02]*det3_123_134
373  + rhs[M03]*det3_123_124 - rhs[M04]*det3_123_123;
374  const Scalar det4_0124_0123 = rhs[M00]*det3_124_123 - rhs[M01]*det3_124_023
375  + rhs[M02]*det3_124_013 - rhs[M03]*det3_124_012;
376  const Scalar det4_0124_0124 = rhs[M00]*det3_124_124 - rhs[M01]*det3_124_024
377  + rhs[M02]*det3_124_014 - rhs[M04]*det3_124_012;
378  const Scalar det4_0124_0134 = rhs[M00]*det3_124_134 - rhs[M01]*det3_124_034
379  + rhs[M03]*det3_124_014 - rhs[M04]*det3_124_013;
380  const Scalar det4_0124_0234 = rhs[M00]*det3_124_234 - rhs[M02]*det3_124_034
381  + rhs[M03]*det3_124_024 - rhs[M04]*det3_124_023;
382  const Scalar det4_0124_1234 = rhs[M01]*det3_124_234 - rhs[M02]*det3_124_134
383  + rhs[M03]*det3_124_124 - rhs[M04]*det3_124_123;
384  const Scalar det4_0134_0123 = rhs[M00]*det3_134_123 - rhs[M01]*det3_134_023
385  + rhs[M02]*det3_134_013 - rhs[M03]*det3_134_012;
386  const Scalar det4_0134_0124 = rhs[M00]*det3_134_124 - rhs[M01]*det3_134_024
387  + rhs[M02]*det3_134_014 - rhs[M04]*det3_134_012;
388  const Scalar det4_0134_0134 = rhs[M00]*det3_134_134 - rhs[M01]*det3_134_034
389  + rhs[M03]*det3_134_014 - rhs[M04]*det3_134_013;
390  const Scalar det4_0134_0234 = rhs[M00]*det3_134_234 - rhs[M02]*det3_134_034
391  + rhs[M03]*det3_134_024 - rhs[M04]*det3_134_023;
392  const Scalar det4_0134_1234 = rhs[M01]*det3_134_234 - rhs[M02]*det3_134_134
393  + rhs[M03]*det3_134_124 - rhs[M04]*det3_134_123;
394  const Scalar det4_0234_0123 = rhs[M00]*det3_234_123 - rhs[M01]*det3_234_023
395  + rhs[M02]*det3_234_013 - rhs[M03]*det3_234_012;
396  const Scalar det4_0234_0124 = rhs[M00]*det3_234_124 - rhs[M01]*det3_234_024
397  + rhs[M02]*det3_234_014 - rhs[M04]*det3_234_012;
398  const Scalar det4_0234_0134 = rhs[M00]*det3_234_134 - rhs[M01]*det3_234_034
399  + rhs[M03]*det3_234_014 - rhs[M04]*det3_234_013;
400  const Scalar det4_0234_0234 = rhs[M00]*det3_234_234 - rhs[M02]*det3_234_034
401  + rhs[M03]*det3_234_024 - rhs[M04]*det3_234_023;
402  const Scalar det4_0234_1234 = rhs[M01]*det3_234_234 - rhs[M02]*det3_234_134
403  + rhs[M03]*det3_234_124 - rhs[M04]*det3_234_123;
404  const Scalar det4_1234_0123 = rhs[M10]*det3_234_123 - rhs[M11]*det3_234_023
405  + rhs[M12]*det3_234_013 - rhs[M13]*det3_234_012;
406  const Scalar det4_1234_0124 = rhs[M10]*det3_234_124 - rhs[M11]*det3_234_024
407  + rhs[M12]*det3_234_014 - rhs[M14]*det3_234_012;
408  const Scalar det4_1234_0134 = rhs[M10]*det3_234_134 - rhs[M11]*det3_234_034
409  + rhs[M13]*det3_234_014 - rhs[M14]*det3_234_013;
410  const Scalar det4_1234_0234 = rhs[M10]*det3_234_234 - rhs[M12]*det3_234_034
411  + rhs[M13]*det3_234_024 - rhs[M14]*det3_234_023;
412  const Scalar det4_1234_1234 = rhs[M11]*det3_234_234 - rhs[M12]*det3_234_134
413  + rhs[M13]*det3_234_124 - rhs[M14]*det3_234_123;
414 
415  // Find the 5x5 det:
416 
417  const Scalar det = rhs[M00]*det4_1234_1234 - rhs[M01]*det4_1234_0234 + rhs[M02]*det4_1234_0134
418  - rhs[M03]*det4_1234_0124 + rhs[M04]*det4_1234_0123;
419 
420 // if (determ)
421 // *determ = det;
422 
423  if ( det == 0 ) {
424  //Error("Inv5x5","matrix is singular");
425  //m.Invalidate();
426  return false;
427  }
428 
429  const Scalar oneOverDet = 1.0f / det;
430  const Scalar mn1OverDet = - oneOverDet;
431 
432  rhs[M00] = det4_1234_1234 * oneOverDet;
433  rhs[M01] = det4_0234_1234 * mn1OverDet;
434  rhs[M02] = det4_0134_1234 * oneOverDet;
435  rhs[M03] = det4_0124_1234 * mn1OverDet;
436  rhs[M04] = det4_0123_1234 * oneOverDet;
437 
438  rhs[M10] = det4_1234_0234 * mn1OverDet;
439  rhs[M11] = det4_0234_0234 * oneOverDet;
440  rhs[M12] = det4_0134_0234 * mn1OverDet;
441  rhs[M13] = det4_0124_0234 * oneOverDet;
442  rhs[M14] = det4_0123_0234 * mn1OverDet;
443 
444  rhs[M20] = det4_1234_0134 * oneOverDet;
445  rhs[M21] = det4_0234_0134 * mn1OverDet;
446  rhs[M22] = det4_0134_0134 * oneOverDet;
447  rhs[M23] = det4_0124_0134 * mn1OverDet;
448  rhs[M24] = det4_0123_0134 * oneOverDet;
449 
450  rhs[M30] = det4_1234_0124 * mn1OverDet;
451  rhs[M31] = det4_0234_0124 * oneOverDet;
452  rhs[M32] = det4_0134_0124 * mn1OverDet;
453  rhs[M33] = det4_0124_0124 * oneOverDet;
454  rhs[M34] = det4_0123_0124 * mn1OverDet;
455 
456  rhs[M40] = det4_1234_0123 * oneOverDet;
457  rhs[M41] = det4_0234_0123 * mn1OverDet;
458  rhs[M42] = det4_0134_0123 * oneOverDet;
459  rhs[M43] = det4_0124_0123 * mn1OverDet;
460  rhs[M44] = det4_0123_0123 * oneOverDet;
461 
462  return true;
463 }
464 
465 
466  } // namespace Math
467 
468 } // namespace ROOT
469 
470 
471 
472 
473 // undef macros to avoid conflicts
474 
475 // 4x4 indices
476 
477 #undef F00
478 #undef F01
479 #undef F02
480 #undef F03
481 
482 #undef F10
483 #undef F11
484 #undef F12
485 #undef F13
486 
487 #undef F20
488 #undef F21
489 #undef F22
490 #undef F23
491 
492 #undef F30
493 #undef F31
494 #undef F32
495 #undef F33
496 
497 // undef 5x5 indices
498 
499 #undef M00
500 #undef M01
501 #undef M02
502 #undef M03
503 #undef M04
504 
505 #undef M10
506 #undef M11
507 #undef M12
508 #undef M13
509 #undef M14
510 
511 #undef M20
512 #undef M21
513 #undef M22
514 #undef M23
515 #undef M24
516 
517 #undef M30
518 #undef M31
519 #undef M32
520 #undef M33
521 #undef M34
522 
523 #undef M40
524 #undef M41
525 #undef M42
526 #undef M43
527 #undef M44
528 
529 
530 #endif
#define M21
#define M11
#define M00
#define M43
#define M01
#define M32
#define M31
#define F03
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
#define M14
#define M40
#define F20
#define M41
static bool Dinv(MatrixRep &rhs)
Definition: Dinv.h:148
#define M13
#define M12
#define F10
#define F01
#define F12
#define F31
#define M23
#define F22
#define F30
#define F32
#define M02
#define M33
#define M30
#define M42
#define M10
#define M22
auto * t1
Definition: textangle.C:20
#define M20
#define M04
static constexpr double s
#define M03
Namespace for new Math classes and functions.
#define F21
#define F13
#define M24
#define M44
#define F00
#define F02
#define F33
#define F11
Rotation3D::Scalar Scalar
#define F23
#define M34