Logo ROOT  

Reference Guide
TMatrixTSymCramerInv.cxx
Go to the documentation of this file.
1// @(#)root/base:$Id$
2// Authors: Fons Rademakers, Eddy Offermann Oct 2004
3
4/*************************************************************************
5 * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class TMatrixTSymCramerInv
13 \ingroup Matrix
14
15 TMatrixTSymCramerInv
16
17 Encapsulate templates of Cramer Inversion routines.
18
19 The 4x4, 5x5 and 6x6 are adapted from routines written by
20 Mark Fischler and Steven Haywood as part of the CLHEP package
21
22 Although for sizes <= 6x6 the Cramer Inversion has a gain in speed
23 compared to factorization schemes (like LU) , one pays a price in
24 accuracy .
25
26 For Example:
27~~~
28 H * H^-1 = U, where H is a 5x5 Hilbert matrix
29 U is a 5x5 Unity matrix
30
31 LU : |U_jk| < 10e-13 for j!=k
32 Cramer: |U_jk| < 10e-7 for j!=k
33~~~
34
35 however Cramer algorithm is about 10 (!) times faster
36*/
37
38#include "TMatrixTSymCramerInv.h"
39
40#include "TError.h"
41
42#if !defined(R__SOLARIS) && !defined(R__ACC) && !defined(R__FBSD)
43NamespaceImp(TMatrixTSymCramerInv);
44#endif
45
46////////////////////////////////////////////////////////////////////////////////
47
48template<class Element>
49Bool_t TMatrixTSymCramerInv::Inv2x2(TMatrixTSym<Element> &m,Double_t *determ)
50{
51 if (m.GetNrows() != 2) {
52 Error("Inv2x2","matrix should be square 2x2");
53 return kFALSE;
54 }
55
56 Element *pM = m.GetMatrixArray();
57
58 const Double_t det = pM[0] * pM[3] - pM[1] * pM[1];
59
60 if (determ)
61 *determ = det;
62
63 if ( det == 0 ) {
64 Error("Inv2x2","matrix is singular");
65 return kFALSE;
66 }
67
68 const Double_t tmp1 = pM[3] / det;
69 pM[3] = pM[0] / det;
70 pM[2] = pM[1] = -pM[1] / det;
71 pM[0] = tmp1;
72
73 return kTRUE;
74}
75
76////////////////////////////////////////////////////////////////////////////////
77
78template<class Element>
79Bool_t TMatrixTSymCramerInv::Inv3x3(TMatrixTSym<Element> &m,Double_t *determ)
80{
81 if (m.GetNrows() != 3) {
82 Error("Inv3x3","matrix should be square 3x3");
83 return kFALSE;
84 }
85
86 Element *pM = m.GetMatrixArray();
87
88 const Double_t c00 = pM[4] * pM[8] - pM[5] * pM[5];
89 const Double_t c01 = pM[5] * pM[2] - pM[1] * pM[8];
90 const Double_t c02 = pM[1] * pM[5] - pM[4] * pM[2];
91 const Double_t c11 = pM[8] * pM[0] - pM[2] * pM[2];
92 const Double_t c12 = pM[2] * pM[1] - pM[5] * pM[0];
93 const Double_t c22 = pM[0] * pM[4] - pM[1] * pM[1];
94
95 const Double_t t0 = TMath::Abs(pM[0]);
96 const Double_t t1 = TMath::Abs(pM[1]);
97 const Double_t t2 = TMath::Abs(pM[2]);
98
99 Double_t det;
100 Double_t tmp;
101
102 if (t0 >= t1) {
103 if (t2 >= t0) {
104 tmp = pM[2];
105 det = c12*c01-c11*c02;
106 } else {
107 tmp = pM[0];
108 det = c11*c22-c12*c12;
109 }
110 } else if (t2 >= t1) {
111 tmp = pM[2];
112 det = c12*c01-c11*c02;
113 } else {
114 tmp = pM[1];
115 det = c02*c12-c01*c22;
116 }
117
118 if ( det == 0 || tmp == 0) {
119 Error("Inv3x3","matrix is singular");
120 return kFALSE;
121 }
122
123 Double_t s = tmp/det;
124 if (determ)
125 *determ = 1./s;
126
127 pM[0] = s*c00;
128 pM[1] = s*c01;
129 pM[2] = s*c02;
130 pM[3] = pM[1];
131 pM[4] = s*c11;
132 pM[5] = s*c12;
133 pM[6] = pM[2];
134 pM[7] = pM[5];
135 pM[8] = s*c22;
136
137 return kTRUE;
138}
139
140// SFij are indices for a 4x4 symmetric matrix.
141
142#define SF00 0
143#define SF01 1
144#define SF02 2
145#define SF03 3
146
147#define SF10 1
148#define SF11 5
149#define SF12 6
150#define SF13 7
151
152#define SF20 2
153#define SF21 6
154#define SF22 10
155#define SF23 11
156
157#define SF30 3
158#define SF31 7
159#define SF32 11
160#define SF33 15
161
162////////////////////////////////////////////////////////////////////////////////
163
164template<class Element>
165Bool_t TMatrixTSymCramerInv::Inv4x4(TMatrixTSym<Element> &m,Double_t *determ)
166{
167 if (m.GetNrows() != 4) {
168 Error("Inv4x4","matrix should be square 4x4");
169 return kFALSE;
170 }
171
172 Element *pM = m.GetMatrixArray();
173
174 // Find all NECESSARY 2x2 dets: (14 of them)
175
176 const Double_t mDet2_12_01 = pM[SF10]*pM[SF21] - pM[SF11]*pM[SF20];
177 const Double_t mDet2_12_02 = pM[SF10]*pM[SF22] - pM[SF12]*pM[SF20];
178 const Double_t mDet2_12_12 = pM[SF11]*pM[SF22] - pM[SF12]*pM[SF21];
179 const Double_t mDet2_13_01 = pM[SF10]*pM[SF31] - pM[SF11]*pM[SF30];
180 const Double_t mDet2_13_02 = pM[SF10]*pM[SF32] - pM[SF12]*pM[SF30];
181 const Double_t mDet2_13_03 = pM[SF10]*pM[SF33] - pM[SF13]*pM[SF30];
182 const Double_t mDet2_13_12 = pM[SF11]*pM[SF32] - pM[SF12]*pM[SF31];
183 const Double_t mDet2_13_13 = pM[SF11]*pM[SF33] - pM[SF13]*pM[SF31];
184 const Double_t mDet2_23_01 = pM[SF20]*pM[SF31] - pM[SF21]*pM[SF30];
185 const Double_t mDet2_23_02 = pM[SF20]*pM[SF32] - pM[SF22]*pM[SF30];
186 const Double_t mDet2_23_03 = pM[SF20]*pM[SF33] - pM[SF23]*pM[SF30];
187 const Double_t mDet2_23_12 = pM[SF21]*pM[SF32] - pM[SF22]*pM[SF31];
188 const Double_t mDet2_23_13 = pM[SF21]*pM[SF33] - pM[SF23]*pM[SF31];
189 const Double_t mDet2_23_23 = pM[SF22]*pM[SF33] - pM[SF23]*pM[SF32];
190
191 // SFind all NECESSSFRY 3x3 dets: (10 of them)
192
193 const Double_t mDet3_012_012 = pM[SF00]*mDet2_12_12 - pM[SF01]*mDet2_12_02
194 + pM[SF02]*mDet2_12_01;
195 const Double_t mDet3_013_012 = pM[SF00]*mDet2_13_12 - pM[SF01]*mDet2_13_02
196 + pM[SF02]*mDet2_13_01;
197 const Double_t mDet3_013_013 = pM[SF00]*mDet2_13_13 - pM[SF01]*mDet2_13_03
198 + pM[SF03]*mDet2_13_01;
199 const Double_t mDet3_023_012 = pM[SF00]*mDet2_23_12 - pM[SF01]*mDet2_23_02
200 + pM[SF02]*mDet2_23_01;
201 const Double_t mDet3_023_013 = pM[SF00]*mDet2_23_13 - pM[SF01]*mDet2_23_03
202 + pM[SF03]*mDet2_23_01;
203 const Double_t mDet3_023_023 = pM[SF00]*mDet2_23_23 - pM[SF02]*mDet2_23_03
204 + pM[SF03]*mDet2_23_02;
205 const Double_t mDet3_123_012 = pM[SF10]*mDet2_23_12 - pM[SF11]*mDet2_23_02
206 + pM[SF12]*mDet2_23_01;
207 const Double_t mDet3_123_013 = pM[SF10]*mDet2_23_13 - pM[SF11]*mDet2_23_03
208 + pM[SF13]*mDet2_23_01;
209 const Double_t mDet3_123_023 = pM[SF10]*mDet2_23_23 - pM[SF12]*mDet2_23_03
210 + pM[SF13]*mDet2_23_02;
211 const Double_t mDet3_123_123 = pM[SF11]*mDet2_23_23 - pM[SF12]*mDet2_23_13
212 + pM[SF13]*mDet2_23_12;
213
214 // Find the 4x4 det:
215
216 const Double_t det = pM[SF00]*mDet3_123_123 - pM[SF01]*mDet3_123_023
217 + pM[SF02]*mDet3_123_013 - pM[SF03]*mDet3_123_012;
218
219 if (determ)
220 *determ = det;
221
222 if ( det == 0 ) {
223 Error("Inv4x4","matrix is singular");
224 return kFALSE;
225 }
226
227 const Double_t oneOverDet = 1.0/det;
228 const Double_t mn1OverDet = - oneOverDet;
229
230 pM[SF00] = mDet3_123_123 * oneOverDet;
231 pM[SF01] = mDet3_123_023 * mn1OverDet;
232 pM[SF02] = mDet3_123_013 * oneOverDet;
233 pM[SF03] = mDet3_123_012 * mn1OverDet;
234
235 pM[SF11] = mDet3_023_023 * oneOverDet;
236 pM[SF12] = mDet3_023_013 * mn1OverDet;
237 pM[SF13] = mDet3_023_012 * oneOverDet;
238
239 pM[SF22] = mDet3_013_013 * oneOverDet;
240 pM[SF23] = mDet3_013_012 * mn1OverDet;
241
242 pM[SF33] = mDet3_012_012 * oneOverDet;
243
244 for (Int_t irow = 0; irow < 4; irow++) {
245 const Int_t rowOff1 = irow*4;
246 for (Int_t icol = 0; icol < irow; icol++) {
247 const Int_t rowOff2 = icol*4;
248 pM[rowOff1+icol] = pM[rowOff2+irow];
249 }
250 }
251
252 return kTRUE;
253}
254
255// Mij are indices for a 5x5 matrix.
256
257#define SM00 0
258#define SM01 1
259#define SM02 2
260#define SM03 3
261#define SM04 4
262
263#define SM10 1
264#define SM11 6
265#define SM12 7
266#define SM13 8
267#define SM14 9
268
269#define SM20 2
270#define SM21 7
271#define SM22 12
272#define SM23 13
273#define SM24 14
274
275#define SM30 3
276#define SM31 8
277#define SM32 13
278#define SM33 18
279#define SM34 19
280
281#define SM40 4
282#define SM41 9
283#define SM42 14
284#define SM43 19
285#define SM44 24
286
287////////////////////////////////////////////////////////////////////////////////
288
289template<class Element>
290Bool_t TMatrixTSymCramerInv::Inv5x5(TMatrixTSym<Element> &m,Double_t *determ)
291{
292 if (m.GetNrows() != 5) {
293 Error("Inv5x5","matrix should be square 5x5");
294 return kFALSE;
295 }
296
297 Element *pM = m.GetMatrixArray();
298
299 // Find all NECESSARY 2x2 dets: (25 of them)
300
301 const Double_t mDet2_23_01 = pM[SM20]*pM[SM31] - pM[SM21]*pM[SM30];
302 const Double_t mDet2_23_02 = pM[SM20]*pM[SM32] - pM[SM22]*pM[SM30];
303 const Double_t mDet2_23_03 = pM[SM20]*pM[SM33] - pM[SM23]*pM[SM30];
304 const Double_t mDet2_23_12 = pM[SM21]*pM[SM32] - pM[SM22]*pM[SM31];
305 const Double_t mDet2_23_13 = pM[SM21]*pM[SM33] - pM[SM23]*pM[SM31];
306 const Double_t mDet2_23_23 = pM[SM22]*pM[SM33] - pM[SM23]*pM[SM32];
307 const Double_t mDet2_24_01 = pM[SM20]*pM[SM41] - pM[SM21]*pM[SM40];
308 const Double_t mDet2_24_02 = pM[SM20]*pM[SM42] - pM[SM22]*pM[SM40];
309 const Double_t mDet2_24_03 = pM[SM20]*pM[SM43] - pM[SM23]*pM[SM40];
310 const Double_t mDet2_24_04 = pM[SM20]*pM[SM44] - pM[SM24]*pM[SM40];
311 const Double_t mDet2_24_12 = pM[SM21]*pM[SM42] - pM[SM22]*pM[SM41];
312 const Double_t mDet2_24_13 = pM[SM21]*pM[SM43] - pM[SM23]*pM[SM41];
313 const Double_t mDet2_24_14 = pM[SM21]*pM[SM44] - pM[SM24]*pM[SM41];
314 const Double_t mDet2_24_23 = pM[SM22]*pM[SM43] - pM[SM23]*pM[SM42];
315 const Double_t mDet2_24_24 = pM[SM22]*pM[SM44] - pM[SM24]*pM[SM42];
316 const Double_t mDet2_34_01 = pM[SM30]*pM[SM41] - pM[SM31]*pM[SM40];
317 const Double_t mDet2_34_02 = pM[SM30]*pM[SM42] - pM[SM32]*pM[SM40];
318 const Double_t mDet2_34_03 = pM[SM30]*pM[SM43] - pM[SM33]*pM[SM40];
319 const Double_t mDet2_34_04 = pM[SM30]*pM[SM44] - pM[SM34]*pM[SM40];
320 const Double_t mDet2_34_12 = pM[SM31]*pM[SM42] - pM[SM32]*pM[SM41];
321 const Double_t mDet2_34_13 = pM[SM31]*pM[SM43] - pM[SM33]*pM[SM41];
322 const Double_t mDet2_34_14 = pM[SM31]*pM[SM44] - pM[SM34]*pM[SM41];
323 const Double_t mDet2_34_23 = pM[SM32]*pM[SM43] - pM[SM33]*pM[SM42];
324 const Double_t mDet2_34_24 = pM[SM32]*pM[SM44] - pM[SM34]*pM[SM42];
325 const Double_t mDet2_34_34 = pM[SM33]*pM[SM44] - pM[SM34]*pM[SM43];
326
327 // Find all NECESSARY 3x3 dets: (30 of them)
328
329 const Double_t mDet3_123_012 = pM[SM10]*mDet2_23_12 - pM[SM11]*mDet2_23_02 + pM[SM12]*mDet2_23_01;
330 const Double_t mDet3_123_013 = pM[SM10]*mDet2_23_13 - pM[SM11]*mDet2_23_03 + pM[SM13]*mDet2_23_01;
331 const Double_t mDet3_123_023 = pM[SM10]*mDet2_23_23 - pM[SM12]*mDet2_23_03 + pM[SM13]*mDet2_23_02;
332 const Double_t mDet3_123_123 = pM[SM11]*mDet2_23_23 - pM[SM12]*mDet2_23_13 + pM[SM13]*mDet2_23_12;
333 const Double_t mDet3_124_012 = pM[SM10]*mDet2_24_12 - pM[SM11]*mDet2_24_02 + pM[SM12]*mDet2_24_01;
334 const Double_t mDet3_124_013 = pM[SM10]*mDet2_24_13 - pM[SM11]*mDet2_24_03 + pM[SM13]*mDet2_24_01;
335 const Double_t mDet3_124_014 = pM[SM10]*mDet2_24_14 - pM[SM11]*mDet2_24_04 + pM[SM14]*mDet2_24_01;
336 const Double_t mDet3_124_023 = pM[SM10]*mDet2_24_23 - pM[SM12]*mDet2_24_03 + pM[SM13]*mDet2_24_02;
337 const Double_t mDet3_124_024 = pM[SM10]*mDet2_24_24 - pM[SM12]*mDet2_24_04 + pM[SM14]*mDet2_24_02;
338 const Double_t mDet3_124_123 = pM[SM11]*mDet2_24_23 - pM[SM12]*mDet2_24_13 + pM[SM13]*mDet2_24_12;
339 const Double_t mDet3_124_124 = pM[SM11]*mDet2_24_24 - pM[SM12]*mDet2_24_14 + pM[SM14]*mDet2_24_12;
340 const Double_t mDet3_134_012 = pM[SM10]*mDet2_34_12 - pM[SM11]*mDet2_34_02 + pM[SM12]*mDet2_34_01;
341 const Double_t mDet3_134_013 = pM[SM10]*mDet2_34_13 - pM[SM11]*mDet2_34_03 + pM[SM13]*mDet2_34_01;
342 const Double_t mDet3_134_014 = pM[SM10]*mDet2_34_14 - pM[SM11]*mDet2_34_04 + pM[SM14]*mDet2_34_01;
343 const Double_t mDet3_134_023 = pM[SM10]*mDet2_34_23 - pM[SM12]*mDet2_34_03 + pM[SM13]*mDet2_34_02;
344 const Double_t mDet3_134_024 = pM[SM10]*mDet2_34_24 - pM[SM12]*mDet2_34_04 + pM[SM14]*mDet2_34_02;
345 const Double_t mDet3_134_034 = pM[SM10]*mDet2_34_34 - pM[SM13]*mDet2_34_04 + pM[SM14]*mDet2_34_03;
346 const Double_t mDet3_134_123 = pM[SM11]*mDet2_34_23 - pM[SM12]*mDet2_34_13 + pM[SM13]*mDet2_34_12;
347 const Double_t mDet3_134_124 = pM[SM11]*mDet2_34_24 - pM[SM12]*mDet2_34_14 + pM[SM14]*mDet2_34_12;
348 const Double_t mDet3_134_134 = pM[SM11]*mDet2_34_34 - pM[SM13]*mDet2_34_14 + pM[SM14]*mDet2_34_13;
349 const Double_t mDet3_234_012 = pM[SM20]*mDet2_34_12 - pM[SM21]*mDet2_34_02 + pM[SM22]*mDet2_34_01;
350 const Double_t mDet3_234_013 = pM[SM20]*mDet2_34_13 - pM[SM21]*mDet2_34_03 + pM[SM23]*mDet2_34_01;
351 const Double_t mDet3_234_014 = pM[SM20]*mDet2_34_14 - pM[SM21]*mDet2_34_04 + pM[SM24]*mDet2_34_01;
352 const Double_t mDet3_234_023 = pM[SM20]*mDet2_34_23 - pM[SM22]*mDet2_34_03 + pM[SM23]*mDet2_34_02;
353 const Double_t mDet3_234_024 = pM[SM20]*mDet2_34_24 - pM[SM22]*mDet2_34_04 + pM[SM24]*mDet2_34_02;
354 const Double_t mDet3_234_034 = pM[SM20]*mDet2_34_34 - pM[SM23]*mDet2_34_04 + pM[SM24]*mDet2_34_03;
355 const Double_t mDet3_234_123 = pM[SM21]*mDet2_34_23 - pM[SM22]*mDet2_34_13 + pM[SM23]*mDet2_34_12;
356 const Double_t mDet3_234_124 = pM[SM21]*mDet2_34_24 - pM[SM22]*mDet2_34_14 + pM[SM24]*mDet2_34_12;
357 const Double_t mDet3_234_134 = pM[SM21]*mDet2_34_34 - pM[SM23]*mDet2_34_14 + pM[SM24]*mDet2_34_13;
358 const Double_t mDet3_234_234 = pM[SM22]*mDet2_34_34 - pM[SM23]*mDet2_34_24 + pM[SM24]*mDet2_34_23;
359
360 // Find all NECESSARY 4x4 dets: (15 of them)
361
362 const Double_t mDet4_0123_0123 = pM[SM00]*mDet3_123_123 - pM[SM01]*mDet3_123_023
363 + pM[SM02]*mDet3_123_013 - pM[SM03]*mDet3_123_012;
364 const Double_t mDet4_0124_0123 = pM[SM00]*mDet3_124_123 - pM[SM01]*mDet3_124_023
365 + pM[SM02]*mDet3_124_013 - pM[SM03]*mDet3_124_012;
366 const Double_t mDet4_0124_0124 = pM[SM00]*mDet3_124_124 - pM[SM01]*mDet3_124_024
367 + pM[SM02]*mDet3_124_014 - pM[SM04]*mDet3_124_012;
368 const Double_t mDet4_0134_0123 = pM[SM00]*mDet3_134_123 - pM[SM01]*mDet3_134_023
369 + pM[SM02]*mDet3_134_013 - pM[SM03]*mDet3_134_012;
370 const Double_t mDet4_0134_0124 = pM[SM00]*mDet3_134_124 - pM[SM01]*mDet3_134_024
371 + pM[SM02]*mDet3_134_014 - pM[SM04]*mDet3_134_012;
372 const Double_t mDet4_0134_0134 = pM[SM00]*mDet3_134_134 - pM[SM01]*mDet3_134_034
373 + pM[SM03]*mDet3_134_014 - pM[SM04]*mDet3_134_013;
374 const Double_t mDet4_0234_0123 = pM[SM00]*mDet3_234_123 - pM[SM01]*mDet3_234_023
375 + pM[SM02]*mDet3_234_013 - pM[SM03]*mDet3_234_012;
376 const Double_t mDet4_0234_0124 = pM[SM00]*mDet3_234_124 - pM[SM01]*mDet3_234_024
377 + pM[SM02]*mDet3_234_014 - pM[SM04]*mDet3_234_012;
378 const Double_t mDet4_0234_0134 = pM[SM00]*mDet3_234_134 - pM[SM01]*mDet3_234_034
379 + pM[SM03]*mDet3_234_014 - pM[SM04]*mDet3_234_013;
380 const Double_t mDet4_0234_0234 = pM[SM00]*mDet3_234_234 - pM[SM02]*mDet3_234_034
381 + pM[SM03]*mDet3_234_024 - pM[SM04]*mDet3_234_023;
382 const Double_t mDet4_1234_0123 = pM[SM10]*mDet3_234_123 - pM[SM11]*mDet3_234_023
383 + pM[SM12]*mDet3_234_013 - pM[SM13]*mDet3_234_012;
384 const Double_t mDet4_1234_0124 = pM[SM10]*mDet3_234_124 - pM[SM11]*mDet3_234_024
385 + pM[SM12]*mDet3_234_014 - pM[SM14]*mDet3_234_012;
386 const Double_t mDet4_1234_0134 = pM[SM10]*mDet3_234_134 - pM[SM11]*mDet3_234_034
387 + pM[SM13]*mDet3_234_014 - pM[SM14]*mDet3_234_013;
388 const Double_t mDet4_1234_0234 = pM[SM10]*mDet3_234_234 - pM[SM12]*mDet3_234_034
389 + pM[SM13]*mDet3_234_024 - pM[SM14]*mDet3_234_023;
390 const Double_t mDet4_1234_1234 = pM[SM11]*mDet3_234_234 - pM[SM12]*mDet3_234_134
391 + pM[SM13]*mDet3_234_124 - pM[SM14]*mDet3_234_123;
392
393 // Find the 5x5 det:
394
395 const Double_t det = pM[SM00]*mDet4_1234_1234 - pM[SM01]*mDet4_1234_0234 + pM[SM02]*mDet4_1234_0134
396 - pM[SM03]*mDet4_1234_0124 + pM[SM04]*mDet4_1234_0123;
397 if (determ)
398 *determ = det;
399
400 if ( det == 0 ) {
401 Error("Inv5x5","matrix is singular");
402 return kFALSE;
403 }
404
405 const Double_t oneOverDet = 1.0/det;
406 const Double_t mn1OverDet = - oneOverDet;
407
408 pM[SM00] = mDet4_1234_1234 * oneOverDet;
409 pM[SM01] = mDet4_1234_0234 * mn1OverDet;
410 pM[SM02] = mDet4_1234_0134 * oneOverDet;
411 pM[SM03] = mDet4_1234_0124 * mn1OverDet;
412 pM[SM04] = mDet4_1234_0123 * oneOverDet;
413
414 pM[SM11] = mDet4_0234_0234 * oneOverDet;
415 pM[SM12] = mDet4_0234_0134 * mn1OverDet;
416 pM[SM13] = mDet4_0234_0124 * oneOverDet;
417 pM[SM14] = mDet4_0234_0123 * mn1OverDet;
418
419 pM[SM22] = mDet4_0134_0134 * oneOverDet;
420 pM[SM23] = mDet4_0134_0124 * mn1OverDet;
421 pM[SM24] = mDet4_0134_0123 * oneOverDet;
422
423 pM[SM33] = mDet4_0124_0124 * oneOverDet;
424 pM[SM34] = mDet4_0124_0123 * mn1OverDet;
425
426 pM[SM44] = mDet4_0123_0123 * oneOverDet;
427
428 for (Int_t irow = 0; irow < 5; irow++) {
429 const Int_t rowOff1 = irow*5;
430 for (Int_t icol = 0; icol < irow; icol++) {
431 const Int_t rowOff2 = icol*5;
432 pM[rowOff1+icol] = pM[rowOff2+irow];
433 }
434 }
435
436 return kTRUE;
437}
438
439// Aij are indices for a 6x6 symmetric matrix.
440
441#define SA00 0
442#define SA01 1
443#define SA02 2
444#define SA03 3
445#define SA04 4
446#define SA05 5
447
448#define SA10 1
449#define SA11 7
450#define SA12 8
451#define SA13 9
452#define SA14 10
453#define SA15 11
454
455#define SA20 2
456#define SA21 8
457#define SA22 14
458#define SA23 15
459#define SA24 16
460#define SA25 17
461
462#define SA30 3
463#define SA31 9
464#define SA32 15
465#define SA33 21
466#define SA34 22
467#define SA35 23
468
469#define SA40 4
470#define SA41 10
471#define SA42 16
472#define SA43 22
473#define SA44 28
474#define SA45 29
475
476#define SA50 5
477#define SA51 11
478#define SA52 17
479#define SA53 23
480#define SA54 29
481#define SA55 35
482
483////////////////////////////////////////////////////////////////////////////////
484
485template<class Element>
486Bool_t TMatrixTSymCramerInv::Inv6x6(TMatrixTSym<Element> &m,Double_t *determ)
487{
488 if (m.GetNrows() != 6 || m.GetNcols() != 6 || m.GetRowLwb() != m.GetColLwb()) {
489 Error("Inv6x6","matrix should be square 6x6");
490 return kFALSE;
491 }
492
493 Element *pM = m.GetMatrixArray();
494
495 // Find all NECESSSARY 2x2 dets: (39 of them)
496
497 const Double_t mDet2_34_01 = pM[SA30]*pM[SA41] - pM[SA31]*pM[SA40];
498 const Double_t mDet2_34_02 = pM[SA30]*pM[SA42] - pM[SA32]*pM[SA40];
499 const Double_t mDet2_34_03 = pM[SA30]*pM[SA43] - pM[SA33]*pM[SA40];
500 const Double_t mDet2_34_04 = pM[SA30]*pM[SA44] - pM[SA34]*pM[SA40];
501 const Double_t mDet2_34_12 = pM[SA31]*pM[SA42] - pM[SA32]*pM[SA41];
502 const Double_t mDet2_34_13 = pM[SA31]*pM[SA43] - pM[SA33]*pM[SA41];
503 const Double_t mDet2_34_14 = pM[SA31]*pM[SA44] - pM[SA34]*pM[SA41];
504 const Double_t mDet2_34_23 = pM[SA32]*pM[SA43] - pM[SA33]*pM[SA42];
505 const Double_t mDet2_34_24 = pM[SA32]*pM[SA44] - pM[SA34]*pM[SA42];
506 const Double_t mDet2_34_34 = pM[SA33]*pM[SA44] - pM[SA34]*pM[SA43];
507 const Double_t mDet2_35_01 = pM[SA30]*pM[SA51] - pM[SA31]*pM[SA50];
508 const Double_t mDet2_35_02 = pM[SA30]*pM[SA52] - pM[SA32]*pM[SA50];
509 const Double_t mDet2_35_03 = pM[SA30]*pM[SA53] - pM[SA33]*pM[SA50];
510 const Double_t mDet2_35_04 = pM[SA30]*pM[SA54] - pM[SA34]*pM[SA50];
511 const Double_t mDet2_35_05 = pM[SA30]*pM[SA55] - pM[SA35]*pM[SA50];
512 const Double_t mDet2_35_12 = pM[SA31]*pM[SA52] - pM[SA32]*pM[SA51];
513 const Double_t mDet2_35_13 = pM[SA31]*pM[SA53] - pM[SA33]*pM[SA51];
514 const Double_t mDet2_35_14 = pM[SA31]*pM[SA54] - pM[SA34]*pM[SA51];
515 const Double_t mDet2_35_15 = pM[SA31]*pM[SA55] - pM[SA35]*pM[SA51];
516 const Double_t mDet2_35_23 = pM[SA32]*pM[SA53] - pM[SA33]*pM[SA52];
517 const Double_t mDet2_35_24 = pM[SA32]*pM[SA54] - pM[SA34]*pM[SA52];
518 const Double_t mDet2_35_25 = pM[SA32]*pM[SA55] - pM[SA35]*pM[SA52];
519 const Double_t mDet2_35_34 = pM[SA33]*pM[SA54] - pM[SA34]*pM[SA53];
520 const Double_t mDet2_35_35 = pM[SA33]*pM[SA55] - pM[SA35]*pM[SA53];
521 const Double_t mDet2_45_01 = pM[SA40]*pM[SA51] - pM[SA41]*pM[SA50];
522 const Double_t mDet2_45_02 = pM[SA40]*pM[SA52] - pM[SA42]*pM[SA50];
523 const Double_t mDet2_45_03 = pM[SA40]*pM[SA53] - pM[SA43]*pM[SA50];
524 const Double_t mDet2_45_04 = pM[SA40]*pM[SA54] - pM[SA44]*pM[SA50];
525 const Double_t mDet2_45_05 = pM[SA40]*pM[SA55] - pM[SA45]*pM[SA50];
526 const Double_t mDet2_45_12 = pM[SA41]*pM[SA52] - pM[SA42]*pM[SA51];
527 const Double_t mDet2_45_13 = pM[SA41]*pM[SA53] - pM[SA43]*pM[SA51];
528 const Double_t mDet2_45_14 = pM[SA41]*pM[SA54] - pM[SA44]*pM[SA51];
529 const Double_t mDet2_45_15 = pM[SA41]*pM[SA55] - pM[SA45]*pM[SA51];
530 const Double_t mDet2_45_23 = pM[SA42]*pM[SA53] - pM[SA43]*pM[SA52];
531 const Double_t mDet2_45_24 = pM[SA42]*pM[SA54] - pM[SA44]*pM[SA52];
532 const Double_t mDet2_45_25 = pM[SA42]*pM[SA55] - pM[SA45]*pM[SA52];
533 const Double_t mDet2_45_34 = pM[SA43]*pM[SA54] - pM[SA44]*pM[SA53];
534 const Double_t mDet2_45_35 = pM[SA43]*pM[SA55] - pM[SA45]*pM[SA53];
535 const Double_t mDet2_45_45 = pM[SA44]*pM[SA55] - pM[SA45]*pM[SA54];
536
537 // Find all NECESSSARY 3x3 dets: (65 of them)
538
539 const Double_t mDet3_234_012 = pM[SA20]*mDet2_34_12 - pM[SA21]*mDet2_34_02 + pM[SA22]*mDet2_34_01;
540 const Double_t mDet3_234_013 = pM[SA20]*mDet2_34_13 - pM[SA21]*mDet2_34_03 + pM[SA23]*mDet2_34_01;
541 const Double_t mDet3_234_014 = pM[SA20]*mDet2_34_14 - pM[SA21]*mDet2_34_04 + pM[SA24]*mDet2_34_01;
542 const Double_t mDet3_234_023 = pM[SA20]*mDet2_34_23 - pM[SA22]*mDet2_34_03 + pM[SA23]*mDet2_34_02;
543 const Double_t mDet3_234_024 = pM[SA20]*mDet2_34_24 - pM[SA22]*mDet2_34_04 + pM[SA24]*mDet2_34_02;
544 const Double_t mDet3_234_034 = pM[SA20]*mDet2_34_34 - pM[SA23]*mDet2_34_04 + pM[SA24]*mDet2_34_03;
545 const Double_t mDet3_234_123 = pM[SA21]*mDet2_34_23 - pM[SA22]*mDet2_34_13 + pM[SA23]*mDet2_34_12;
546 const Double_t mDet3_234_124 = pM[SA21]*mDet2_34_24 - pM[SA22]*mDet2_34_14 + pM[SA24]*mDet2_34_12;
547 const Double_t mDet3_234_134 = pM[SA21]*mDet2_34_34 - pM[SA23]*mDet2_34_14 + pM[SA24]*mDet2_34_13;
548 const Double_t mDet3_234_234 = pM[SA22]*mDet2_34_34 - pM[SA23]*mDet2_34_24 + pM[SA24]*mDet2_34_23;
549 const Double_t mDet3_235_012 = pM[SA20]*mDet2_35_12 - pM[SA21]*mDet2_35_02 + pM[SA22]*mDet2_35_01;
550 const Double_t mDet3_235_013 = pM[SA20]*mDet2_35_13 - pM[SA21]*mDet2_35_03 + pM[SA23]*mDet2_35_01;
551 const Double_t mDet3_235_014 = pM[SA20]*mDet2_35_14 - pM[SA21]*mDet2_35_04 + pM[SA24]*mDet2_35_01;
552 const Double_t mDet3_235_015 = pM[SA20]*mDet2_35_15 - pM[SA21]*mDet2_35_05 + pM[SA25]*mDet2_35_01;
553 const Double_t mDet3_235_023 = pM[SA20]*mDet2_35_23 - pM[SA22]*mDet2_35_03 + pM[SA23]*mDet2_35_02;
554 const Double_t mDet3_235_024 = pM[SA20]*mDet2_35_24 - pM[SA22]*mDet2_35_04 + pM[SA24]*mDet2_35_02;
555 const Double_t mDet3_235_025 = pM[SA20]*mDet2_35_25 - pM[SA22]*mDet2_35_05 + pM[SA25]*mDet2_35_02;
556 const Double_t mDet3_235_034 = pM[SA20]*mDet2_35_34 - pM[SA23]*mDet2_35_04 + pM[SA24]*mDet2_35_03;
557 const Double_t mDet3_235_035 = pM[SA20]*mDet2_35_35 - pM[SA23]*mDet2_35_05 + pM[SA25]*mDet2_35_03;
558 const Double_t mDet3_235_123 = pM[SA21]*mDet2_35_23 - pM[SA22]*mDet2_35_13 + pM[SA23]*mDet2_35_12;
559 const Double_t mDet3_235_124 = pM[SA21]*mDet2_35_24 - pM[SA22]*mDet2_35_14 + pM[SA24]*mDet2_35_12;
560 const Double_t mDet3_235_125 = pM[SA21]*mDet2_35_25 - pM[SA22]*mDet2_35_15 + pM[SA25]*mDet2_35_12;
561 const Double_t mDet3_235_134 = pM[SA21]*mDet2_35_34 - pM[SA23]*mDet2_35_14 + pM[SA24]*mDet2_35_13;
562 const Double_t mDet3_235_135 = pM[SA21]*mDet2_35_35 - pM[SA23]*mDet2_35_15 + pM[SA25]*mDet2_35_13;
563 const Double_t mDet3_235_234 = pM[SA22]*mDet2_35_34 - pM[SA23]*mDet2_35_24 + pM[SA24]*mDet2_35_23;
564 const Double_t mDet3_235_235 = pM[SA22]*mDet2_35_35 - pM[SA23]*mDet2_35_25 + pM[SA25]*mDet2_35_23;
565 const Double_t mDet3_245_012 = pM[SA20]*mDet2_45_12 - pM[SA21]*mDet2_45_02 + pM[SA22]*mDet2_45_01;
566 const Double_t mDet3_245_013 = pM[SA20]*mDet2_45_13 - pM[SA21]*mDet2_45_03 + pM[SA23]*mDet2_45_01;
567 const Double_t mDet3_245_014 = pM[SA20]*mDet2_45_14 - pM[SA21]*mDet2_45_04 + pM[SA24]*mDet2_45_01;
568 const Double_t mDet3_245_015 = pM[SA20]*mDet2_45_15 - pM[SA21]*mDet2_45_05 + pM[SA25]*mDet2_45_01;
569 const Double_t mDet3_245_023 = pM[SA20]*mDet2_45_23 - pM[SA22]*mDet2_45_03 + pM[SA23]*mDet2_45_02;
570 const Double_t mDet3_245_024 = pM[SA20]*mDet2_45_24 - pM[SA22]*mDet2_45_04 + pM[SA24]*mDet2_45_02;
571 const Double_t mDet3_245_025 = pM[SA20]*mDet2_45_25 - pM[SA22]*mDet2_45_05 + pM[SA25]*mDet2_45_02;
572 const Double_t mDet3_245_034 = pM[SA20]*mDet2_45_34 - pM[SA23]*mDet2_45_04 + pM[SA24]*mDet2_45_03;
573 const Double_t mDet3_245_035 = pM[SA20]*mDet2_45_35 - pM[SA23]*mDet2_45_05 + pM[SA25]*mDet2_45_03;
574 const Double_t mDet3_245_045 = pM[SA20]*mDet2_45_45 - pM[SA24]*mDet2_45_05 + pM[SA25]*mDet2_45_04;
575 const Double_t mDet3_245_123 = pM[SA21]*mDet2_45_23 - pM[SA22]*mDet2_45_13 + pM[SA23]*mDet2_45_12;
576 const Double_t mDet3_245_124 = pM[SA21]*mDet2_45_24 - pM[SA22]*mDet2_45_14 + pM[SA24]*mDet2_45_12;
577 const Double_t mDet3_245_125 = pM[SA21]*mDet2_45_25 - pM[SA22]*mDet2_45_15 + pM[SA25]*mDet2_45_12;
578 const Double_t mDet3_245_134 = pM[SA21]*mDet2_45_34 - pM[SA23]*mDet2_45_14 + pM[SA24]*mDet2_45_13;
579 const Double_t mDet3_245_135 = pM[SA21]*mDet2_45_35 - pM[SA23]*mDet2_45_15 + pM[SA25]*mDet2_45_13;
580 const Double_t mDet3_245_145 = pM[SA21]*mDet2_45_45 - pM[SA24]*mDet2_45_15 + pM[SA25]*mDet2_45_14;
581 const Double_t mDet3_245_234 = pM[SA22]*mDet2_45_34 - pM[SA23]*mDet2_45_24 + pM[SA24]*mDet2_45_23;
582 const Double_t mDet3_245_235 = pM[SA22]*mDet2_45_35 - pM[SA23]*mDet2_45_25 + pM[SA25]*mDet2_45_23;
583 const Double_t mDet3_245_245 = pM[SA22]*mDet2_45_45 - pM[SA24]*mDet2_45_25 + pM[SA25]*mDet2_45_24;
584 const Double_t mDet3_345_012 = pM[SA30]*mDet2_45_12 - pM[SA31]*mDet2_45_02 + pM[SA32]*mDet2_45_01;
585 const Double_t mDet3_345_013 = pM[SA30]*mDet2_45_13 - pM[SA31]*mDet2_45_03 + pM[SA33]*mDet2_45_01;
586 const Double_t mDet3_345_014 = pM[SA30]*mDet2_45_14 - pM[SA31]*mDet2_45_04 + pM[SA34]*mDet2_45_01;
587 const Double_t mDet3_345_015 = pM[SA30]*mDet2_45_15 - pM[SA31]*mDet2_45_05 + pM[SA35]*mDet2_45_01;
588 const Double_t mDet3_345_023 = pM[SA30]*mDet2_45_23 - pM[SA32]*mDet2_45_03 + pM[SA33]*mDet2_45_02;
589 const Double_t mDet3_345_024 = pM[SA30]*mDet2_45_24 - pM[SA32]*mDet2_45_04 + pM[SA34]*mDet2_45_02;
590 const Double_t mDet3_345_025 = pM[SA30]*mDet2_45_25 - pM[SA32]*mDet2_45_05 + pM[SA35]*mDet2_45_02;
591 const Double_t mDet3_345_034 = pM[SA30]*mDet2_45_34 - pM[SA33]*mDet2_45_04 + pM[SA34]*mDet2_45_03;
592 const Double_t mDet3_345_035 = pM[SA30]*mDet2_45_35 - pM[SA33]*mDet2_45_05 + pM[SA35]*mDet2_45_03;
593 const Double_t mDet3_345_045 = pM[SA30]*mDet2_45_45 - pM[SA34]*mDet2_45_05 + pM[SA35]*mDet2_45_04;
594 const Double_t mDet3_345_123 = pM[SA31]*mDet2_45_23 - pM[SA32]*mDet2_45_13 + pM[SA33]*mDet2_45_12;
595 const Double_t mDet3_345_124 = pM[SA31]*mDet2_45_24 - pM[SA32]*mDet2_45_14 + pM[SA34]*mDet2_45_12;
596 const Double_t mDet3_345_125 = pM[SA31]*mDet2_45_25 - pM[SA32]*mDet2_45_15 + pM[SA35]*mDet2_45_12;
597 const Double_t mDet3_345_134 = pM[SA31]*mDet2_45_34 - pM[SA33]*mDet2_45_14 + pM[SA34]*mDet2_45_13;
598 const Double_t mDet3_345_135 = pM[SA31]*mDet2_45_35 - pM[SA33]*mDet2_45_15 + pM[SA35]*mDet2_45_13;
599 const Double_t mDet3_345_145 = pM[SA31]*mDet2_45_45 - pM[SA34]*mDet2_45_15 + pM[SA35]*mDet2_45_14;
600 const Double_t mDet3_345_234 = pM[SA32]*mDet2_45_34 - pM[SA33]*mDet2_45_24 + pM[SA34]*mDet2_45_23;
601 const Double_t mDet3_345_235 = pM[SA32]*mDet2_45_35 - pM[SA33]*mDet2_45_25 + pM[SA35]*mDet2_45_23;
602 const Double_t mDet3_345_245 = pM[SA32]*mDet2_45_45 - pM[SA34]*mDet2_45_25 + pM[SA35]*mDet2_45_24;
603 const Double_t mDet3_345_345 = pM[SA33]*mDet2_45_45 - pM[SA34]*mDet2_45_35 + pM[SA35]*mDet2_45_34;
604
605 // Find all NECESSSARY 4x4 dets: (55 of them)
606
607 const Double_t mDet4_1234_0123 = pM[SA10]*mDet3_234_123 - pM[SA11]*mDet3_234_023
608 + pM[SA12]*mDet3_234_013 - pM[SA13]*mDet3_234_012;
609 const Double_t mDet4_1234_0124 = pM[SA10]*mDet3_234_124 - pM[SA11]*mDet3_234_024
610 + pM[SA12]*mDet3_234_014 - pM[SA14]*mDet3_234_012;
611 const Double_t mDet4_1234_0134 = pM[SA10]*mDet3_234_134 - pM[SA11]*mDet3_234_034
612 + pM[SA13]*mDet3_234_014 - pM[SA14]*mDet3_234_013;
613 const Double_t mDet4_1234_0234 = pM[SA10]*mDet3_234_234 - pM[SA12]*mDet3_234_034
614 + pM[SA13]*mDet3_234_024 - pM[SA14]*mDet3_234_023;
615 const Double_t mDet4_1234_1234 = pM[SA11]*mDet3_234_234 - pM[SA12]*mDet3_234_134
616 + pM[SA13]*mDet3_234_124 - pM[SA14]*mDet3_234_123;
617 const Double_t mDet4_1235_0123 = pM[SA10]*mDet3_235_123 - pM[SA11]*mDet3_235_023
618 + pM[SA12]*mDet3_235_013 - pM[SA13]*mDet3_235_012;
619 const Double_t mDet4_1235_0124 = pM[SA10]*mDet3_235_124 - pM[SA11]*mDet3_235_024
620 + pM[SA12]*mDet3_235_014 - pM[SA14]*mDet3_235_012;
621 const Double_t mDet4_1235_0125 = pM[SA10]*mDet3_235_125 - pM[SA11]*mDet3_235_025
622 + pM[SA12]*mDet3_235_015 - pM[SA15]*mDet3_235_012;
623 const Double_t mDet4_1235_0134 = pM[SA10]*mDet3_235_134 - pM[SA11]*mDet3_235_034
624 + pM[SA13]*mDet3_235_014 - pM[SA14]*mDet3_235_013;
625 const Double_t mDet4_1235_0135 = pM[SA10]*mDet3_235_135 - pM[SA11]*mDet3_235_035
626 + pM[SA13]*mDet3_235_015 - pM[SA15]*mDet3_235_013;
627 const Double_t mDet4_1235_0234 = pM[SA10]*mDet3_235_234 - pM[SA12]*mDet3_235_034
628 + pM[SA13]*mDet3_235_024 - pM[SA14]*mDet3_235_023;
629 const Double_t mDet4_1235_0235 = pM[SA10]*mDet3_235_235 - pM[SA12]*mDet3_235_035
630 + pM[SA13]*mDet3_235_025 - pM[SA15]*mDet3_235_023;
631 const Double_t mDet4_1235_1234 = pM[SA11]*mDet3_235_234 - pM[SA12]*mDet3_235_134
632 + pM[SA13]*mDet3_235_124 - pM[SA14]*mDet3_235_123;
633 const Double_t mDet4_1235_1235 = pM[SA11]*mDet3_235_235 - pM[SA12]*mDet3_235_135
634 + pM[SA13]*mDet3_235_125 - pM[SA15]*mDet3_235_123;
635 const Double_t mDet4_1245_0123 = pM[SA10]*mDet3_245_123 - pM[SA11]*mDet3_245_023
636 + pM[SA12]*mDet3_245_013 - pM[SA13]*mDet3_245_012;
637 const Double_t mDet4_1245_0124 = pM[SA10]*mDet3_245_124 - pM[SA11]*mDet3_245_024
638 + pM[SA12]*mDet3_245_014 - pM[SA14]*mDet3_245_012;
639 const Double_t mDet4_1245_0125 = pM[SA10]*mDet3_245_125 - pM[SA11]*mDet3_245_025
640 + pM[SA12]*mDet3_245_015 - pM[SA15]*mDet3_245_012;
641 const Double_t mDet4_1245_0134 = pM[SA10]*mDet3_245_134 - pM[SA11]*mDet3_245_034
642 + pM[SA13]*mDet3_245_014 - pM[SA14]*mDet3_245_013;
643 const Double_t mDet4_1245_0135 = pM[SA10]*mDet3_245_135 - pM[SA11]*mDet3_245_035
644 + pM[SA13]*mDet3_245_015 - pM[SA15]*mDet3_245_013;
645 const Double_t mDet4_1245_0145 = pM[SA10]*mDet3_245_145 - pM[SA11]*mDet3_245_045
646 + pM[SA14]*mDet3_245_015 - pM[SA15]*mDet3_245_014;
647 const Double_t mDet4_1245_0234 = pM[SA10]*mDet3_245_234 - pM[SA12]*mDet3_245_034
648 + pM[SA13]*mDet3_245_024 - pM[SA14]*mDet3_245_023;
649 const Double_t mDet4_1245_0235 = pM[SA10]*mDet3_245_235 - pM[SA12]*mDet3_245_035
650 + pM[SA13]*mDet3_245_025 - pM[SA15]*mDet3_245_023;
651 const Double_t mDet4_1245_0245 = pM[SA10]*mDet3_245_245 - pM[SA12]*mDet3_245_045
652 + pM[SA14]*mDet3_245_025 - pM[SA15]*mDet3_245_024;
653 const Double_t mDet4_1245_1234 = pM[SA11]*mDet3_245_234 - pM[SA12]*mDet3_245_134
654 + pM[SA13]*mDet3_245_124 - pM[SA14]*mDet3_245_123;
655 const Double_t mDet4_1245_1235 = pM[SA11]*mDet3_245_235 - pM[SA12]*mDet3_245_135
656 + pM[SA13]*mDet3_245_125 - pM[SA15]*mDet3_245_123;
657 const Double_t mDet4_1245_1245 = pM[SA11]*mDet3_245_245 - pM[SA12]*mDet3_245_145
658 + pM[SA14]*mDet3_245_125 - pM[SA15]*mDet3_245_124;
659 const Double_t mDet4_1345_0123 = pM[SA10]*mDet3_345_123 - pM[SA11]*mDet3_345_023
660 + pM[SA12]*mDet3_345_013 - pM[SA13]*mDet3_345_012;
661 const Double_t mDet4_1345_0124 = pM[SA10]*mDet3_345_124 - pM[SA11]*mDet3_345_024
662 + pM[SA12]*mDet3_345_014 - pM[SA14]*mDet3_345_012;
663 const Double_t mDet4_1345_0125 = pM[SA10]*mDet3_345_125 - pM[SA11]*mDet3_345_025
664 + pM[SA12]*mDet3_345_015 - pM[SA15]*mDet3_345_012;
665 const Double_t mDet4_1345_0134 = pM[SA10]*mDet3_345_134 - pM[SA11]*mDet3_345_034
666 + pM[SA13]*mDet3_345_014 - pM[SA14]*mDet3_345_013;
667 const Double_t mDet4_1345_0135 = pM[SA10]*mDet3_345_135 - pM[SA11]*mDet3_345_035
668 + pM[SA13]*mDet3_345_015 - pM[SA15]*mDet3_345_013;
669 const Double_t mDet4_1345_0145 = pM[SA10]*mDet3_345_145 - pM[SA11]*mDet3_345_045
670 + pM[SA14]*mDet3_345_015 - pM[SA15]*mDet3_345_014;
671 const Double_t mDet4_1345_0234 = pM[SA10]*mDet3_345_234 - pM[SA12]*mDet3_345_034
672 + pM[SA13]*mDet3_345_024 - pM[SA14]*mDet3_345_023;
673 const Double_t mDet4_1345_0235 = pM[SA10]*mDet3_345_235 - pM[SA12]*mDet3_345_035
674 + pM[SA13]*mDet3_345_025 - pM[SA15]*mDet3_345_023;
675 const Double_t mDet4_1345_0245 = pM[SA10]*mDet3_345_245 - pM[SA12]*mDet3_345_045
676 + pM[SA14]*mDet3_345_025 - pM[SA15]*mDet3_345_024;
677 const Double_t mDet4_1345_0345 = pM[SA10]*mDet3_345_345 - pM[SA13]*mDet3_345_045
678 + pM[SA14]*mDet3_345_035 - pM[SA15]*mDet3_345_034;
679 const Double_t mDet4_1345_1234 = pM[SA11]*mDet3_345_234 - pM[SA12]*mDet3_345_134
680 + pM[SA13]*mDet3_345_124 - pM[SA14]*mDet3_345_123;
681 const Double_t mDet4_1345_1235 = pM[SA11]*mDet3_345_235 - pM[SA12]*mDet3_345_135
682 + pM[SA13]*mDet3_345_125 - pM[SA15]*mDet3_345_123;
683 const Double_t mDet4_1345_1245 = pM[SA11]*mDet3_345_245 - pM[SA12]*mDet3_345_145
684 + pM[SA14]*mDet3_345_125 - pM[SA15]*mDet3_345_124;
685 const Double_t mDet4_1345_1345 = pM[SA11]*mDet3_345_345 - pM[SA13]*mDet3_345_145
686 + pM[SA14]*mDet3_345_135 - pM[SA15]*mDet3_345_134;
687 const Double_t mDet4_2345_0123 = pM[SA20]*mDet3_345_123 - pM[SA21]*mDet3_345_023
688 + pM[SA22]*mDet3_345_013 - pM[SA23]*mDet3_345_012;
689 const Double_t mDet4_2345_0124 = pM[SA20]*mDet3_345_124 - pM[SA21]*mDet3_345_024
690 + pM[SA22]*mDet3_345_014 - pM[SA24]*mDet3_345_012;
691 const Double_t mDet4_2345_0125 = pM[SA20]*mDet3_345_125 - pM[SA21]*mDet3_345_025
692 + pM[SA22]*mDet3_345_015 - pM[SA25]*mDet3_345_012;
693 const Double_t mDet4_2345_0134 = pM[SA20]*mDet3_345_134 - pM[SA21]*mDet3_345_034
694 + pM[SA23]*mDet3_345_014 - pM[SA24]*mDet3_345_013;
695 const Double_t mDet4_2345_0135 = pM[SA20]*mDet3_345_135 - pM[SA21]*mDet3_345_035
696 + pM[SA23]*mDet3_345_015 - pM[SA25]*mDet3_345_013;
697 const Double_t mDet4_2345_0145 = pM[SA20]*mDet3_345_145 - pM[SA21]*mDet3_345_045
698 + pM[SA24]*mDet3_345_015 - pM[SA25]*mDet3_345_014;
699 const Double_t mDet4_2345_0234 = pM[SA20]*mDet3_345_234 - pM[SA22]*mDet3_345_034
700 + pM[SA23]*mDet3_345_024 - pM[SA24]*mDet3_345_023;
701 const Double_t mDet4_2345_0235 = pM[SA20]*mDet3_345_235 - pM[SA22]*mDet3_345_035
702 + pM[SA23]*mDet3_345_025 - pM[SA25]*mDet3_345_023;
703 const Double_t mDet4_2345_0245 = pM[SA20]*mDet3_345_245 - pM[SA22]*mDet3_345_045
704 + pM[SA24]*mDet3_345_025 - pM[SA25]*mDet3_345_024;
705 const Double_t mDet4_2345_0345 = pM[SA20]*mDet3_345_345 - pM[SA23]*mDet3_345_045
706 + pM[SA24]*mDet3_345_035 - pM[SA25]*mDet3_345_034;
707 const Double_t mDet4_2345_1234 = pM[SA21]*mDet3_345_234 - pM[SA22]*mDet3_345_134
708 + pM[SA23]*mDet3_345_124 - pM[SA24]*mDet3_345_123;
709 const Double_t mDet4_2345_1235 = pM[SA21]*mDet3_345_235 - pM[SA22]*mDet3_345_135
710 + pM[SA23]*mDet3_345_125 - pM[SA25]*mDet3_345_123;
711 const Double_t mDet4_2345_1245 = pM[SA21]*mDet3_345_245 - pM[SA22]*mDet3_345_145
712 + pM[SA24]*mDet3_345_125 - pM[SA25]*mDet3_345_124;
713 const Double_t mDet4_2345_1345 = pM[SA21]*mDet3_345_345 - pM[SA23]*mDet3_345_145
714 + pM[SA24]*mDet3_345_135 - pM[SA25]*mDet3_345_134;
715 const Double_t mDet4_2345_2345 = pM[SA22]*mDet3_345_345 - pM[SA23]*mDet3_345_245
716 + pM[SA24]*mDet3_345_235 - pM[SA25]*mDet3_345_234;
717
718 // Find all NECESSSARY 5x5 dets: (19 of them)
719
720 const Double_t mDet5_01234_01234 = pM[SA00]*mDet4_1234_1234 - pM[SA01]*mDet4_1234_0234
721 + pM[SA02]*mDet4_1234_0134 - pM[SA03]*mDet4_1234_0124 + pM[SA04]*mDet4_1234_0123;
722 const Double_t mDet5_01235_01234 = pM[SA00]*mDet4_1235_1234 - pM[SA01]*mDet4_1235_0234
723 + pM[SA02]*mDet4_1235_0134 - pM[SA03]*mDet4_1235_0124 + pM[SA04]*mDet4_1235_0123;
724 const Double_t mDet5_01235_01235 = pM[SA00]*mDet4_1235_1235 - pM[SA01]*mDet4_1235_0235
725 + pM[SA02]*mDet4_1235_0135 - pM[SA03]*mDet4_1235_0125 + pM[SA05]*mDet4_1235_0123;
726 const Double_t mDet5_01245_01234 = pM[SA00]*mDet4_1245_1234 - pM[SA01]*mDet4_1245_0234
727 + pM[SA02]*mDet4_1245_0134 - pM[SA03]*mDet4_1245_0124 + pM[SA04]*mDet4_1245_0123;
728 const Double_t mDet5_01245_01235 = pM[SA00]*mDet4_1245_1235 - pM[SA01]*mDet4_1245_0235
729 + pM[SA02]*mDet4_1245_0135 - pM[SA03]*mDet4_1245_0125 + pM[SA05]*mDet4_1245_0123;
730 const Double_t mDet5_01245_01245 = pM[SA00]*mDet4_1245_1245 - pM[SA01]*mDet4_1245_0245
731 + pM[SA02]*mDet4_1245_0145 - pM[SA04]*mDet4_1245_0125 + pM[SA05]*mDet4_1245_0124;
732 const Double_t mDet5_01345_01234 = pM[SA00]*mDet4_1345_1234 - pM[SA01]*mDet4_1345_0234
733 + pM[SA02]*mDet4_1345_0134 - pM[SA03]*mDet4_1345_0124 + pM[SA04]*mDet4_1345_0123;
734 const Double_t mDet5_01345_01235 = pM[SA00]*mDet4_1345_1235 - pM[SA01]*mDet4_1345_0235
735 + pM[SA02]*mDet4_1345_0135 - pM[SA03]*mDet4_1345_0125 + pM[SA05]*mDet4_1345_0123;
736 const Double_t mDet5_01345_01245 = pM[SA00]*mDet4_1345_1245 - pM[SA01]*mDet4_1345_0245
737 + pM[SA02]*mDet4_1345_0145 - pM[SA04]*mDet4_1345_0125 + pM[SA05]*mDet4_1345_0124;
738 const Double_t mDet5_01345_01345 = pM[SA00]*mDet4_1345_1345 - pM[SA01]*mDet4_1345_0345
739 + pM[SA03]*mDet4_1345_0145 - pM[SA04]*mDet4_1345_0135 + pM[SA05]*mDet4_1345_0134;
740 const Double_t mDet5_02345_01234 = pM[SA00]*mDet4_2345_1234 - pM[SA01]*mDet4_2345_0234
741 + pM[SA02]*mDet4_2345_0134 - pM[SA03]*mDet4_2345_0124 + pM[SA04]*mDet4_2345_0123;
742 const Double_t mDet5_02345_01235 = pM[SA00]*mDet4_2345_1235 - pM[SA01]*mDet4_2345_0235
743 + pM[SA02]*mDet4_2345_0135 - pM[SA03]*mDet4_2345_0125 + pM[SA05]*mDet4_2345_0123;
744 const Double_t mDet5_02345_01245 = pM[SA00]*mDet4_2345_1245 - pM[SA01]*mDet4_2345_0245
745 + pM[SA02]*mDet4_2345_0145 - pM[SA04]*mDet4_2345_0125 + pM[SA05]*mDet4_2345_0124;
746 const Double_t mDet5_02345_01345 = pM[SA00]*mDet4_2345_1345 - pM[SA01]*mDet4_2345_0345
747 + pM[SA03]*mDet4_2345_0145 - pM[SA04]*mDet4_2345_0135 + pM[SA05]*mDet4_2345_0134;
748 const Double_t mDet5_02345_02345 = pM[SA00]*mDet4_2345_2345 - pM[SA02]*mDet4_2345_0345
749 + pM[SA03]*mDet4_2345_0245 - pM[SA04]*mDet4_2345_0235 + pM[SA05]*mDet4_2345_0234;
750 const Double_t mDet5_12345_01234 = pM[SA10]*mDet4_2345_1234 - pM[SA11]*mDet4_2345_0234
751 + pM[SA12]*mDet4_2345_0134 - pM[SA13]*mDet4_2345_0124 + pM[SA14]*mDet4_2345_0123;
752 const Double_t mDet5_12345_01235 = pM[SA10]*mDet4_2345_1235 - pM[SA11]*mDet4_2345_0235
753 + pM[SA12]*mDet4_2345_0135 - pM[SA13]*mDet4_2345_0125 + pM[SA15]*mDet4_2345_0123;
754 const Double_t mDet5_12345_01245 = pM[SA10]*mDet4_2345_1245 - pM[SA11]*mDet4_2345_0245
755 + pM[SA12]*mDet4_2345_0145 - pM[SA14]*mDet4_2345_0125 + pM[SA15]*mDet4_2345_0124;
756 const Double_t mDet5_12345_01345 = pM[SA10]*mDet4_2345_1345 - pM[SA11]*mDet4_2345_0345
757 + pM[SA13]*mDet4_2345_0145 - pM[SA14]*mDet4_2345_0135 + pM[SA15]*mDet4_2345_0134;
758 const Double_t mDet5_12345_02345 = pM[SA10]*mDet4_2345_2345 - pM[SA12]*mDet4_2345_0345
759 + pM[SA13]*mDet4_2345_0245 - pM[SA14]*mDet4_2345_0235 + pM[SA15]*mDet4_2345_0234;
760 const Double_t mDet5_12345_12345 = pM[SA11]*mDet4_2345_2345 - pM[SA12]*mDet4_2345_1345
761 + pM[SA13]*mDet4_2345_1245 - pM[SA14]*mDet4_2345_1235 + pM[SA15]*mDet4_2345_1234;
762
763 // Find the determinant
764
765 const Double_t det = pM[SA00]*mDet5_12345_12345 - pM[SA01]*mDet5_12345_02345 + pM[SA02]*mDet5_12345_01345
766 - pM[SA03]*mDet5_12345_01245 + pM[SA04]*mDet5_12345_01235 - pM[SA05]*mDet5_12345_01234;
767
768 if (determ)
769 *determ = det;
770
771 if ( det == 0 ) {
772 Error("Inv6x6","matrix is singular");
773 return kFALSE;
774 }
775
776 const Double_t oneOverDet = 1.0/det;
777 const Double_t mn1OverDet = - oneOverDet;
778
779 pM[SA00] = mDet5_12345_12345*oneOverDet;
780 pM[SA01] = mDet5_12345_02345*mn1OverDet;
781 pM[SA02] = mDet5_12345_01345*oneOverDet;
782 pM[SA03] = mDet5_12345_01245*mn1OverDet;
783 pM[SA04] = mDet5_12345_01235*oneOverDet;
784 pM[SA05] = mDet5_12345_01234*mn1OverDet;
785
786 pM[SA11] = mDet5_02345_02345*oneOverDet;
787 pM[SA12] = mDet5_02345_01345*mn1OverDet;
788 pM[SA13] = mDet5_02345_01245*oneOverDet;
789 pM[SA14] = mDet5_02345_01235*mn1OverDet;
790 pM[SA15] = mDet5_02345_01234*oneOverDet;
791
792 pM[SA22] = mDet5_01345_01345*oneOverDet;
793 pM[SA23] = mDet5_01345_01245*mn1OverDet;
794 pM[SA24] = mDet5_01345_01235*oneOverDet;
795 pM[SA25] = mDet5_01345_01234*mn1OverDet;
796
797 pM[SA33] = mDet5_01245_01245*oneOverDet;
798 pM[SA34] = mDet5_01245_01235*mn1OverDet;
799 pM[SA35] = mDet5_01245_01234*oneOverDet;
800
801 pM[SA44] = mDet5_01235_01235*oneOverDet;
802 pM[SA45] = mDet5_01235_01234*mn1OverDet;
803
804 pM[SA55] = mDet5_01234_01234*oneOverDet;
805
806 for (Int_t irow = 0; irow < 6; irow++) {
807 const Int_t rowOff1 = irow*6;
808 for (Int_t icol = 0; icol < irow; icol++) {
809 const Int_t rowOff2 = icol*6;
810 pM[rowOff1+icol] = pM[rowOff2+irow];
811 }
812 }
813
814 return kTRUE;
815}
816
817#include "TMatrixFSymfwd.h"
818
819template Bool_t TMatrixTSymCramerInv::Inv2x2<Float_t>(TMatrixFSym&,Double_t*);
820template Bool_t TMatrixTSymCramerInv::Inv3x3<Float_t>(TMatrixFSym&,Double_t*);
821template Bool_t TMatrixTSymCramerInv::Inv4x4<Float_t>(TMatrixFSym&,Double_t*);
822template Bool_t TMatrixTSymCramerInv::Inv5x5<Float_t>(TMatrixFSym&,Double_t*);
823template Bool_t TMatrixTSymCramerInv::Inv6x6<Float_t>(TMatrixFSym&,Double_t*);
824
825#include "TMatrixDSymfwd.h"
826
827template Bool_t TMatrixTSymCramerInv::Inv2x2<Double_t>(TMatrixDSym&,Double_t*);
828template Bool_t TMatrixTSymCramerInv::Inv3x3<Double_t>(TMatrixDSym&,Double_t*);
829template Bool_t TMatrixTSymCramerInv::Inv4x4<Double_t>(TMatrixDSym&,Double_t*);
830template Bool_t TMatrixTSymCramerInv::Inv5x5<Double_t>(TMatrixDSym&,Double_t*);
831template Bool_t TMatrixTSymCramerInv::Inv6x6<Double_t>(TMatrixDSym&,Double_t*);
Int_t
int Int_t
Definition: RtypesCore.h:41
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:88
Bool_t
bool Bool_t
Definition: RtypesCore.h:59
Double_t
double Double_t
Definition: RtypesCore.h:55
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:87
NamespaceImp
#define NamespaceImp(name)
Definition: Rtypes.h:380
Error
void Error(const char *location, const char *msgfmt,...)
SM33
#define SM33
Definition: TMatrixTSymCramerInv.cxx:278
SA15
#define SA15
Definition: TMatrixTSymCramerInv.cxx:453
SA25
#define SA25
Definition: TMatrixTSymCramerInv.cxx:460
SA31
#define SA31
Definition: TMatrixTSymCramerInv.cxx:463
SA53
#define SA53
Definition: TMatrixTSymCramerInv.cxx:479
SM12
#define SM12
Definition: TMatrixTSymCramerInv.cxx:265
SM13
#define SM13
Definition: TMatrixTSymCramerInv.cxx:266
SM40
#define SM40
Definition: TMatrixTSymCramerInv.cxx:281
SA45
#define SA45
Definition: TMatrixTSymCramerInv.cxx:474
SA05
#define SA05
Definition: TMatrixTSymCramerInv.cxx:446
SF03
#define SF03
Definition: TMatrixTSymCramerInv.cxx:145
SA13
#define SA13
Definition: TMatrixTSymCramerInv.cxx:451
SA43
#define SA43
Definition: TMatrixTSymCramerInv.cxx:472
SM14
#define SM14
Definition: TMatrixTSymCramerInv.cxx:267
SA54
#define SA54
Definition: TMatrixTSymCramerInv.cxx:480
SA01
#define SA01
Definition: TMatrixTSymCramerInv.cxx:442
SF02
#define SF02
Definition: TMatrixTSymCramerInv.cxx:144
SA14
#define SA14
Definition: TMatrixTSymCramerInv.cxx:452
SF31
#define SF31
Definition: TMatrixTSymCramerInv.cxx:158
SA21
#define SA21
Definition: TMatrixTSymCramerInv.cxx:456
SF00
#define SF00
Definition: TMatrixTSymCramerInv.cxx:142
SF13
#define SF13
Definition: TMatrixTSymCramerInv.cxx:150
SA24
#define SA24
Definition: TMatrixTSymCramerInv.cxx:459
SA02
#define SA02
Definition: TMatrixTSymCramerInv.cxx:443
SF22
#define SF22
Definition: TMatrixTSymCramerInv.cxx:154
SA52
#define SA52
Definition: TMatrixTSymCramerInv.cxx:478
SF32
#define SF32
Definition: TMatrixTSymCramerInv.cxx:159
SM23
#define SM23
Definition: TMatrixTSymCramerInv.cxx:272
SF10
#define SF10
Definition: TMatrixTSymCramerInv.cxx:147
SM22
#define SM22
Definition: TMatrixTSymCramerInv.cxx:271
SA22
#define SA22
Definition: TMatrixTSymCramerInv.cxx:457
SF33
#define SF33
Definition: TMatrixTSymCramerInv.cxx:160
SM10
#define SM10
Definition: TMatrixTSymCramerInv.cxx:263
SM00
#define SM00
Definition: TMatrixTSymCramerInv.cxx:257
SF21
#define SF21
Definition: TMatrixTSymCramerInv.cxx:153
SA33
#define SA33
Definition: TMatrixTSymCramerInv.cxx:465
SA40
#define SA40
Definition: TMatrixTSymCramerInv.cxx:469
SM30
#define SM30
Definition: TMatrixTSymCramerInv.cxx:275
SF12
#define SF12
Definition: TMatrixTSymCramerInv.cxx:149
SA41
#define SA41
Definition: TMatrixTSymCramerInv.cxx:470
SM01
#define SM01
Definition: TMatrixTSymCramerInv.cxx:258
SF23
#define SF23
Definition: TMatrixTSymCramerInv.cxx:155
SM02
#define SM02
Definition: TMatrixTSymCramerInv.cxx:259
SF20
#define SF20
Definition: TMatrixTSymCramerInv.cxx:152
SM41
#define SM41
Definition: TMatrixTSymCramerInv.cxx:282
SA04
#define SA04
Definition: TMatrixTSymCramerInv.cxx:445
SA20
#define SA20
Definition: TMatrixTSymCramerInv.cxx:455
SM21
#define SM21
Definition: TMatrixTSymCramerInv.cxx:270
SA00
#define SA00
Definition: TMatrixTSymCramerInv.cxx:441
SA23
#define SA23
Definition: TMatrixTSymCramerInv.cxx:458
SM32
#define SM32
Definition: TMatrixTSymCramerInv.cxx:277
SM11
#define SM11
Definition: TMatrixTSymCramerInv.cxx:264
SA30
#define SA30
Definition: TMatrixTSymCramerInv.cxx:462
SM20
#define SM20
Definition: TMatrixTSymCramerInv.cxx:269
SM44
#define SM44
Definition: TMatrixTSymCramerInv.cxx:285
SA11
#define SA11
Definition: TMatrixTSymCramerInv.cxx:449
SA03
#define SA03
Definition: TMatrixTSymCramerInv.cxx:444
SA10
#define SA10
Definition: TMatrixTSymCramerInv.cxx:448
SM04
#define SM04
Definition: TMatrixTSymCramerInv.cxx:261
SF30
#define SF30
Definition: TMatrixTSymCramerInv.cxx:157
SF11
#define SF11
Definition: TMatrixTSymCramerInv.cxx:148
SA35
#define SA35
Definition: TMatrixTSymCramerInv.cxx:467
SM43
#define SM43
Definition: TMatrixTSymCramerInv.cxx:284
SA12
#define SA12
Definition: TMatrixTSymCramerInv.cxx:450
SA34
#define SA34
Definition: TMatrixTSymCramerInv.cxx:466
SM31
#define SM31
Definition: TMatrixTSymCramerInv.cxx:276
SM03
#define SM03
Definition: TMatrixTSymCramerInv.cxx:260
SM34
#define SM34
Definition: TMatrixTSymCramerInv.cxx:279
SM42
#define SM42
Definition: TMatrixTSymCramerInv.cxx:283
SA51
#define SA51
Definition: TMatrixTSymCramerInv.cxx:477
SA32
#define SA32
Definition: TMatrixTSymCramerInv.cxx:464
SA42
#define SA42
Definition: TMatrixTSymCramerInv.cxx:471
SA55
#define SA55
Definition: TMatrixTSymCramerInv.cxx:481
SA44
#define SA44
Definition: TMatrixTSymCramerInv.cxx:473
SA50
#define SA50
Definition: TMatrixTSymCramerInv.cxx:476
SF01
#define SF01
Definition: TMatrixTSymCramerInv.cxx:143
SM24
#define SM24
Definition: TMatrixTSymCramerInv.cxx:273
TMatrixTSym
TMatrixTSym.
Definition: TMatrixTSym.h:34
TGeant4Unit::s
static constexpr double s
Definition: TGeant4SystemOfUnits.h:162
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
TMatrixTSymCramerInv::Inv6x6
Bool_t Inv6x6(TMatrixTSym< Element > &m, Double_t *determ)
TMatrixTSymCramerInv::Inv3x3
Bool_t Inv3x3(TMatrixTSym< Element > &m, Double_t *determ)
TMatrixTSymCramerInv::Inv4x4
Bool_t Inv4x4(TMatrixTSym< Element > &m, Double_t *determ)
TMatrixTSymCramerInv::Inv2x2
Bool_t Inv2x2(TMatrixTSym< Element > &m, Double_t *determ)
TMatrixTSymCramerInv::Inv5x5
Bool_t Inv5x5(TMatrixTSym< Element > &m, Double_t *determ)
m
auto * m
Definition: textangle.C:8
t1
auto * t1
Definition: textangle.C:20