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