Logo ROOT   6.18/05
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
30namespace ROOT {
31
32 namespace Math {
33
34
35
36//==============================================================================
37// Inversion for 3x3 matrices
38//==============================================================================
39
40/**
41 Inversion for a 3x3 matrix
42 */
43template <class MatrixRep>
44bool 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 */
134template <class MatrixRep>
135bool 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 */
276template <class MatrixRep>
277bool 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
F01
#define F01
Definition: CramerInversion.icc:112
F22
#define F22
Definition: CramerInversion.icc:123
F00
#define F00
Definition: CramerInversion.icc:111
M00
#define M00
Definition: CramerInversion.icc:242
F02
#define F02
Definition: CramerInversion.icc:113
M34
#define M34
Definition: CramerInversion.icc:264
M03
#define M03
Definition: CramerInversion.icc:245
F31
#define F31
Definition: CramerInversion.icc:127
F32
#define F32
Definition: CramerInversion.icc:128
M10
#define M10
Definition: CramerInversion.icc:248
F21
#define F21
Definition: CramerInversion.icc:122
F03
#define F03
Definition: CramerInversion.icc:114
M41
#define M41
Definition: CramerInversion.icc:267
M20
#define M20
Definition: CramerInversion.icc:254
M42
#define M42
Definition: CramerInversion.icc:268
M13
#define M13
Definition: CramerInversion.icc:251
F30
#define F30
Definition: CramerInversion.icc:126
M33
#define M33
Definition: CramerInversion.icc:263
M04
#define M04
Definition: CramerInversion.icc:246
M23
#define M23
Definition: CramerInversion.icc:257
M11
#define M11
Definition: CramerInversion.icc:249
M44
#define M44
Definition: CramerInversion.icc:270
F11
#define F11
Definition: CramerInversion.icc:117
M31
#define M31
Definition: CramerInversion.icc:261
F10
#define F10
Definition: CramerInversion.icc:116
M02
#define M02
Definition: CramerInversion.icc:244
M21
#define M21
Definition: CramerInversion.icc:255
M01
#define M01
Definition: CramerInversion.icc:243
F33
#define F33
Definition: CramerInversion.icc:129
M24
#define M24
Definition: CramerInversion.icc:258
F20
#define F20
Definition: CramerInversion.icc:121
M12
#define M12
Definition: CramerInversion.icc:250
M22
#define M22
Definition: CramerInversion.icc:256
F12
#define F12
Definition: CramerInversion.icc:118
M43
#define M43
Definition: CramerInversion.icc:269
F23
#define F23
Definition: CramerInversion.icc:124
M14
#define M14
Definition: CramerInversion.icc:252
M30
#define M30
Definition: CramerInversion.icc:260
M40
#define M40
Definition: CramerInversion.icc:266
M32
#define M32
Definition: CramerInversion.icc:262
F13
#define F13
Definition: CramerInversion.icc:119
ROOT::Math::FastInverter::Dinv
static bool Dinv(MatrixRep &rhs)
Definition: Dinv.h:148
Math
Namespace for new Math classes and functions.
ROOT::Math::Scalar
Rotation3D::Scalar Scalar
Definition: Rotation3DxAxial.cxx:69
ROOT
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
t1
auto * t1
Definition: textangle.C:20