16 #include "CLHEP/Matrix/defs.h"
17 #include "CLHEP/Matrix/SymMatrix.h"
21 double HepSymMatrix::posDefFraction5x5 = 1.0;
22 double HepSymMatrix::posDefFraction6x6 = 1.0;
23 double HepSymMatrix::adjustment5x5 = 0.0;
24 double HepSymMatrix::adjustment6x6 = 0.0;
25 const double HepSymMatrix::CHOLESKY_THRESHOLD_5x5 = .5;
26 const double HepSymMatrix::CHOLESKY_THRESHOLD_6x6 = .2;
27 const double HepSymMatrix::CHOLESKY_CREEP_5x5 = .005;
28 const double HepSymMatrix::CHOLESKY_CREEP_6x6 = .002;
77 void HepSymMatrix::invert5(
int & ifail) {
78 if (posDefFraction5x5 >= CHOLESKY_THRESHOLD_5x5) {
80 posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
85 if (posDefFraction5x5 + adjustment5x5 >= CHOLESKY_THRESHOLD_5x5) {
87 posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
94 adjustment5x5 += CHOLESKY_CREEP_5x5;
100 void HepSymMatrix::invert6(
int & ifail) {
101 if (posDefFraction6x6 >= CHOLESKY_THRESHOLD_6x6) {
103 posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
108 if (posDefFraction6x6 + adjustment6x6 >= CHOLESKY_THRESHOLD_6x6) {
110 posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
117 adjustment6x6 += CHOLESKY_CREEP_6x6;
158 double Det3_123_012 = m[
A10]*Det2_23_12 - m[
A11]*Det2_23_02
160 double Det3_123_013 = m[
A10]*Det2_23_13 - m[
A11]*Det2_23_03
162 double Det3_123_023 = m[
A10]*Det2_23_23 - m[
A12]*Det2_23_03
164 double Det3_123_123 = m[
A11]*Det2_23_23 - m[
A12]*Det2_23_13
166 double Det3_124_012 = m[
A10]*Det2_24_12 - m[
A11]*Det2_24_02
168 double Det3_124_013 = m[
A10]*Det2_24_13 - m[
A11]*Det2_24_03
170 double Det3_124_014 = m[
A10]*Det2_24_14 - m[
A11]*Det2_24_04
172 double Det3_124_023 = m[
A10]*Det2_24_23 - m[
A12]*Det2_24_03
174 double Det3_124_024 = m[
A10]*Det2_24_24 - m[
A12]*Det2_24_04
176 double Det3_124_123 = m[
A11]*Det2_24_23 - m[
A12]*Det2_24_13
178 double Det3_124_124 = m[
A11]*Det2_24_24 - m[
A12]*Det2_24_14
180 double Det3_134_012 = m[
A10]*Det2_34_12 - m[
A11]*Det2_34_02
182 double Det3_134_013 = m[
A10]*Det2_34_13 - m[
A11]*Det2_34_03
184 double Det3_134_014 = m[
A10]*Det2_34_14 - m[
A11]*Det2_34_04
186 double Det3_134_023 = m[
A10]*Det2_34_23 - m[
A12]*Det2_34_03
188 double Det3_134_024 = m[
A10]*Det2_34_24 - m[
A12]*Det2_34_04
190 double Det3_134_034 = m[
A10]*Det2_34_34 - m[
A13]*Det2_34_04
192 double Det3_134_123 = m[
A11]*Det2_34_23 - m[
A12]*Det2_34_13
194 double Det3_134_124 = m[
A11]*Det2_34_24 - m[
A12]*Det2_34_14
196 double Det3_134_134 = m[
A11]*Det2_34_34 - m[
A13]*Det2_34_14
198 double Det3_234_012 = m[
A20]*Det2_34_12 - m[
A21]*Det2_34_02
200 double Det3_234_013 = m[
A20]*Det2_34_13 - m[
A21]*Det2_34_03
202 double Det3_234_014 = m[
A20]*Det2_34_14 - m[
A21]*Det2_34_04
204 double Det3_234_023 = m[
A20]*Det2_34_23 - m[
A22]*Det2_34_03
206 double Det3_234_024 = m[
A20]*Det2_34_24 - m[
A22]*Det2_34_04
208 double Det3_234_034 = m[
A20]*Det2_34_34 - m[
A23]*Det2_34_04
210 double Det3_234_123 = m[
A21]*Det2_34_23 - m[
A22]*Det2_34_13
212 double Det3_234_124 = m[
A21]*Det2_34_24 - m[
A22]*Det2_34_14
214 double Det3_234_134 = m[
A21]*Det2_34_34 - m[
A23]*Det2_34_14
216 double Det3_234_234 = m[
A22]*Det2_34_34 - m[
A23]*Det2_34_24
221 double Det4_0123_0123 = m[
A00]*Det3_123_123 - m[
A01]*Det3_123_023
222 + m[
A02]*Det3_123_013 - m[
A03]*Det3_123_012;
223 double Det4_0124_0123 = m[
A00]*Det3_124_123 - m[
A01]*Det3_124_023
224 + m[
A02]*Det3_124_013 - m[
A03]*Det3_124_012;
225 double Det4_0124_0124 = m[
A00]*Det3_124_124 - m[
A01]*Det3_124_024
226 + m[
A02]*Det3_124_014 - m[
A04]*Det3_124_012;
227 double Det4_0134_0123 = m[
A00]*Det3_134_123 - m[
A01]*Det3_134_023
228 + m[
A02]*Det3_134_013 - m[
A03]*Det3_134_012;
229 double Det4_0134_0124 = m[
A00]*Det3_134_124 - m[
A01]*Det3_134_024
230 + m[
A02]*Det3_134_014 - m[
A04]*Det3_134_012;
231 double Det4_0134_0134 = m[
A00]*Det3_134_134 - m[
A01]*Det3_134_034
232 + m[
A03]*Det3_134_014 - m[
A04]*Det3_134_013;
233 double Det4_0234_0123 = m[
A00]*Det3_234_123 - m[
A01]*Det3_234_023
234 + m[
A02]*Det3_234_013 - m[
A03]*Det3_234_012;
235 double Det4_0234_0124 = m[
A00]*Det3_234_124 - m[
A01]*Det3_234_024
236 + m[
A02]*Det3_234_014 - m[
A04]*Det3_234_012;
237 double Det4_0234_0134 = m[
A00]*Det3_234_134 - m[
A01]*Det3_234_034
238 + m[
A03]*Det3_234_014 - m[
A04]*Det3_234_013;
239 double Det4_0234_0234 = m[
A00]*Det3_234_234 - m[
A02]*Det3_234_034
240 + m[
A03]*Det3_234_024 - m[
A04]*Det3_234_023;
241 double Det4_1234_0123 = m[
A10]*Det3_234_123 - m[
A11]*Det3_234_023
242 + m[
A12]*Det3_234_013 - m[
A13]*Det3_234_012;
243 double Det4_1234_0124 = m[
A10]*Det3_234_124 - m[
A11]*Det3_234_024
244 + m[
A12]*Det3_234_014 - m[
A14]*Det3_234_012;
245 double Det4_1234_0134 = m[
A10]*Det3_234_134 - m[
A11]*Det3_234_034
246 + m[
A13]*Det3_234_014 - m[
A14]*Det3_234_013;
247 double Det4_1234_0234 = m[
A10]*Det3_234_234 - m[
A12]*Det3_234_034
248 + m[
A13]*Det3_234_024 - m[
A14]*Det3_234_023;
249 double Det4_1234_1234 = m[
A11]*Det3_234_234 - m[
A12]*Det3_234_134
250 + m[
A13]*Det3_234_124 - m[
A14]*Det3_234_123;
254 double det = m[
A00]*Det4_1234_1234
255 - m[
A01]*Det4_1234_0234
256 + m[
A02]*Det4_1234_0134
257 - m[
A03]*Det4_1234_0124
258 + m[
A04]*Det4_1234_0123;
261 #ifdef SINGULAR_DIAGNOSTICS
262 std::cerr <<
"Kramer's rule inversion of a singular 5x5 matrix: "
269 double oneOverDet = 1.0/det;
270 double mn1OverDet = - oneOverDet;
272 m[
A00] = Det4_1234_1234 * oneOverDet;
273 m[
A01] = Det4_1234_0234 * mn1OverDet;
274 m[
A02] = Det4_1234_0134 * oneOverDet;
275 m[
A03] = Det4_1234_0124 * mn1OverDet;
276 m[
A04] = Det4_1234_0123 * oneOverDet;
278 m[
A11] = Det4_0234_0234 * oneOverDet;
279 m[
A12] = Det4_0234_0134 * mn1OverDet;
280 m[
A13] = Det4_0234_0124 * oneOverDet;
281 m[
A14] = Det4_0234_0123 * mn1OverDet;
283 m[
A22] = Det4_0134_0134 * oneOverDet;
284 m[
A23] = Det4_0134_0124 * mn1OverDet;
285 m[
A24] = Det4_0134_0123 * oneOverDet;
287 m[
A33] = Det4_0124_0124 * oneOverDet;
288 m[
A34] = Det4_0124_0123 * mn1OverDet;
290 m[
A44] = Det4_0123_0123 * oneOverDet;
343 double Det3_234_012 = m[
A20]*Det2_34_12 - m[
A21]*Det2_34_02
345 double Det3_234_013 = m[
A20]*Det2_34_13 - m[
A21]*Det2_34_03
347 double Det3_234_014 = m[
A20]*Det2_34_14 - m[
A21]*Det2_34_04
349 double Det3_234_023 = m[
A20]*Det2_34_23 - m[
A22]*Det2_34_03
351 double Det3_234_024 = m[
A20]*Det2_34_24 - m[
A22]*Det2_34_04
353 double Det3_234_034 = m[
A20]*Det2_34_34 - m[
A23]*Det2_34_04
355 double Det3_234_123 = m[
A21]*Det2_34_23 - m[
A22]*Det2_34_13
357 double Det3_234_124 = m[
A21]*Det2_34_24 - m[
A22]*Det2_34_14
359 double Det3_234_134 = m[
A21]*Det2_34_34 - m[
A23]*Det2_34_14
361 double Det3_234_234 = m[
A22]*Det2_34_34 - m[
A23]*Det2_34_24
363 double Det3_235_012 = m[
A20]*Det2_35_12 - m[
A21]*Det2_35_02
365 double Det3_235_013 = m[
A20]*Det2_35_13 - m[
A21]*Det2_35_03
367 double Det3_235_014 = m[
A20]*Det2_35_14 - m[
A21]*Det2_35_04
369 double Det3_235_015 = m[
A20]*Det2_35_15 - m[
A21]*Det2_35_05
371 double Det3_235_023 = m[
A20]*Det2_35_23 - m[
A22]*Det2_35_03
373 double Det3_235_024 = m[
A20]*Det2_35_24 - m[
A22]*Det2_35_04
375 double Det3_235_025 = m[
A20]*Det2_35_25 - m[
A22]*Det2_35_05
377 double Det3_235_034 = m[
A20]*Det2_35_34 - m[
A23]*Det2_35_04
379 double Det3_235_035 = m[
A20]*Det2_35_35 - m[
A23]*Det2_35_05
381 double Det3_235_123 = m[
A21]*Det2_35_23 - m[
A22]*Det2_35_13
383 double Det3_235_124 = m[
A21]*Det2_35_24 - m[
A22]*Det2_35_14
385 double Det3_235_125 = m[
A21]*Det2_35_25 - m[
A22]*Det2_35_15
387 double Det3_235_134 = m[
A21]*Det2_35_34 - m[
A23]*Det2_35_14
389 double Det3_235_135 = m[
A21]*Det2_35_35 - m[
A23]*Det2_35_15
391 double Det3_235_234 = m[
A22]*Det2_35_34 - m[
A23]*Det2_35_24
393 double Det3_235_235 = m[
A22]*Det2_35_35 - m[
A23]*Det2_35_25
395 double Det3_245_012 = m[
A20]*Det2_45_12 - m[
A21]*Det2_45_02
397 double Det3_245_013 = m[
A20]*Det2_45_13 - m[
A21]*Det2_45_03
399 double Det3_245_014 = m[
A20]*Det2_45_14 - m[
A21]*Det2_45_04
401 double Det3_245_015 = m[
A20]*Det2_45_15 - m[
A21]*Det2_45_05
403 double Det3_245_023 = m[
A20]*Det2_45_23 - m[
A22]*Det2_45_03
405 double Det3_245_024 = m[
A20]*Det2_45_24 - m[
A22]*Det2_45_04
407 double Det3_245_025 = m[
A20]*Det2_45_25 - m[
A22]*Det2_45_05
409 double Det3_245_034 = m[
A20]*Det2_45_34 - m[
A23]*Det2_45_04
411 double Det3_245_035 = m[
A20]*Det2_45_35 - m[
A23]*Det2_45_05
413 double Det3_245_045 = m[
A20]*Det2_45_45 - m[
A24]*Det2_45_05
415 double Det3_245_123 = m[
A21]*Det2_45_23 - m[
A22]*Det2_45_13
417 double Det3_245_124 = m[
A21]*Det2_45_24 - m[
A22]*Det2_45_14
419 double Det3_245_125 = m[
A21]*Det2_45_25 - m[
A22]*Det2_45_15
421 double Det3_245_134 = m[
A21]*Det2_45_34 - m[
A23]*Det2_45_14
423 double Det3_245_135 = m[
A21]*Det2_45_35 - m[
A23]*Det2_45_15
425 double Det3_245_145 = m[
A21]*Det2_45_45 - m[
A24]*Det2_45_15
427 double Det3_245_234 = m[
A22]*Det2_45_34 - m[
A23]*Det2_45_24
429 double Det3_245_235 = m[
A22]*Det2_45_35 - m[
A23]*Det2_45_25
431 double Det3_245_245 = m[
A22]*Det2_45_45 - m[
A24]*Det2_45_25
433 double Det3_345_012 = m[
A30]*Det2_45_12 - m[
A31]*Det2_45_02
435 double Det3_345_013 = m[
A30]*Det2_45_13 - m[
A31]*Det2_45_03
437 double Det3_345_014 = m[
A30]*Det2_45_14 - m[
A31]*Det2_45_04
439 double Det3_345_015 = m[
A30]*Det2_45_15 - m[
A31]*Det2_45_05
441 double Det3_345_023 = m[
A30]*Det2_45_23 - m[
A32]*Det2_45_03
443 double Det3_345_024 = m[
A30]*Det2_45_24 - m[
A32]*Det2_45_04
445 double Det3_345_025 = m[
A30]*Det2_45_25 - m[
A32]*Det2_45_05
447 double Det3_345_034 = m[
A30]*Det2_45_34 - m[
A33]*Det2_45_04
449 double Det3_345_035 = m[
A30]*Det2_45_35 - m[
A33]*Det2_45_05
451 double Det3_345_045 = m[
A30]*Det2_45_45 - m[
A34]*Det2_45_05
453 double Det3_345_123 = m[
A31]*Det2_45_23 - m[
A32]*Det2_45_13
455 double Det3_345_124 = m[
A31]*Det2_45_24 - m[
A32]*Det2_45_14
457 double Det3_345_125 = m[
A31]*Det2_45_25 - m[
A32]*Det2_45_15
459 double Det3_345_134 = m[
A31]*Det2_45_34 - m[
A33]*Det2_45_14
461 double Det3_345_135 = m[
A31]*Det2_45_35 - m[
A33]*Det2_45_15
463 double Det3_345_145 = m[
A31]*Det2_45_45 - m[
A34]*Det2_45_15
465 double Det3_345_234 = m[
A32]*Det2_45_34 - m[
A33]*Det2_45_24
467 double Det3_345_235 = m[
A32]*Det2_45_35 - m[
A33]*Det2_45_25
469 double Det3_345_245 = m[
A32]*Det2_45_45 - m[
A34]*Det2_45_25
471 double Det3_345_345 = m[
A33]*Det2_45_45 - m[
A34]*Det2_45_35
476 double Det4_1234_0123 = m[
A10]*Det3_234_123 - m[
A11]*Det3_234_023
477 + m[
A12]*Det3_234_013 - m[
A13]*Det3_234_012;
478 double Det4_1234_0124 = m[
A10]*Det3_234_124 - m[
A11]*Det3_234_024
479 + m[
A12]*Det3_234_014 - m[
A14]*Det3_234_012;
480 double Det4_1234_0134 = m[
A10]*Det3_234_134 - m[
A11]*Det3_234_034
481 + m[
A13]*Det3_234_014 - m[
A14]*Det3_234_013;
482 double Det4_1234_0234 = m[
A10]*Det3_234_234 - m[
A12]*Det3_234_034
483 + m[
A13]*Det3_234_024 - m[
A14]*Det3_234_023;
484 double Det4_1234_1234 = m[
A11]*Det3_234_234 - m[
A12]*Det3_234_134
485 + m[
A13]*Det3_234_124 - m[
A14]*Det3_234_123;
486 double Det4_1235_0123 = m[
A10]*Det3_235_123 - m[
A11]*Det3_235_023
487 + m[
A12]*Det3_235_013 - m[
A13]*Det3_235_012;
488 double Det4_1235_0124 = m[
A10]*Det3_235_124 - m[
A11]*Det3_235_024
489 + m[
A12]*Det3_235_014 - m[
A14]*Det3_235_012;
490 double Det4_1235_0125 = m[
A10]*Det3_235_125 - m[
A11]*Det3_235_025
491 + m[
A12]*Det3_235_015 - m[
A15]*Det3_235_012;
492 double Det4_1235_0134 = m[
A10]*Det3_235_134 - m[
A11]*Det3_235_034
493 + m[
A13]*Det3_235_014 - m[
A14]*Det3_235_013;
494 double Det4_1235_0135 = m[
A10]*Det3_235_135 - m[
A11]*Det3_235_035
495 + m[
A13]*Det3_235_015 - m[
A15]*Det3_235_013;
496 double Det4_1235_0234 = m[
A10]*Det3_235_234 - m[
A12]*Det3_235_034
497 + m[
A13]*Det3_235_024 - m[
A14]*Det3_235_023;
498 double Det4_1235_0235 = m[
A10]*Det3_235_235 - m[
A12]*Det3_235_035
499 + m[
A13]*Det3_235_025 - m[
A15]*Det3_235_023;
500 double Det4_1235_1234 = m[
A11]*Det3_235_234 - m[
A12]*Det3_235_134
501 + m[
A13]*Det3_235_124 - m[
A14]*Det3_235_123;
502 double Det4_1235_1235 = m[
A11]*Det3_235_235 - m[
A12]*Det3_235_135
503 + m[
A13]*Det3_235_125 - m[
A15]*Det3_235_123;
504 double Det4_1245_0123 = m[
A10]*Det3_245_123 - m[
A11]*Det3_245_023
505 + m[
A12]*Det3_245_013 - m[
A13]*Det3_245_012;
506 double Det4_1245_0124 = m[
A10]*Det3_245_124 - m[
A11]*Det3_245_024
507 + m[
A12]*Det3_245_014 - m[
A14]*Det3_245_012;
508 double Det4_1245_0125 = m[
A10]*Det3_245_125 - m[
A11]*Det3_245_025
509 + m[
A12]*Det3_245_015 - m[
A15]*Det3_245_012;
510 double Det4_1245_0134 = m[
A10]*Det3_245_134 - m[
A11]*Det3_245_034
511 + m[
A13]*Det3_245_014 - m[
A14]*Det3_245_013;
512 double Det4_1245_0135 = m[
A10]*Det3_245_135 - m[
A11]*Det3_245_035
513 + m[
A13]*Det3_245_015 - m[
A15]*Det3_245_013;
514 double Det4_1245_0145 = m[
A10]*Det3_245_145 - m[
A11]*Det3_245_045
515 + m[
A14]*Det3_245_015 - m[
A15]*Det3_245_014;
516 double Det4_1245_0234 = m[
A10]*Det3_245_234 - m[
A12]*Det3_245_034
517 + m[
A13]*Det3_245_024 - m[
A14]*Det3_245_023;
518 double Det4_1245_0235 = m[
A10]*Det3_245_235 - m[
A12]*Det3_245_035
519 + m[
A13]*Det3_245_025 - m[
A15]*Det3_245_023;
520 double Det4_1245_0245 = m[
A10]*Det3_245_245 - m[
A12]*Det3_245_045
521 + m[
A14]*Det3_245_025 - m[
A15]*Det3_245_024;
522 double Det4_1245_1234 = m[
A11]*Det3_245_234 - m[
A12]*Det3_245_134
523 + m[
A13]*Det3_245_124 - m[
A14]*Det3_245_123;
524 double Det4_1245_1235 = m[
A11]*Det3_245_235 - m[
A12]*Det3_245_135
525 + m[
A13]*Det3_245_125 - m[
A15]*Det3_245_123;
526 double Det4_1245_1245 = m[
A11]*Det3_245_245 - m[
A12]*Det3_245_145
527 + m[
A14]*Det3_245_125 - m[
A15]*Det3_245_124;
528 double Det4_1345_0123 = m[
A10]*Det3_345_123 - m[
A11]*Det3_345_023
529 + m[
A12]*Det3_345_013 - m[
A13]*Det3_345_012;
530 double Det4_1345_0124 = m[
A10]*Det3_345_124 - m[
A11]*Det3_345_024
531 + m[
A12]*Det3_345_014 - m[
A14]*Det3_345_012;
532 double Det4_1345_0125 = m[
A10]*Det3_345_125 - m[
A11]*Det3_345_025
533 + m[
A12]*Det3_345_015 - m[
A15]*Det3_345_012;
534 double Det4_1345_0134 = m[
A10]*Det3_345_134 - m[
A11]*Det3_345_034
535 + m[
A13]*Det3_345_014 - m[
A14]*Det3_345_013;
536 double Det4_1345_0135 = m[
A10]*Det3_345_135 - m[
A11]*Det3_345_035
537 + m[
A13]*Det3_345_015 - m[
A15]*Det3_345_013;
538 double Det4_1345_0145 = m[
A10]*Det3_345_145 - m[
A11]*Det3_345_045
539 + m[
A14]*Det3_345_015 - m[
A15]*Det3_345_014;
540 double Det4_1345_0234 = m[
A10]*Det3_345_234 - m[
A12]*Det3_345_034
541 + m[
A13]*Det3_345_024 - m[
A14]*Det3_345_023;
542 double Det4_1345_0235 = m[
A10]*Det3_345_235 - m[
A12]*Det3_345_035
543 + m[
A13]*Det3_345_025 - m[
A15]*Det3_345_023;
544 double Det4_1345_0245 = m[
A10]*Det3_345_245 - m[
A12]*Det3_345_045
545 + m[
A14]*Det3_345_025 - m[
A15]*Det3_345_024;
546 double Det4_1345_0345 = m[
A10]*Det3_345_345 - m[
A13]*Det3_345_045
547 + m[
A14]*Det3_345_035 - m[
A15]*Det3_345_034;
548 double Det4_1345_1234 = m[
A11]*Det3_345_234 - m[
A12]*Det3_345_134
549 + m[
A13]*Det3_345_124 - m[
A14]*Det3_345_123;
550 double Det4_1345_1235 = m[
A11]*Det3_345_235 - m[
A12]*Det3_345_135
551 + m[
A13]*Det3_345_125 - m[
A15]*Det3_345_123;
552 double Det4_1345_1245 = m[
A11]*Det3_345_245 - m[
A12]*Det3_345_145
553 + m[
A14]*Det3_345_125 - m[
A15]*Det3_345_124;
554 double Det4_1345_1345 = m[
A11]*Det3_345_345 - m[
A13]*Det3_345_145
555 + m[
A14]*Det3_345_135 - m[
A15]*Det3_345_134;
556 double Det4_2345_0123 = m[
A20]*Det3_345_123 - m[
A21]*Det3_345_023
557 + m[
A22]*Det3_345_013 - m[
A23]*Det3_345_012;
558 double Det4_2345_0124 = m[
A20]*Det3_345_124 - m[
A21]*Det3_345_024
559 + m[
A22]*Det3_345_014 - m[
A24]*Det3_345_012;
560 double Det4_2345_0125 = m[
A20]*Det3_345_125 - m[
A21]*Det3_345_025
561 + m[
A22]*Det3_345_015 - m[
A25]*Det3_345_012;
562 double Det4_2345_0134 = m[
A20]*Det3_345_134 - m[
A21]*Det3_345_034
563 + m[
A23]*Det3_345_014 - m[
A24]*Det3_345_013;
564 double Det4_2345_0135 = m[
A20]*Det3_345_135 - m[
A21]*Det3_345_035
565 + m[
A23]*Det3_345_015 - m[
A25]*Det3_345_013;
566 double Det4_2345_0145 = m[
A20]*Det3_345_145 - m[
A21]*Det3_345_045
567 + m[
A24]*Det3_345_015 - m[
A25]*Det3_345_014;
568 double Det4_2345_0234 = m[
A20]*Det3_345_234 - m[
A22]*Det3_345_034
569 + m[
A23]*Det3_345_024 - m[
A24]*Det3_345_023;
570 double Det4_2345_0235 = m[
A20]*Det3_345_235 - m[
A22]*Det3_345_035
571 + m[
A23]*Det3_345_025 - m[
A25]*Det3_345_023;
572 double Det4_2345_0245 = m[
A20]*Det3_345_245 - m[
A22]*Det3_345_045
573 + m[
A24]*Det3_345_025 - m[
A25]*Det3_345_024;
574 double Det4_2345_0345 = m[
A20]*Det3_345_345 - m[
A23]*Det3_345_045
575 + m[
A24]*Det3_345_035 - m[
A25]*Det3_345_034;
576 double Det4_2345_1234 = m[
A21]*Det3_345_234 - m[
A22]*Det3_345_134
577 + m[
A23]*Det3_345_124 - m[
A24]*Det3_345_123;
578 double Det4_2345_1235 = m[
A21]*Det3_345_235 - m[
A22]*Det3_345_135
579 + m[
A23]*Det3_345_125 - m[
A25]*Det3_345_123;
580 double Det4_2345_1245 = m[
A21]*Det3_345_245 - m[
A22]*Det3_345_145
581 + m[
A24]*Det3_345_125 - m[
A25]*Det3_345_124;
582 double Det4_2345_1345 = m[
A21]*Det3_345_345 - m[
A23]*Det3_345_145
583 + m[
A24]*Det3_345_135 - m[
A25]*Det3_345_134;
584 double Det4_2345_2345 = m[
A22]*Det3_345_345 - m[
A23]*Det3_345_245
585 + m[
A24]*Det3_345_235 - m[
A25]*Det3_345_234;
589 double Det5_01234_01234 = m[
A00]*Det4_1234_1234 - m[
A01]*Det4_1234_0234
590 + m[
A02]*Det4_1234_0134 - m[
A03]*Det4_1234_0124 + m[
A04]*Det4_1234_0123;
591 double Det5_01235_01234 = m[
A00]*Det4_1235_1234 - m[
A01]*Det4_1235_0234
592 + m[
A02]*Det4_1235_0134 - m[
A03]*Det4_1235_0124 + m[
A04]*Det4_1235_0123;
593 double Det5_01235_01235 = m[
A00]*Det4_1235_1235 - m[
A01]*Det4_1235_0235
594 + m[
A02]*Det4_1235_0135 - m[
A03]*Det4_1235_0125 + m[
A05]*Det4_1235_0123;
595 double Det5_01245_01234 = m[
A00]*Det4_1245_1234 - m[
A01]*Det4_1245_0234
596 + m[
A02]*Det4_1245_0134 - m[
A03]*Det4_1245_0124 + m[
A04]*Det4_1245_0123;
597 double Det5_01245_01235 = m[
A00]*Det4_1245_1235 - m[
A01]*Det4_1245_0235
598 + m[
A02]*Det4_1245_0135 - m[
A03]*Det4_1245_0125 + m[
A05]*Det4_1245_0123;
599 double Det5_01245_01245 = m[
A00]*Det4_1245_1245 - m[
A01]*Det4_1245_0245
600 + m[
A02]*Det4_1245_0145 - m[
A04]*Det4_1245_0125 + m[
A05]*Det4_1245_0124;
601 double Det5_01345_01234 = m[
A00]*Det4_1345_1234 - m[
A01]*Det4_1345_0234
602 + m[
A02]*Det4_1345_0134 - m[
A03]*Det4_1345_0124 + m[
A04]*Det4_1345_0123;
603 double Det5_01345_01235 = m[
A00]*Det4_1345_1235 - m[
A01]*Det4_1345_0235
604 + m[
A02]*Det4_1345_0135 - m[
A03]*Det4_1345_0125 + m[
A05]*Det4_1345_0123;
605 double Det5_01345_01245 = m[
A00]*Det4_1345_1245 - m[
A01]*Det4_1345_0245
606 + m[
A02]*Det4_1345_0145 - m[
A04]*Det4_1345_0125 + m[
A05]*Det4_1345_0124;
607 double Det5_01345_01345 = m[
A00]*Det4_1345_1345 - m[
A01]*Det4_1345_0345
608 + m[
A03]*Det4_1345_0145 - m[
A04]*Det4_1345_0135 + m[
A05]*Det4_1345_0134;
609 double Det5_02345_01234 = m[
A00]*Det4_2345_1234 - m[
A01]*Det4_2345_0234
610 + m[
A02]*Det4_2345_0134 - m[
A03]*Det4_2345_0124 + m[
A04]*Det4_2345_0123;
611 double Det5_02345_01235 = m[
A00]*Det4_2345_1235 - m[
A01]*Det4_2345_0235
612 + m[
A02]*Det4_2345_0135 - m[
A03]*Det4_2345_0125 + m[
A05]*Det4_2345_0123;
613 double Det5_02345_01245 = m[
A00]*Det4_2345_1245 - m[
A01]*Det4_2345_0245
614 + m[
A02]*Det4_2345_0145 - m[
A04]*Det4_2345_0125 + m[
A05]*Det4_2345_0124;
615 double Det5_02345_01345 = m[
A00]*Det4_2345_1345 - m[
A01]*Det4_2345_0345
616 + m[
A03]*Det4_2345_0145 - m[
A04]*Det4_2345_0135 + m[
A05]*Det4_2345_0134;
617 double Det5_02345_02345 = m[
A00]*Det4_2345_2345 - m[
A02]*Det4_2345_0345
618 + m[
A03]*Det4_2345_0245 - m[
A04]*Det4_2345_0235 + m[
A05]*Det4_2345_0234;
619 double Det5_12345_01234 = m[
A10]*Det4_2345_1234 - m[
A11]*Det4_2345_0234
620 + m[
A12]*Det4_2345_0134 - m[
A13]*Det4_2345_0124 + m[
A14]*Det4_2345_0123;
621 double Det5_12345_01235 = m[
A10]*Det4_2345_1235 - m[
A11]*Det4_2345_0235
622 + m[
A12]*Det4_2345_0135 - m[
A13]*Det4_2345_0125 + m[
A15]*Det4_2345_0123;
623 double Det5_12345_01245 = m[
A10]*Det4_2345_1245 - m[
A11]*Det4_2345_0245
624 + m[
A12]*Det4_2345_0145 - m[
A14]*Det4_2345_0125 + m[
A15]*Det4_2345_0124;
625 double Det5_12345_01345 = m[
A10]*Det4_2345_1345 - m[
A11]*Det4_2345_0345
626 + m[
A13]*Det4_2345_0145 - m[
A14]*Det4_2345_0135 + m[
A15]*Det4_2345_0134;
627 double Det5_12345_02345 = m[
A10]*Det4_2345_2345 - m[
A12]*Det4_2345_0345
628 + m[
A13]*Det4_2345_0245 - m[
A14]*Det4_2345_0235 + m[
A15]*Det4_2345_0234;
629 double Det5_12345_12345 = m[
A11]*Det4_2345_2345 - m[
A12]*Det4_2345_1345
630 + m[
A13]*Det4_2345_1245 - m[
A14]*Det4_2345_1235 + m[
A15]*Det4_2345_1234;
634 double det = m[
A00]*Det5_12345_12345
635 - m[
A01]*Det5_12345_02345
636 + m[
A02]*Det5_12345_01345
637 - m[
A03]*Det5_12345_01245
638 + m[
A04]*Det5_12345_01235
639 - m[
A05]*Det5_12345_01234;
642 #ifdef SINGULAR_DIAGNOSTICS
643 std::cerr <<
"Kramer's rule inversion of a singular 6x6 matrix: "
650 double oneOverDet = 1.0/det;
651 double mn1OverDet = - oneOverDet;
653 m[
A00] = Det5_12345_12345*oneOverDet;
654 m[
A01] = Det5_12345_02345*mn1OverDet;
655 m[
A02] = Det5_12345_01345*oneOverDet;
656 m[
A03] = Det5_12345_01245*mn1OverDet;
657 m[
A04] = Det5_12345_01235*oneOverDet;
658 m[
A05] = Det5_12345_01234*mn1OverDet;
660 m[
A11] = Det5_02345_02345*oneOverDet;
661 m[
A12] = Det5_02345_01345*mn1OverDet;
662 m[
A13] = Det5_02345_01245*oneOverDet;
663 m[
A14] = Det5_02345_01235*mn1OverDet;
664 m[
A15] = Det5_02345_01234*oneOverDet;
666 m[
A22] = Det5_01345_01345*oneOverDet;
667 m[
A23] = Det5_01345_01245*mn1OverDet;
668 m[
A24] = Det5_01345_01235*oneOverDet;
669 m[
A25] = Det5_01345_01234*mn1OverDet;
671 m[
A33] = Det5_01245_01245*oneOverDet;
672 m[
A34] = Det5_01245_01235*mn1OverDet;
673 m[
A35] = Det5_01245_01234*oneOverDet;
675 m[
A44] = Det5_01235_01235*oneOverDet;
676 m[
A45] = Det5_01235_01234*mn1OverDet;
678 m[
A55] = Det5_01234_01234*oneOverDet;
697 double h30, h31, h32;
698 double h40, h41, h42, h43;
700 double h00, h11, h22, h33, h44;
705 double g30, g31, g32;
706 double g40, g41, g42, g43;
716 if (h00 <= 0)
return;
717 h00 = 1.0 / sqrt(h00);
726 h11 = m[
A11] - (g10 * g10);
727 if (h11 <= 0)
return;
728 h11 = 1.0 / sqrt(h11);
733 g21 = (m[
A21] - (g10 * g20)) * h11;
734 g31 = (m[
A31] - (g10 * g30)) * h11;
735 g41 = (m[
A41] - (g10 * g40)) * h11;
739 h22 = m[
A22] - (g20 * g20) - (g21 * g21);
740 if (h22 <= 0)
return;
741 h22 = 1.0 / sqrt(h22);
746 g32 = (m[
A32] - (g20 * g30) - (g21 * g31)) * h22;
747 g42 = (m[
A42] - (g20 * g40) - (g21 * g41)) * h22;
751 h33 = m[
A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
752 if (h33 <= 0)
return;
753 h33 = 1.0 / sqrt(h33);
757 g43 = (m[
A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
761 h44 = m[
A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
762 if (h44 <= 0)
return;
763 h44 = 1.0 / sqrt(h44);
770 h43 = -h33 * g43 * h44;
771 h32 = -h22 * g32 * h33;
772 h42 = -h22 * (g32 * h43 + g42 * h44);
773 h21 = -h11 * g21 * h22;
774 h31 = -h11 * (g21 * h32 + g31 * h33);
775 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
776 h10 = -h00 * g10 * h11;
777 h20 = -h00 * (g10 * h21 + g20 * h22);
778 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
779 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
784 m[
A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40;
785 m[
A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41;
786 m[
A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41;
787 m[
A02] = h20 * h22 + h30 * h32 + h40 * h42;
788 m[
A12] = h21 * h22 + h31 * h32 + h41 * h42;
789 m[
A22] = h22 * h22 + h32 * h32 + h42 * h42;
790 m[
A03] = h30 * h33 + h40 * h43;
791 m[
A13] = h31 * h33 + h41 * h43;
792 m[
A23] = h32 * h33 + h42 * h43;
793 m[
A33] = h33 * h33 + h43 * h43;
820 double h30, h31, h32;
821 double h40, h41, h42, h43;
822 double h50, h51, h52, h53, h54;
824 double h00, h11, h22, h33, h44, h55;
829 double g30, g31, g32;
830 double g40, g41, g42, g43;
831 double g50, g51, g52, g53, g54;
841 if (h00 <= 0)
return;
842 h00 = 1.0 / sqrt(h00);
852 h11 = m[
A11] - (g10 * g10);
853 if (h11 <= 0)
return;
854 h11 = 1.0 / sqrt(h11);
859 g21 = (m[
A21] - (g10 * g20)) * h11;
860 g31 = (m[
A31] - (g10 * g30)) * h11;
861 g41 = (m[
A41] - (g10 * g40)) * h11;
862 g51 = (m[
A51] - (g10 * g50)) * h11;
866 h22 = m[
A22] - (g20 * g20) - (g21 * g21);
867 if (h22 <= 0)
return;
868 h22 = 1.0 / sqrt(h22);
873 g32 = (m[
A32] - (g20 * g30) - (g21 * g31)) * h22;
874 g42 = (m[
A42] - (g20 * g40) - (g21 * g41)) * h22;
875 g52 = (m[
A52] - (g20 * g50) - (g21 * g51)) * h22;
879 h33 = m[
A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
880 if (h33 <= 0)
return;
881 h33 = 1.0 / sqrt(h33);
886 g43 = (m[
A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
887 g53 = (m[
A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33;
891 h44 = m[
A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
892 if (h44 <= 0)
return;
893 h44 = 1.0 / sqrt(h44);
897 g54 = (m[
A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;
901 h55 = m[
A55] - (g50*g50) - (g51*g51) - (g52*g52) - (g53*g53) - (g54*g54);
902 if (h55 <= 0)
return;
903 h55 = 1.0 / sqrt(h55);
910 h54 = -h44 * g54 * h55;
911 h43 = -h33 * g43 * h44;
912 h53 = -h33 * (g43 * h54 + g53 * h55);
913 h32 = -h22 * g32 * h33;
914 h42 = -h22 * (g32 * h43 + g42 * h44);
915 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * h55);
916 h21 = -h11 * g21 * h22;
917 h31 = -h11 * (g21 * h32 + g31 * h33);
918 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
919 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55);
920 h10 = -h00 * g10 * h11;
921 h20 = -h00 * (g10 * h21 + g20 * h22);
922 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
923 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
924 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55);
929 m[
A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40 + h50*h50;
930 m[
A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51;
931 m[
A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51;
932 m[
A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52;
933 m[
A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52;
934 m[
A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52;
935 m[
A03] = h30 * h33 + h40 * h43 + h50 * h53;
936 m[
A13] = h31 * h33 + h41 * h43 + h51 * h53;
937 m[
A23] = h32 * h33 + h42 * h43 + h52 * h53;
938 m[
A33] = h33 * h33 + h43 * h43 + h53 * h53;
939 m[
A04] = h40 * h44 + h50 * h54;
940 m[
A14] = h41 * h44 + h51 * h54;
941 m[
A24] = h42 * h44 + h52 * h54;
942 m[
A34] = h43 * h44 + h53 * h54;
943 m[
A44] = h44 * h44 + h54 * h54;
957 void HepSymMatrix::invert4 (
int & ifail) {
980 double Det3_012_012 = m[
A00]*Det2_12_12 - m[
A01]*Det2_12_02
982 double Det3_013_012 = m[
A00]*Det2_13_12 - m[
A01]*Det2_13_02
984 double Det3_013_013 = m[
A00]*Det2_13_13 - m[
A01]*Det2_13_03
986 double Det3_023_012 = m[
A00]*Det2_23_12 - m[
A01]*Det2_23_02
988 double Det3_023_013 = m[
A00]*Det2_23_13 - m[
A01]*Det2_23_03
990 double Det3_023_023 = m[
A00]*Det2_23_23 - m[
A02]*Det2_23_03
992 double Det3_123_012 = m[
A10]*Det2_23_12 - m[
A11]*Det2_23_02
994 double Det3_123_013 = m[
A10]*Det2_23_13 - m[
A11]*Det2_23_03
996 double Det3_123_023 = m[
A10]*Det2_23_23 - m[
A12]*Det2_23_03
998 double Det3_123_123 = m[
A11]*Det2_23_23 - m[
A12]*Det2_23_13
1003 double det = m[
A00]*Det3_123_123
1004 - m[
A01]*Det3_123_023
1005 + m[
A02]*Det3_123_013
1006 - m[
A03]*Det3_123_012;
1009 #ifdef SINGULAR_DIAGNOSTICS
1010 std::cerr <<
"Kramer's rule inversion of a singular 4x4 matrix: "
1017 double oneOverDet = 1.0/det;
1018 double mn1OverDet = - oneOverDet;
1020 m[
A00] = Det3_123_123 * oneOverDet;
1021 m[
A01] = Det3_123_023 * mn1OverDet;
1022 m[
A02] = Det3_123_013 * oneOverDet;
1023 m[
A03] = Det3_123_012 * mn1OverDet;
1026 m[
A11] = Det3_023_023 * oneOverDet;
1027 m[
A12] = Det3_023_013 * mn1OverDet;
1028 m[
A13] = Det3_023_012 * oneOverDet;
1030 m[
A22] = Det3_013_013 * oneOverDet;
1031 m[
A23] = Det3_013_012 * mn1OverDet;
1033 m[
A33] = Det3_012_012 * oneOverDet;
void invertHaywood5(int &ifail)
void invertCholesky6(int &ifail)
void invertHaywood4(int &ifail)
void invertCholesky5(int &ifail)
void invertHaywood6(int &ifail)