CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

MatrixInvert.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
4 //
5 // This software written by Mark Fischler and Steven Haywood
6 //
7 
8 #ifdef GNUPRAGMA
9 #pragma implementation
10 #endif
11 
12 #include <string.h>
13 
14 #include "CLHEP/Matrix/defs.h"
15 #include "CLHEP/Matrix/Matrix.h"
16 
17 namespace CLHEP {
18 
19 // Aij are indices for a 6x6 matrix.
20 // Mij are indices for a 5x5 matrix.
21 // Fij are indices for a 4x4 matrix.
22 
23 #define A00 0
24 #define A01 1
25 #define A02 2
26 #define A03 3
27 #define A04 4
28 #define A05 5
29 
30 #define A10 6
31 #define A11 7
32 #define A12 8
33 #define A13 9
34 #define A14 10
35 #define A15 11
36 
37 #define A20 12
38 #define A21 13
39 #define A22 14
40 #define A23 15
41 #define A24 16
42 #define A25 17
43 
44 #define A30 18
45 #define A31 19
46 #define A32 20
47 #define A33 21
48 #define A34 22
49 #define A35 23
50 
51 #define A40 24
52 #define A41 25
53 #define A42 26
54 #define A43 27
55 #define A44 28
56 #define A45 29
57 
58 #define A50 30
59 #define A51 31
60 #define A52 32
61 #define A53 33
62 #define A54 34
63 #define A55 35
64 
65 #define M00 0
66 #define M01 1
67 #define M02 2
68 #define M03 3
69 #define M04 4
70 
71 #define M10 5
72 #define M11 6
73 #define M12 7
74 #define M13 8
75 #define M14 9
76 
77 #define M20 10
78 #define M21 11
79 #define M22 12
80 #define M23 13
81 #define M24 14
82 
83 #define M30 15
84 #define M31 16
85 #define M32 17
86 #define M33 18
87 #define M34 19
88 
89 #define M40 20
90 #define M41 21
91 #define M42 22
92 #define M43 23
93 #define M44 24
94 
95 #define F00 0
96 #define F01 1
97 #define F02 2
98 #define F03 3
99 
100 #define F10 4
101 #define F11 5
102 #define F12 6
103 #define F13 7
104 
105 #define F20 8
106 #define F21 9
107 #define F22 10
108 #define F23 11
109 
110 #define F30 12
111 #define F31 13
112 #define F32 14
113 #define F33 15
114 
115 
116 void HepMatrix::invertHaywood4 (int & ifail) {
117 
118  ifail = 0;
119 
120  // Find all NECESSARY 2x2 dets: (18 of them)
121 
122  double Det2_12_01 = m[F10]*m[F21] - m[F11]*m[F20];
123  double Det2_12_02 = m[F10]*m[F22] - m[F12]*m[F20];
124  double Det2_12_03 = m[F10]*m[F23] - m[F13]*m[F20]; //
125  double Det2_12_13 = m[F11]*m[F23] - m[F13]*m[F21]; //
126  double Det2_12_23 = m[F12]*m[F23] - m[F13]*m[F22]; //
127  double Det2_12_12 = m[F11]*m[F22] - m[F12]*m[F21];
128  double Det2_13_01 = m[F10]*m[F31] - m[F11]*m[F30];
129  double Det2_13_02 = m[F10]*m[F32] - m[F12]*m[F30];
130  double Det2_13_03 = m[F10]*m[F33] - m[F13]*m[F30];
131  double Det2_13_12 = m[F11]*m[F32] - m[F12]*m[F31];
132  double Det2_13_13 = m[F11]*m[F33] - m[F13]*m[F31];
133  double Det2_13_23 = m[F12]*m[F33] - m[F13]*m[F32]; //
134  double Det2_23_01 = m[F20]*m[F31] - m[F21]*m[F30];
135  double Det2_23_02 = m[F20]*m[F32] - m[F22]*m[F30];
136  double Det2_23_03 = m[F20]*m[F33] - m[F23]*m[F30];
137  double Det2_23_12 = m[F21]*m[F32] - m[F22]*m[F31];
138  double Det2_23_13 = m[F21]*m[F33] - m[F23]*m[F31];
139  double Det2_23_23 = m[F22]*m[F33] - m[F23]*m[F32];
140 
141  // Find all NECESSARY 3x3 dets: (16 of them)
142 
143  double Det3_012_012 = m[F00]*Det2_12_12 - m[F01]*Det2_12_02
144  + m[F02]*Det2_12_01;
145  double Det3_012_013 = m[F00]*Det2_12_13 - m[F01]*Det2_12_03
146  + m[F03]*Det2_12_01; //
147  double Det3_012_023 = m[F00]*Det2_12_23 - m[F02]*Det2_12_03
148  + m[F03]*Det2_12_02; //
149  double Det3_012_123 = m[F01]*Det2_12_23 - m[F02]*Det2_12_13
150  + m[F03]*Det2_12_12; //
151  double Det3_013_012 = m[F00]*Det2_13_12 - m[F01]*Det2_13_02
152  + m[F02]*Det2_13_01;
153  double Det3_013_013 = m[F00]*Det2_13_13 - m[F01]*Det2_13_03
154  + m[F03]*Det2_13_01;
155  double Det3_013_023 = m[F00]*Det2_13_23 - m[F02]*Det2_13_03
156  + m[F03]*Det2_13_02; //
157  double Det3_013_123 = m[F01]*Det2_13_23 - m[F02]*Det2_13_13
158  + m[F03]*Det2_13_12; //
159  double Det3_023_012 = m[F00]*Det2_23_12 - m[F01]*Det2_23_02
160  + m[F02]*Det2_23_01;
161  double Det3_023_013 = m[F00]*Det2_23_13 - m[F01]*Det2_23_03
162  + m[F03]*Det2_23_01;
163  double Det3_023_023 = m[F00]*Det2_23_23 - m[F02]*Det2_23_03
164  + m[F03]*Det2_23_02;
165  double Det3_023_123 = m[F01]*Det2_23_23 - m[F02]*Det2_23_13
166  + m[F03]*Det2_23_12; //
167  double Det3_123_012 = m[F10]*Det2_23_12 - m[F11]*Det2_23_02
168  + m[F12]*Det2_23_01;
169  double Det3_123_013 = m[F10]*Det2_23_13 - m[F11]*Det2_23_03
170  + m[F13]*Det2_23_01;
171  double Det3_123_023 = m[F10]*Det2_23_23 - m[F12]*Det2_23_03
172  + m[F13]*Det2_23_02;
173  double Det3_123_123 = m[F11]*Det2_23_23 - m[F12]*Det2_23_13
174  + m[F13]*Det2_23_12;
175 
176  // Find the 4x4 det:
177 
178  double det = m[F00]*Det3_123_123
179  - m[F01]*Det3_123_023
180  + m[F02]*Det3_123_013
181  - m[F03]*Det3_123_012;
182 
183  if ( det == 0 ) {
184 #ifdef SINGULAR_DIAGNOSTICS
185  std::cerr << "Kramer's rule inversion of a singular 4x4 matrix: "
186  << *this << "\n";
187 #endif
188  ifail = 1;
189  return;
190  }
191 
192  double oneOverDet = 1.0/det;
193  double mn1OverDet = - oneOverDet;
194 
195  m[F00] = Det3_123_123 * oneOverDet;
196  m[F01] = Det3_023_123 * mn1OverDet;
197  m[F02] = Det3_013_123 * oneOverDet;
198  m[F03] = Det3_012_123 * mn1OverDet;
199 
200  m[F10] = Det3_123_023 * mn1OverDet;
201  m[F11] = Det3_023_023 * oneOverDet;
202  m[F12] = Det3_013_023 * mn1OverDet;
203  m[F13] = Det3_012_023 * oneOverDet;
204 
205  m[F20] = Det3_123_013 * oneOverDet;
206  m[F21] = Det3_023_013 * mn1OverDet;
207  m[F22] = Det3_013_013 * oneOverDet;
208  m[F23] = Det3_012_013 * mn1OverDet;
209 
210  m[F30] = Det3_123_012 * mn1OverDet;
211  m[F31] = Det3_023_012 * oneOverDet;
212  m[F32] = Det3_013_012 * mn1OverDet;
213  m[F33] = Det3_012_012 * oneOverDet;
214 
215  return;
216 }
217 
218 
219 
220 void HepMatrix::invertHaywood5 (int & ifail) {
221 
222  ifail = 0;
223 
224  // Find all NECESSARY 2x2 dets: (30 of them)
225 
226  double Det2_23_01 = m[M20]*m[M31] - m[M21]*m[M30];
227  double Det2_23_02 = m[M20]*m[M32] - m[M22]*m[M30];
228  double Det2_23_03 = m[M20]*m[M33] - m[M23]*m[M30];
229  double Det2_23_04 = m[M20]*m[M34] - m[M24]*m[M30];
230  double Det2_23_12 = m[M21]*m[M32] - m[M22]*m[M31]; //
231  double Det2_23_13 = m[M21]*m[M33] - m[M23]*m[M31];
232  double Det2_23_14 = m[M21]*m[M34] - m[M24]*m[M31]; //
233  double Det2_23_23 = m[M22]*m[M33] - m[M23]*m[M32];
234  double Det2_23_24 = m[M22]*m[M34] - m[M24]*m[M32]; //
235  double Det2_23_34 = m[M23]*m[M34] - m[M24]*m[M33]; //
236  double Det2_24_01 = m[M20]*m[M41] - m[M21]*m[M40];
237  double Det2_24_02 = m[M20]*m[M42] - m[M22]*m[M40];
238  double Det2_24_03 = m[M20]*m[M43] - m[M23]*m[M40];
239  double Det2_24_04 = m[M20]*m[M44] - m[M24]*m[M40];
240  double Det2_24_12 = m[M21]*m[M42] - m[M22]*m[M41];
241  double Det2_24_13 = m[M21]*m[M43] - m[M23]*m[M41];
242  double Det2_24_14 = m[M21]*m[M44] - m[M24]*m[M41];
243  double Det2_24_23 = m[M22]*m[M43] - m[M23]*m[M42];
244  double Det2_24_24 = m[M22]*m[M44] - m[M24]*m[M42];
245  double Det2_24_34 = m[M23]*m[M44] - m[M24]*m[M43]; //
246  double Det2_34_01 = m[M30]*m[M41] - m[M31]*m[M40];
247  double Det2_34_02 = m[M30]*m[M42] - m[M32]*m[M40];
248  double Det2_34_03 = m[M30]*m[M43] - m[M33]*m[M40];
249  double Det2_34_04 = m[M30]*m[M44] - m[M34]*m[M40];
250  double Det2_34_12 = m[M31]*m[M42] - m[M32]*m[M41];
251  double Det2_34_13 = m[M31]*m[M43] - m[M33]*m[M41];
252  double Det2_34_14 = m[M31]*m[M44] - m[M34]*m[M41];
253  double Det2_34_23 = m[M32]*m[M43] - m[M33]*m[M42];
254  double Det2_34_24 = m[M32]*m[M44] - m[M34]*m[M42];
255  double Det2_34_34 = m[M33]*m[M44] - m[M34]*m[M43];
256 
257  // Find all NECESSARY 3x3 dets: (40 of them)
258 
259  double Det3_123_012 = m[M10]*Det2_23_12 - m[M11]*Det2_23_02
260  + m[M12]*Det2_23_01;
261  double Det3_123_013 = m[M10]*Det2_23_13 - m[M11]*Det2_23_03
262  + m[M13]*Det2_23_01;
263  double Det3_123_014 = m[M10]*Det2_23_14 - m[M11]*Det2_23_04
264  + m[M14]*Det2_23_01; //
265  double Det3_123_023 = m[M10]*Det2_23_23 - m[M12]*Det2_23_03
266  + m[M13]*Det2_23_02;
267  double Det3_123_024 = m[M10]*Det2_23_24 - m[M12]*Det2_23_04
268  + m[M14]*Det2_23_02; //
269  double Det3_123_034 = m[M10]*Det2_23_34 - m[M13]*Det2_23_04
270  + m[M14]*Det2_23_03; //
271  double Det3_123_123 = m[M11]*Det2_23_23 - m[M12]*Det2_23_13
272  + m[M13]*Det2_23_12;
273  double Det3_123_124 = m[M11]*Det2_23_24 - m[M12]*Det2_23_14
274  + m[M14]*Det2_23_12; //
275  double Det3_123_134 = m[M11]*Det2_23_34 - m[M13]*Det2_23_14
276  + m[M14]*Det2_23_13; //
277  double Det3_123_234 = m[M12]*Det2_23_34 - m[M13]*Det2_23_24
278  + m[M14]*Det2_23_23; //
279  double Det3_124_012 = m[M10]*Det2_24_12 - m[M11]*Det2_24_02
280  + m[M12]*Det2_24_01;
281  double Det3_124_013 = m[M10]*Det2_24_13 - m[M11]*Det2_24_03
282  + m[M13]*Det2_24_01;
283  double Det3_124_014 = m[M10]*Det2_24_14 - m[M11]*Det2_24_04
284  + m[M14]*Det2_24_01;
285  double Det3_124_023 = m[M10]*Det2_24_23 - m[M12]*Det2_24_03
286  + m[M13]*Det2_24_02;
287  double Det3_124_024 = m[M10]*Det2_24_24 - m[M12]*Det2_24_04
288  + m[M14]*Det2_24_02;
289  double Det3_124_034 = m[M10]*Det2_24_34 - m[M13]*Det2_24_04
290  + m[M14]*Det2_24_03; //
291  double Det3_124_123 = m[M11]*Det2_24_23 - m[M12]*Det2_24_13
292  + m[M13]*Det2_24_12;
293  double Det3_124_124 = m[M11]*Det2_24_24 - m[M12]*Det2_24_14
294  + m[M14]*Det2_24_12;
295  double Det3_124_134 = m[M11]*Det2_24_34 - m[M13]*Det2_24_14
296  + m[M14]*Det2_24_13; //
297  double Det3_124_234 = m[M12]*Det2_24_34 - m[M13]*Det2_24_24
298  + m[M14]*Det2_24_23; //
299  double Det3_134_012 = m[M10]*Det2_34_12 - m[M11]*Det2_34_02
300  + m[M12]*Det2_34_01;
301  double Det3_134_013 = m[M10]*Det2_34_13 - m[M11]*Det2_34_03
302  + m[M13]*Det2_34_01;
303  double Det3_134_014 = m[M10]*Det2_34_14 - m[M11]*Det2_34_04
304  + m[M14]*Det2_34_01;
305  double Det3_134_023 = m[M10]*Det2_34_23 - m[M12]*Det2_34_03
306  + m[M13]*Det2_34_02;
307  double Det3_134_024 = m[M10]*Det2_34_24 - m[M12]*Det2_34_04
308  + m[M14]*Det2_34_02;
309  double Det3_134_034 = m[M10]*Det2_34_34 - m[M13]*Det2_34_04
310  + m[M14]*Det2_34_03;
311  double Det3_134_123 = m[M11]*Det2_34_23 - m[M12]*Det2_34_13
312  + m[M13]*Det2_34_12;
313  double Det3_134_124 = m[M11]*Det2_34_24 - m[M12]*Det2_34_14
314  + m[M14]*Det2_34_12;
315  double Det3_134_134 = m[M11]*Det2_34_34 - m[M13]*Det2_34_14
316  + m[M14]*Det2_34_13;
317  double Det3_134_234 = m[M12]*Det2_34_34 - m[M13]*Det2_34_24
318  + m[M14]*Det2_34_23; //
319  double Det3_234_012 = m[M20]*Det2_34_12 - m[M21]*Det2_34_02
320  + m[M22]*Det2_34_01;
321  double Det3_234_013 = m[M20]*Det2_34_13 - m[M21]*Det2_34_03
322  + m[M23]*Det2_34_01;
323  double Det3_234_014 = m[M20]*Det2_34_14 - m[M21]*Det2_34_04
324  + m[M24]*Det2_34_01;
325  double Det3_234_023 = m[M20]*Det2_34_23 - m[M22]*Det2_34_03
326  + m[M23]*Det2_34_02;
327  double Det3_234_024 = m[M20]*Det2_34_24 - m[M22]*Det2_34_04
328  + m[M24]*Det2_34_02;
329  double Det3_234_034 = m[M20]*Det2_34_34 - m[M23]*Det2_34_04
330  + m[M24]*Det2_34_03;
331  double Det3_234_123 = m[M21]*Det2_34_23 - m[M22]*Det2_34_13
332  + m[M23]*Det2_34_12;
333  double Det3_234_124 = m[M21]*Det2_34_24 - m[M22]*Det2_34_14
334  + m[M24]*Det2_34_12;
335  double Det3_234_134 = m[M21]*Det2_34_34 - m[M23]*Det2_34_14
336  + m[M24]*Det2_34_13;
337  double Det3_234_234 = m[M22]*Det2_34_34 - m[M23]*Det2_34_24
338  + m[M24]*Det2_34_23;
339 
340  // Find all NECESSARY 4x4 dets: (25 of them)
341 
342  double Det4_0123_0123 = m[M00]*Det3_123_123 - m[M01]*Det3_123_023
343  + m[M02]*Det3_123_013 - m[M03]*Det3_123_012;
344  double Det4_0123_0124 = m[M00]*Det3_123_124 - m[M01]*Det3_123_024
345  + m[M02]*Det3_123_014 - m[M04]*Det3_123_012; //
346  double Det4_0123_0134 = m[M00]*Det3_123_134 - m[M01]*Det3_123_034
347  + m[M03]*Det3_123_014 - m[M04]*Det3_123_013; //
348  double Det4_0123_0234 = m[M00]*Det3_123_234 - m[M02]*Det3_123_034
349  + m[M03]*Det3_123_024 - m[M04]*Det3_123_023; //
350  double Det4_0123_1234 = m[M01]*Det3_123_234 - m[M02]*Det3_123_134
351  + m[M03]*Det3_123_124 - m[M04]*Det3_123_123; //
352  double Det4_0124_0123 = m[M00]*Det3_124_123 - m[M01]*Det3_124_023
353  + m[M02]*Det3_124_013 - m[M03]*Det3_124_012;
354  double Det4_0124_0124 = m[M00]*Det3_124_124 - m[M01]*Det3_124_024
355  + m[M02]*Det3_124_014 - m[M04]*Det3_124_012;
356  double Det4_0124_0134 = m[M00]*Det3_124_134 - m[M01]*Det3_124_034
357  + m[M03]*Det3_124_014 - m[M04]*Det3_124_013; //
358  double Det4_0124_0234 = m[M00]*Det3_124_234 - m[M02]*Det3_124_034
359  + m[M03]*Det3_124_024 - m[M04]*Det3_124_023; //
360  double Det4_0124_1234 = m[M01]*Det3_124_234 - m[M02]*Det3_124_134
361  + m[M03]*Det3_124_124 - m[M04]*Det3_124_123; //
362  double Det4_0134_0123 = m[M00]*Det3_134_123 - m[M01]*Det3_134_023
363  + m[M02]*Det3_134_013 - m[M03]*Det3_134_012;
364  double Det4_0134_0124 = m[M00]*Det3_134_124 - m[M01]*Det3_134_024
365  + m[M02]*Det3_134_014 - m[M04]*Det3_134_012;
366  double Det4_0134_0134 = m[M00]*Det3_134_134 - m[M01]*Det3_134_034
367  + m[M03]*Det3_134_014 - m[M04]*Det3_134_013;
368  double Det4_0134_0234 = m[M00]*Det3_134_234 - m[M02]*Det3_134_034
369  + m[M03]*Det3_134_024 - m[M04]*Det3_134_023; //
370  double Det4_0134_1234 = m[M01]*Det3_134_234 - m[M02]*Det3_134_134
371  + m[M03]*Det3_134_124 - m[M04]*Det3_134_123; //
372  double Det4_0234_0123 = m[M00]*Det3_234_123 - m[M01]*Det3_234_023
373  + m[M02]*Det3_234_013 - m[M03]*Det3_234_012;
374  double Det4_0234_0124 = m[M00]*Det3_234_124 - m[M01]*Det3_234_024
375  + m[M02]*Det3_234_014 - m[M04]*Det3_234_012;
376  double Det4_0234_0134 = m[M00]*Det3_234_134 - m[M01]*Det3_234_034
377  + m[M03]*Det3_234_014 - m[M04]*Det3_234_013;
378  double Det4_0234_0234 = m[M00]*Det3_234_234 - m[M02]*Det3_234_034
379  + m[M03]*Det3_234_024 - m[M04]*Det3_234_023;
380  double Det4_0234_1234 = m[M01]*Det3_234_234 - m[M02]*Det3_234_134
381  + m[M03]*Det3_234_124 - m[M04]*Det3_234_123; //
382  double Det4_1234_0123 = m[M10]*Det3_234_123 - m[M11]*Det3_234_023
383  + m[M12]*Det3_234_013 - m[M13]*Det3_234_012;
384  double Det4_1234_0124 = m[M10]*Det3_234_124 - m[M11]*Det3_234_024
385  + m[M12]*Det3_234_014 - m[M14]*Det3_234_012;
386  double Det4_1234_0134 = m[M10]*Det3_234_134 - m[M11]*Det3_234_034
387  + m[M13]*Det3_234_014 - m[M14]*Det3_234_013;
388  double Det4_1234_0234 = m[M10]*Det3_234_234 - m[M12]*Det3_234_034
389  + m[M13]*Det3_234_024 - m[M14]*Det3_234_023;
390  double Det4_1234_1234 = m[M11]*Det3_234_234 - m[M12]*Det3_234_134
391  + m[M13]*Det3_234_124 - m[M14]*Det3_234_123;
392 
393  // Find the 5x5 det:
394 
395  double det = m[M00]*Det4_1234_1234
396  - m[M01]*Det4_1234_0234
397  + m[M02]*Det4_1234_0134
398  - m[M03]*Det4_1234_0124
399  + m[M04]*Det4_1234_0123;
400 
401  if ( det == 0 ) {
402 #ifdef SINGULAR_DIAGNOSTICS
403  std::cerr << "Kramer's rule inversion of a singular 5x5 matrix: "
404  << *this << "\n";
405 #endif
406  ifail = 1;
407  return;
408  }
409 
410  double oneOverDet = 1.0/det;
411  double mn1OverDet = - oneOverDet;
412 
413  m[M00] = Det4_1234_1234 * oneOverDet;
414  m[M01] = Det4_0234_1234 * mn1OverDet;
415  m[M02] = Det4_0134_1234 * oneOverDet;
416  m[M03] = Det4_0124_1234 * mn1OverDet;
417  m[M04] = Det4_0123_1234 * oneOverDet;
418 
419  m[M10] = Det4_1234_0234 * mn1OverDet;
420  m[M11] = Det4_0234_0234 * oneOverDet;
421  m[M12] = Det4_0134_0234 * mn1OverDet;
422  m[M13] = Det4_0124_0234 * oneOverDet;
423  m[M14] = Det4_0123_0234 * mn1OverDet;
424 
425  m[M20] = Det4_1234_0134 * oneOverDet;
426  m[M21] = Det4_0234_0134 * mn1OverDet;
427  m[M22] = Det4_0134_0134 * oneOverDet;
428  m[M23] = Det4_0124_0134 * mn1OverDet;
429  m[M24] = Det4_0123_0134 * oneOverDet;
430 
431  m[M30] = Det4_1234_0124 * mn1OverDet;
432  m[M31] = Det4_0234_0124 * oneOverDet;
433  m[M32] = Det4_0134_0124 * mn1OverDet;
434  m[M33] = Det4_0124_0124 * oneOverDet;
435  m[M34] = Det4_0123_0124 * mn1OverDet;
436 
437  m[M40] = Det4_1234_0123 * oneOverDet;
438  m[M41] = Det4_0234_0123 * mn1OverDet;
439  m[M42] = Det4_0134_0123 * oneOverDet;
440  m[M43] = Det4_0124_0123 * mn1OverDet;
441  m[M44] = Det4_0123_0123 * oneOverDet;
442 
443  return;
444 }
445 
446 void HepMatrix::invertHaywood6 (int & ifail) {
447 
448  ifail = 0;
449 
450  // Find all NECESSARY 2x2 dets: (45 of them)
451 
452  double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
453  double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
454  double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
455  double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
456  double Det2_34_05 = m[A30]*m[A45] - m[A35]*m[A40]; //
457  double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
458  double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
459  double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
460  double Det2_34_15 = m[A31]*m[A45] - m[A35]*m[A41]; //
461  double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
462  double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
463  double Det2_34_25 = m[A32]*m[A45] - m[A35]*m[A42]; //
464  double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
465  double Det2_34_35 = m[A33]*m[A45] - m[A35]*m[A43]; //
466  double Det2_34_45 = m[A34]*m[A45] - m[A35]*m[A44]; //
467  double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
468  double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
469  double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
470  double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
471  double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
472  double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
473  double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
474  double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
475  double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
476  double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
477  double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
478  double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
479  double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
480  double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
481  double Det2_35_45 = m[A34]*m[A55] - m[A35]*m[A54]; //
482  double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
483  double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
484  double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
485  double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
486  double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
487  double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
488  double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
489  double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
490  double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
491  double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
492  double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
493  double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
494  double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
495  double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
496  double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
497 
498  // Find all NECESSARY 3x3 dets: (80 of them)
499 
500  double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
501  + m[A22]*Det2_34_01;
502  double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
503  + m[A23]*Det2_34_01;
504  double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
505  + m[A24]*Det2_34_01;
506  double Det3_234_015 = m[A20]*Det2_34_15 - m[A21]*Det2_34_05
507  + m[A25]*Det2_34_01; //
508  double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
509  + m[A23]*Det2_34_02;
510  double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
511  + m[A24]*Det2_34_02;
512  double Det3_234_025 = m[A20]*Det2_34_25 - m[A22]*Det2_34_05
513  + m[A25]*Det2_34_02; //
514  double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
515  + m[A24]*Det2_34_03;
516  double Det3_234_035 = m[A20]*Det2_34_35 - m[A23]*Det2_34_05
517  + m[A25]*Det2_34_03; //
518  double Det3_234_045 = m[A20]*Det2_34_45 - m[A24]*Det2_34_05
519  + m[A25]*Det2_34_04; //
520  double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
521  + m[A23]*Det2_34_12;
522  double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
523  + m[A24]*Det2_34_12;
524  double Det3_234_125 = m[A21]*Det2_34_25 - m[A22]*Det2_34_15
525  + m[A25]*Det2_34_12; //
526  double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
527  + m[A24]*Det2_34_13;
528  double Det3_234_135 = m[A21]*Det2_34_35 - m[A23]*Det2_34_15
529  + m[A25]*Det2_34_13; //
530  double Det3_234_145 = m[A21]*Det2_34_45 - m[A24]*Det2_34_15
531  + m[A25]*Det2_34_14; //
532  double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
533  + m[A24]*Det2_34_23;
534  double Det3_234_235 = m[A22]*Det2_34_35 - m[A23]*Det2_34_25
535  + m[A25]*Det2_34_23; //
536  double Det3_234_245 = m[A22]*Det2_34_45 - m[A24]*Det2_34_25
537  + m[A25]*Det2_34_24; //
538  double Det3_234_345 = m[A23]*Det2_34_45 - m[A24]*Det2_34_35
539  + m[A25]*Det2_34_34; //
540  double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02
541  + m[A22]*Det2_35_01;
542  double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03
543  + m[A23]*Det2_35_01;
544  double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04
545  + m[A24]*Det2_35_01;
546  double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05
547  + m[A25]*Det2_35_01;
548  double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03
549  + m[A23]*Det2_35_02;
550  double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04
551  + m[A24]*Det2_35_02;
552  double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05
553  + m[A25]*Det2_35_02;
554  double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04
555  + m[A24]*Det2_35_03;
556  double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05
557  + m[A25]*Det2_35_03;
558  double Det3_235_045 = m[A20]*Det2_35_45 - m[A24]*Det2_35_05
559  + m[A25]*Det2_35_04; //
560  double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13
561  + m[A23]*Det2_35_12;
562  double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14
563  + m[A24]*Det2_35_12;
564  double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15
565  + m[A25]*Det2_35_12;
566  double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14
567  + m[A24]*Det2_35_13;
568  double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15
569  + m[A25]*Det2_35_13;
570  double Det3_235_145 = m[A21]*Det2_35_45 - m[A24]*Det2_35_15
571  + m[A25]*Det2_35_14; //
572  double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24
573  + m[A24]*Det2_35_23;
574  double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25
575  + m[A25]*Det2_35_23;
576  double Det3_235_245 = m[A22]*Det2_35_45 - m[A24]*Det2_35_25
577  + m[A25]*Det2_35_24; //
578  double Det3_235_345 = m[A23]*Det2_35_45 - m[A24]*Det2_35_35
579  + m[A25]*Det2_35_34; //
580  double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02
581  + m[A22]*Det2_45_01;
582  double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03
583  + m[A23]*Det2_45_01;
584  double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04
585  + m[A24]*Det2_45_01;
586  double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05
587  + m[A25]*Det2_45_01;
588  double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03
589  + m[A23]*Det2_45_02;
590  double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04
591  + m[A24]*Det2_45_02;
592  double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05
593  + m[A25]*Det2_45_02;
594  double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04
595  + m[A24]*Det2_45_03;
596  double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05
597  + m[A25]*Det2_45_03;
598  double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05
599  + m[A25]*Det2_45_04;
600  double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13
601  + m[A23]*Det2_45_12;
602  double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14
603  + m[A24]*Det2_45_12;
604  double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15
605  + m[A25]*Det2_45_12;
606  double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14
607  + m[A24]*Det2_45_13;
608  double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15
609  + m[A25]*Det2_45_13;
610  double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15
611  + m[A25]*Det2_45_14;
612  double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24
613  + m[A24]*Det2_45_23;
614  double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25
615  + m[A25]*Det2_45_23;
616  double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25
617  + m[A25]*Det2_45_24;
618  double Det3_245_345 = m[A23]*Det2_45_45 - m[A24]*Det2_45_35
619  + m[A25]*Det2_45_34; //
620  double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02
621  + m[A32]*Det2_45_01;
622  double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03
623  + m[A33]*Det2_45_01;
624  double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04
625  + m[A34]*Det2_45_01;
626  double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05
627  + m[A35]*Det2_45_01;
628  double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03
629  + m[A33]*Det2_45_02;
630  double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04
631  + m[A34]*Det2_45_02;
632  double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05
633  + m[A35]*Det2_45_02;
634  double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04
635  + m[A34]*Det2_45_03;
636  double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05
637  + m[A35]*Det2_45_03;
638  double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05
639  + m[A35]*Det2_45_04;
640  double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13
641  + m[A33]*Det2_45_12;
642  double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14
643  + m[A34]*Det2_45_12;
644  double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15
645  + m[A35]*Det2_45_12;
646  double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14
647  + m[A34]*Det2_45_13;
648  double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15
649  + m[A35]*Det2_45_13;
650  double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15
651  + m[A35]*Det2_45_14;
652  double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24
653  + m[A34]*Det2_45_23;
654  double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25
655  + m[A35]*Det2_45_23;
656  double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25
657  + m[A35]*Det2_45_24;
658  double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35
659  + m[A35]*Det2_45_34;
660 
661  // Find all NECESSARY 4x4 dets: (75 of them)
662 
663  double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
664  + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
665  double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
666  + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
667  double Det4_1234_0125 = m[A10]*Det3_234_125 - m[A11]*Det3_234_025
668  + m[A12]*Det3_234_015 - m[A15]*Det3_234_012; //
669  double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
670  + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
671  double Det4_1234_0135 = m[A10]*Det3_234_135 - m[A11]*Det3_234_035
672  + m[A13]*Det3_234_015 - m[A15]*Det3_234_013; //
673  double Det4_1234_0145 = m[A10]*Det3_234_145 - m[A11]*Det3_234_045
674  + m[A14]*Det3_234_015 - m[A15]*Det3_234_014; //
675  double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
676  + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
677  double Det4_1234_0235 = m[A10]*Det3_234_235 - m[A12]*Det3_234_035
678  + m[A13]*Det3_234_025 - m[A15]*Det3_234_023; //
679  double Det4_1234_0245 = m[A10]*Det3_234_245 - m[A12]*Det3_234_045
680  + m[A14]*Det3_234_025 - m[A15]*Det3_234_024; //
681  double Det4_1234_0345 = m[A10]*Det3_234_345 - m[A13]*Det3_234_045
682  + m[A14]*Det3_234_035 - m[A15]*Det3_234_034; //
683  double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
684  + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
685  double Det4_1234_1235 = m[A11]*Det3_234_235 - m[A12]*Det3_234_135
686  + m[A13]*Det3_234_125 - m[A15]*Det3_234_123; //
687  double Det4_1234_1245 = m[A11]*Det3_234_245 - m[A12]*Det3_234_145
688  + m[A14]*Det3_234_125 - m[A15]*Det3_234_124; //
689  double Det4_1234_1345 = m[A11]*Det3_234_345 - m[A13]*Det3_234_145
690  + m[A14]*Det3_234_135 - m[A15]*Det3_234_134; //
691  double Det4_1234_2345 = m[A12]*Det3_234_345 - m[A13]*Det3_234_245
692  + m[A14]*Det3_234_235 - m[A15]*Det3_234_234; //
693  double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023
694  + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
695  double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024
696  + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
697  double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025
698  + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
699  double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034
700  + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
701  double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035
702  + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
703  double Det4_1235_0145 = m[A10]*Det3_235_145 - m[A11]*Det3_235_045
704  + m[A14]*Det3_235_015 - m[A15]*Det3_235_014; //
705  double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034
706  + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
707  double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035
708  + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
709  double Det4_1235_0245 = m[A10]*Det3_235_245 - m[A12]*Det3_235_045
710  + m[A14]*Det3_235_025 - m[A15]*Det3_235_024; //
711  double Det4_1235_0345 = m[A10]*Det3_235_345 - m[A13]*Det3_235_045
712  + m[A14]*Det3_235_035 - m[A15]*Det3_235_034; //
713  double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134
714  + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
715  double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135
716  + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
717  double Det4_1235_1245 = m[A11]*Det3_235_245 - m[A12]*Det3_235_145
718  + m[A14]*Det3_235_125 - m[A15]*Det3_235_124; //
719  double Det4_1235_1345 = m[A11]*Det3_235_345 - m[A13]*Det3_235_145
720  + m[A14]*Det3_235_135 - m[A15]*Det3_235_134; //
721  double Det4_1235_2345 = m[A12]*Det3_235_345 - m[A13]*Det3_235_245
722  + m[A14]*Det3_235_235 - m[A15]*Det3_235_234; //
723  double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023
724  + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
725  double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024
726  + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
727  double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025
728  + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
729  double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034
730  + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
731  double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035
732  + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
733  double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045
734  + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
735  double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034
736  + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
737  double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035
738  + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
739  double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045
740  + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
741  double Det4_1245_0345 = m[A10]*Det3_245_345 - m[A13]*Det3_245_045
742  + m[A14]*Det3_245_035 - m[A15]*Det3_245_034; //
743  double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134
744  + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
745  double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135
746  + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
747  double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145
748  + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
749  double Det4_1245_1345 = m[A11]*Det3_245_345 - m[A13]*Det3_245_145
750  + m[A14]*Det3_245_135 - m[A15]*Det3_245_134; //
751  double Det4_1245_2345 = m[A12]*Det3_245_345 - m[A13]*Det3_245_245
752  + m[A14]*Det3_245_235 - m[A15]*Det3_245_234; //
753  double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023
754  + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
755  double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024
756  + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
757  double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025
758  + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
759  double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034
760  + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
761  double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035
762  + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
763  double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045
764  + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
765  double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034
766  + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
767  double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035
768  + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
769  double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045
770  + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
771  double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045
772  + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
773  double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134
774  + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
775  double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135
776  + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
777  double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145
778  + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
779  double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145
780  + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
781  double Det4_1345_2345 = m[A12]*Det3_345_345 - m[A13]*Det3_345_245
782  + m[A14]*Det3_345_235 - m[A15]*Det3_345_234; //
783  double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023
784  + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
785  double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024
786  + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
787  double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025
788  + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
789  double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034
790  + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
791  double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035
792  + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
793  double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045
794  + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
795  double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034
796  + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
797  double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035
798  + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
799  double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045
800  + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
801  double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045
802  + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
803  double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134
804  + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
805  double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135
806  + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
807  double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145
808  + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
809  double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145
810  + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
811  double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245
812  + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
813 
814  // Find all NECESSARY 5x5 dets: (36 of them)
815 
816  double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234
817  + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
818  double Det5_01234_01235 = m[A00]*Det4_1234_1235 - m[A01]*Det4_1234_0235
819  + m[A02]*Det4_1234_0135 - m[A03]*Det4_1234_0125 + m[A05]*Det4_1234_0123; //
820  double Det5_01234_01245 = m[A00]*Det4_1234_1245 - m[A01]*Det4_1234_0245
821  + m[A02]*Det4_1234_0145 - m[A04]*Det4_1234_0125 + m[A05]*Det4_1234_0124; //
822  double Det5_01234_01345 = m[A00]*Det4_1234_1345 - m[A01]*Det4_1234_0345
823  + m[A03]*Det4_1234_0145 - m[A04]*Det4_1234_0135 + m[A05]*Det4_1234_0134; //
824  double Det5_01234_02345 = m[A00]*Det4_1234_2345 - m[A02]*Det4_1234_0345
825  + m[A03]*Det4_1234_0245 - m[A04]*Det4_1234_0235 + m[A05]*Det4_1234_0234; //
826  double Det5_01234_12345 = m[A01]*Det4_1234_2345 - m[A02]*Det4_1234_1345
827  + m[A03]*Det4_1234_1245 - m[A04]*Det4_1234_1235 + m[A05]*Det4_1234_1234; //
828  double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234
829  + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
830  double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235
831  + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
832  double Det5_01235_01245 = m[A00]*Det4_1235_1245 - m[A01]*Det4_1235_0245
833  + m[A02]*Det4_1235_0145 - m[A04]*Det4_1235_0125 + m[A05]*Det4_1235_0124; //
834  double Det5_01235_01345 = m[A00]*Det4_1235_1345 - m[A01]*Det4_1235_0345
835  + m[A03]*Det4_1235_0145 - m[A04]*Det4_1235_0135 + m[A05]*Det4_1235_0134; //
836  double Det5_01235_02345 = m[A00]*Det4_1235_2345 - m[A02]*Det4_1235_0345
837  + m[A03]*Det4_1235_0245 - m[A04]*Det4_1235_0235 + m[A05]*Det4_1235_0234; //
838  double Det5_01235_12345 = m[A01]*Det4_1235_2345 - m[A02]*Det4_1235_1345
839  + m[A03]*Det4_1235_1245 - m[A04]*Det4_1235_1235 + m[A05]*Det4_1235_1234; //
840  double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234
841  + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
842  double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235
843  + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
844  double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245
845  + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
846  double Det5_01245_01345 = m[A00]*Det4_1245_1345 - m[A01]*Det4_1245_0345
847  + m[A03]*Det4_1245_0145 - m[A04]*Det4_1245_0135 + m[A05]*Det4_1245_0134; //
848  double Det5_01245_02345 = m[A00]*Det4_1245_2345 - m[A02]*Det4_1245_0345
849  + m[A03]*Det4_1245_0245 - m[A04]*Det4_1245_0235 + m[A05]*Det4_1245_0234; //
850  double Det5_01245_12345 = m[A01]*Det4_1245_2345 - m[A02]*Det4_1245_1345
851  + m[A03]*Det4_1245_1245 - m[A04]*Det4_1245_1235 + m[A05]*Det4_1245_1234; //
852  double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234
853  + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
854  double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235
855  + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
856  double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245
857  + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
858  double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345
859  + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
860  double Det5_01345_02345 = m[A00]*Det4_1345_2345 - m[A02]*Det4_1345_0345
861  + m[A03]*Det4_1345_0245 - m[A04]*Det4_1345_0235 + m[A05]*Det4_1345_0234; //
862  double Det5_01345_12345 = m[A01]*Det4_1345_2345 - m[A02]*Det4_1345_1345
863  + m[A03]*Det4_1345_1245 - m[A04]*Det4_1345_1235 + m[A05]*Det4_1345_1234; //
864  double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234
865  + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
866  double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235
867  + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
868  double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245
869  + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
870  double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345
871  + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
872  double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345
873  + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
874  double Det5_02345_12345 = m[A01]*Det4_2345_2345 - m[A02]*Det4_2345_1345
875  + m[A03]*Det4_2345_1245 - m[A04]*Det4_2345_1235 + m[A05]*Det4_2345_1234; //
876  double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234
877  + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
878  double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235
879  + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
880  double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245
881  + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
882  double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345
883  + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
884  double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345
885  + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
886  double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345
887  + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
888 
889  // Find the determinant
890 
891  double det = m[A00]*Det5_12345_12345
892  - m[A01]*Det5_12345_02345
893  + m[A02]*Det5_12345_01345
894  - m[A03]*Det5_12345_01245
895  + m[A04]*Det5_12345_01235
896  - m[A05]*Det5_12345_01234;
897 
898  if ( det == 0 ) {
899 #ifdef SINGULAR_DIAGNOSTICS
900  std::cerr << "Kramer's rule inversion of a singular 6x6 matrix: "
901  << *this << "\n";
902 #endif
903  ifail = 1;
904  return;
905  }
906 
907  double oneOverDet = 1.0/det;
908  double mn1OverDet = - oneOverDet;
909 
910  m[A00] = Det5_12345_12345*oneOverDet;
911  m[A01] = Det5_02345_12345*mn1OverDet;
912  m[A02] = Det5_01345_12345*oneOverDet;
913  m[A03] = Det5_01245_12345*mn1OverDet;
914  m[A04] = Det5_01235_12345*oneOverDet;
915  m[A05] = Det5_01234_12345*mn1OverDet;
916 
917  m[A10] = Det5_12345_02345*mn1OverDet;
918  m[A11] = Det5_02345_02345*oneOverDet;
919  m[A12] = Det5_01345_02345*mn1OverDet;
920  m[A13] = Det5_01245_02345*oneOverDet;
921  m[A14] = Det5_01235_02345*mn1OverDet;
922  m[A15] = Det5_01234_02345*oneOverDet;
923 
924  m[A20] = Det5_12345_01345*oneOverDet;
925  m[A21] = Det5_02345_01345*mn1OverDet;
926  m[A22] = Det5_01345_01345*oneOverDet;
927  m[A23] = Det5_01245_01345*mn1OverDet;
928  m[A24] = Det5_01235_01345*oneOverDet;
929  m[A25] = Det5_01234_01345*mn1OverDet;
930 
931  m[A30] = Det5_12345_01245*mn1OverDet;
932  m[A31] = Det5_02345_01245*oneOverDet;
933  m[A32] = Det5_01345_01245*mn1OverDet;
934  m[A33] = Det5_01245_01245*oneOverDet;
935  m[A34] = Det5_01235_01245*mn1OverDet;
936  m[A35] = Det5_01234_01245*oneOverDet;
937 
938  m[A40] = Det5_12345_01235*oneOverDet;
939  m[A41] = Det5_02345_01235*mn1OverDet;
940  m[A42] = Det5_01345_01235*oneOverDet;
941  m[A43] = Det5_01245_01235*mn1OverDet;
942  m[A44] = Det5_01235_01235*oneOverDet;
943  m[A45] = Det5_01234_01235*mn1OverDet;
944 
945  m[A50] = Det5_12345_01234*mn1OverDet;
946  m[A51] = Det5_02345_01234*oneOverDet;
947  m[A52] = Det5_01345_01234*mn1OverDet;
948  m[A53] = Det5_01245_01234*oneOverDet;
949  m[A54] = Det5_01235_01234*mn1OverDet;
950  m[A55] = Det5_01234_01234*oneOverDet;
951 
952  return;
953 }
954 
955 
956 } // namespace CLHEP
Generated on Sun Jun 17 2012 08:08:26 for CLHEP by   doxygen 1.8.1.1