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