Logo ROOT   6.13/01
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)
43 NamespaceImp(TMatrixTSymCramerInv);
44 #endif
45 
46 ////////////////////////////////////////////////////////////////////////////////
47 
48 template<class Element>
49 Bool_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 
78 template<class Element>
79 Bool_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 
164 template<class Element>
165 Bool_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 
289 template<class Element>
290 Bool_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 
485 template<class Element>
486 Bool_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 
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 #include "TMatrixDSymfwd.h"
826 
827 template Bool_t TMatrixTSymCramerInv::Inv2x2<Double_t>(TMatrixDSym&,Double_t*);
828 template Bool_t TMatrixTSymCramerInv::Inv3x3<Double_t>(TMatrixDSym&,Double_t*);
829 template Bool_t TMatrixTSymCramerInv::Inv4x4<Double_t>(TMatrixDSym&,Double_t*);
830 template Bool_t TMatrixTSymCramerInv::Inv5x5<Double_t>(TMatrixDSym&,Double_t*);
831 template Bool_t TMatrixTSymCramerInv::Inv6x6<Double_t>(TMatrixDSym&,Double_t*);
SM01
#define SM01
Definition: TMatrixTSymCramerInv.cxx:258
SM43
#define SM43
Definition: TMatrixTSymCramerInv.cxx:284
TMatrixTBase::GetRowLwb
Int_t GetRowLwb() const
Definition: TMatrixTBase.h:120
SA21
#define SA21
Definition: TMatrixTSymCramerInv.cxx:456
SA45
#define SA45
Definition: TMatrixTSymCramerInv.cxx:474
SM13
#define SM13
Definition: TMatrixTSymCramerInv.cxx:266
SF00
#define SF00
Definition: TMatrixTSymCramerInv.cxx:142
SM14
#define SM14
Definition: TMatrixTSymCramerInv.cxx:267
SM44
#define SM44
Definition: TMatrixTSymCramerInv.cxx:285
SM04
#define SM04
Definition: TMatrixTSymCramerInv.cxx:261
TMatrixTBase::GetNcols
Int_t GetNcols() const
Definition: TMatrixTBase.h:125
SM23
#define SM23
Definition: TMatrixTSymCramerInv.cxx:272
SM02
#define SM02
Definition: TMatrixTSymCramerInv.cxx:259
SA11
#define SA11
Definition: TMatrixTSymCramerInv.cxx:449
SM30
#define SM30
Definition: TMatrixTSymCramerInv.cxx:275
SA20
#define SA20
Definition: TMatrixTSymCramerInv.cxx:455
SA32
#define SA32
Definition: TMatrixTSymCramerInv.cxx:464
SF01
#define SF01
Definition: TMatrixTSymCramerInv.cxx:143
SA00
#define SA00
Definition: TMatrixTSymCramerInv.cxx:441
SF30
#define SF30
Definition: TMatrixTSymCramerInv.cxx:157
SM20
#define SM20
Definition: TMatrixTSymCramerInv.cxx:269
SA53
#define SA53
Definition: TMatrixTSymCramerInv.cxx:479
SA41
#define SA41
Definition: TMatrixTSymCramerInv.cxx:470
TMatrixTSymCramerInv::Inv6x6
Bool_t Inv6x6(TMatrixTSym< Element > &m, Double_t *determ)
TMatrixTSym::GetMatrixArray
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSym.h:189
SA01
#define SA01
Definition: TMatrixTSymCramerInv.cxx:442
SA40
#define SA40
Definition: TMatrixTSymCramerInv.cxx:469
SA44
#define SA44
Definition: TMatrixTSymCramerInv.cxx:473
SA02
#define SA02
Definition: TMatrixTSymCramerInv.cxx:443
SF23
#define SF23
Definition: TMatrixTSymCramerInv.cxx:155
SM31
#define SM31
Definition: TMatrixTSymCramerInv.cxx:276
SM32
#define SM32
Definition: TMatrixTSymCramerInv.cxx:277
SA42
#define SA42
Definition: TMatrixTSymCramerInv.cxx:471
SM12
#define SM12
Definition: TMatrixTSymCramerInv.cxx:265
SF32
#define SF32
Definition: TMatrixTSymCramerInv.cxx:159
TMatrixTBase::GetColLwb
Int_t GetColLwb() const
Definition: TMatrixTBase.h:123
SA23
#define SA23
Definition: TMatrixTSymCramerInv.cxx:458
SA35
#define SA35
Definition: TMatrixTSymCramerInv.cxx:467
SM03
#define SM03
Definition: TMatrixTSymCramerInv.cxx:260
SA03
#define SA03
Definition: TMatrixTSymCramerInv.cxx:444
SM22
#define SM22
Definition: TMatrixTSymCramerInv.cxx:271
SA25
#define SA25
Definition: TMatrixTSymCramerInv.cxx:460
SM24
#define SM24
Definition: TMatrixTSymCramerInv.cxx:273
TMatrixTSymCramerInv::Inv3x3
Bool_t Inv3x3(TMatrixTSym< Element > &m, Double_t *determ)
TMatrixTSym
TMatrixTSym.
Definition: TMatrixDSymfwd.h:22
SF11
#define SF11
Definition: TMatrixTSymCramerInv.cxx:148
SA22
#define SA22
Definition: TMatrixTSymCramerInv.cxx:457
m
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t Int_t m
Definition: TRolke.cxx:637
SA50
#define SA50
Definition: TMatrixTSymCramerInv.cxx:476
SM34
#define SM34
Definition: TMatrixTSymCramerInv.cxx:279
SF10
#define SF10
Definition: TMatrixTSymCramerInv.cxx:147
TMatrixTSymCramerInv::Inv4x4
Bool_t Inv4x4(TMatrixTSym< Element > &m, Double_t *determ)
SM33
#define SM33
Definition: TMatrixTSymCramerInv.cxx:278
SF33
#define SF33
Definition: TMatrixTSymCramerInv.cxx:160
TMatrixTSymCramerInv::Inv5x5
Bool_t Inv5x5(TMatrixTSym< Element > &m, Double_t *determ)
SA52
#define SA52
Definition: TMatrixTSymCramerInv.cxx:478
SF02
#define SF02
Definition: TMatrixTSymCramerInv.cxx:144
SA54
#define SA54
Definition: TMatrixTSymCramerInv.cxx:480
SF31
#define SF31
Definition: TMatrixTSymCramerInv.cxx:158
TMatrixTBase::GetNrows
Int_t GetNrows() const
Definition: TMatrixTBase.h:122
SF13
#define SF13
Definition: TMatrixTSymCramerInv.cxx:150
SA05
#define SA05
Definition: TMatrixTSymCramerInv.cxx:446
SA13
#define SA13
Definition: TMatrixTSymCramerInv.cxx:451
SM42
#define SM42
Definition: TMatrixTSymCramerInv.cxx:283
SF03
#define SF03
Definition: TMatrixTSymCramerInv.cxx:145
SA15
#define SA15
Definition: TMatrixTSymCramerInv.cxx:453
SF12
#define SF12
Definition: TMatrixTSymCramerInv.cxx:149
SA31
#define SA31
Definition: TMatrixTSymCramerInv.cxx:463
SA30
#define SA30
Definition: TMatrixTSymCramerInv.cxx:462
SF21
#define SF21
Definition: TMatrixTSymCramerInv.cxx:153
SF22
#define SF22
Definition: TMatrixTSymCramerInv.cxx:154
SA33
#define SA33
Definition: TMatrixTSymCramerInv.cxx:465
SM11
#define SM11
Definition: TMatrixTSymCramerInv.cxx:264
SA12
#define SA12
Definition: TMatrixTSymCramerInv.cxx:450
SA10
#define SA10
Definition: TMatrixTSymCramerInv.cxx:448
SM00
#define SM00
Definition: TMatrixTSymCramerInv.cxx:257
SF20
#define SF20
Definition: TMatrixTSymCramerInv.cxx:152
SM10
#define SM10
Definition: TMatrixTSymCramerInv.cxx:263
SM21
#define SM21
Definition: TMatrixTSymCramerInv.cxx:270
SA51
#define SA51
Definition: TMatrixTSymCramerInv.cxx:477
SA04
#define SA04
Definition: TMatrixTSymCramerInv.cxx:445
SA55
#define SA55
Definition: TMatrixTSymCramerInv.cxx:481
TMatrixTSymCramerInv::Inv2x2
Bool_t Inv2x2(TMatrixTSym< Element > &m, Double_t *determ)
SA14
#define SA14
Definition: TMatrixTSymCramerInv.cxx:452
SM41
#define SM41
Definition: TMatrixTSymCramerInv.cxx:282
SA24
#define SA24
Definition: TMatrixTSymCramerInv.cxx:459
SM40
#define SM40
Definition: TMatrixTSymCramerInv.cxx:281
SA43
#define SA43
Definition: TMatrixTSymCramerInv.cxx:472
SA34
#define SA34
Definition: TMatrixTSymCramerInv.cxx:466