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