| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 // GEANT 4 class implementation file
29 // -------------------------------------------
30
31 #include "globals.hh"
32 #include <iostream>
33 #include <cmath>
34
35 #include "G4ErrorSymMatrix.hh"
36 #include "G4ErrorMatrix.hh"
37
38 // Simple operation for all elements
39
40 #define SIMPLE_UOP(OPER)
41 G4ErrorMatrixIter a = m.begin();
42 G4ErrorMatrixIter e = m.begin() + num_size()
43 for(; a < e; a++)
44 (*a) OPER t;
45
46 #define SIMPLE_BOP(OPER)
47 G4ErrorMatrixIter a = m.begin();
48 G4ErrorMatrixConstIter b = mat2.m.begin();
49 G4ErrorMatrixConstIter e = m.begin() + num_s
50 for(; a < e; a++, b++)
51 (*a) OPER(*b);
52
53 #define SIMPLE_TOP(OPER)
54 G4ErrorMatrixConstIter a = mat1.m.begin();
55 G4ErrorMatrixConstIter b = mat2.m.begin();
56 G4ErrorMatrixIter t = mret.m.begin();
57 G4ErrorMatrixConstIter e = mat1.m.begin() +
58 for(; a < e; a++, b++, t++)
59 (*t) = (*a) OPER(*b);
60
61 #define CHK_DIM_2(r1, r2, c1, c2, fun)
62 if(r1 != r2 || c1 != c2)
63 {
64 G4ErrorMatrix::error("Range error in Matri
65 }
66
67 #define CHK_DIM_1(c1, r2, fun)
68 if(c1 != r2)
69 {
70 G4ErrorMatrix::error("Range error in Matri
71 }
72
73 // Constructors. (Default constructors are inl
74
75 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p)
76 : m(p * (p + 1) / 2)
77 , nrow(p)
78 {
79 size = nrow * (nrow + 1) / 2;
80 m.assign(size, 0);
81 }
82
83 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p, G4
84 : m(p * (p + 1) / 2)
85 , nrow(p)
86 {
87 size = nrow * (nrow + 1) / 2;
88
89 m.assign(size, 0);
90 switch(init)
91 {
92 case 0:
93 break;
94
95 case 1:
96 {
97 G4ErrorMatrixIter a = m.begin();
98 for(G4int i = 1; i <= nrow; i++)
99 {
100 *a = 1.0;
101 a += (i + 1);
102 }
103 break;
104 }
105 default:
106 G4ErrorMatrix::error("G4ErrorSymMatrix:
107 }
108 }
109
110 //
111 // Destructor
112 //
113
114 G4ErrorSymMatrix::~G4ErrorSymMatrix() {}
115
116 G4ErrorSymMatrix::G4ErrorSymMatrix(const G4Err
117 : m(mat1.size)
118 , nrow(mat1.nrow)
119 , size(mat1.size)
120 {
121 m = mat1.m;
122 }
123
124 //
125 //
126 // Sub matrix
127 //
128 //
129
130 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int m
131 {
132 G4ErrorSymMatrix mret(max_row - min_row + 1)
133 if(max_row > num_row())
134 {
135 G4ErrorMatrix::error("G4ErrorSymMatrix::su
136 }
137 G4ErrorMatrixIter a = mret.m.begin();
138 G4ErrorMatrixConstIter b1 = m.begin() + (min
139 for(G4int irow = 1; irow <= mret.num_row();
140 {
141 G4ErrorMatrixConstIter b = b1;
142 for(G4int icol = 1; icol <= irow; icol++)
143 {
144 *(a++) = *(b++);
145 }
146 b1 += irow + min_row - 1;
147 }
148 return mret;
149 }
150
151 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int m
152 {
153 G4ErrorSymMatrix mret(max_row - min_row + 1)
154 if(max_row > num_row())
155 {
156 G4ErrorMatrix::error("G4ErrorSymMatrix::su
157 }
158 G4ErrorMatrixIter a = mret.m.begin();
159 G4ErrorMatrixIter b1 = m.begin() + (min_row
160 for(G4int irow = 1; irow <= mret.num_row();
161 {
162 G4ErrorMatrixIter b = b1;
163 for(G4int icol = 1; icol <= irow; icol++)
164 {
165 *(a++) = *(b++);
166 }
167 b1 += irow + min_row - 1;
168 }
169 return mret;
170 }
171
172 void G4ErrorSymMatrix::sub(G4int row, const G4
173 {
174 if(row < 1 || row + mat1.num_row() - 1 > num
175 {
176 G4ErrorMatrix::error("G4ErrorSymMatrix::su
177 }
178 G4ErrorMatrixConstIter a = mat1.m.begin();
179 G4ErrorMatrixIter b1 = m.begin() + (row
180 for(G4int irow = 1; irow <= mat1.num_row();
181 {
182 G4ErrorMatrixIter b = b1;
183 for(G4int icol = 1; icol <= irow; icol++)
184 {
185 *(b++) = *(a++);
186 }
187 b1 += irow + row - 1;
188 }
189 }
190
191 //
192 // Direct sum of two matricies
193 //
194
195 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix&
196 const G4ErrorSymMatrix&
197 {
198 G4ErrorSymMatrix mret(mat1.num_row() + mat2.
199 mret.sub(1, mat1);
200 mret.sub(mat1.num_row() + 1, mat2);
201 return mret;
202 }
203
204 /* -------------------------------------------
205 This section contains support routines for
206 The two argument functions +,-. They call t
207 -------------------------------------------
208
209 G4ErrorSymMatrix G4ErrorSymMatrix::operator-()
210 {
211 G4ErrorSymMatrix mat2(nrow);
212 G4ErrorMatrixConstIter a = m.begin();
213 G4ErrorMatrixIter b = mat2.m.begin();
214 G4ErrorMatrixConstIter e = m.begin() + num_s
215 for(; a < e; a++, b++)
216 {
217 (*b) = -(*a);
218 }
219 return mat2;
220 }
221
222 G4ErrorMatrix operator+(const G4ErrorMatrix& m
223 {
224 G4ErrorMatrix mret(mat1);
225 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma
226 mret += mat2;
227 return mret;
228 }
229
230 G4ErrorMatrix operator+(const G4ErrorSymMatrix
231 {
232 G4ErrorMatrix mret(mat2);
233 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma
234 mret += mat1;
235 return mret;
236 }
237
238 G4ErrorSymMatrix operator+(const G4ErrorSymMat
239 const G4ErrorSymMat
240 {
241 G4ErrorSymMatrix mret(mat1.nrow);
242 CHK_DIM_1(mat1.nrow, mat2.nrow, +);
243 SIMPLE_TOP(+)
244 return mret;
245 }
246
247 //
248 // operator -
249 //
250
251 G4ErrorMatrix operator-(const G4ErrorMatrix& m
252 {
253 G4ErrorMatrix mret(mat1);
254 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma
255 mret -= mat2;
256 return mret;
257 }
258
259 G4ErrorMatrix operator-(const G4ErrorSymMatrix
260 {
261 G4ErrorMatrix mret(mat1);
262 CHK_DIM_2(mat1.num_row(), mat2.num_row(), ma
263 mret -= mat2;
264 return mret;
265 }
266
267 G4ErrorSymMatrix operator-(const G4ErrorSymMat
268 const G4ErrorSymMat
269 {
270 G4ErrorSymMatrix mret(mat1.num_row());
271 CHK_DIM_1(mat1.num_row(), mat2.num_row(), -)
272 SIMPLE_TOP(-)
273 return mret;
274 }
275
276 /* -------------------------------------------
277 This section contains support routines for
278 The two argument functions *,/. They call c
279 -------------------------------------------
280
281 G4ErrorSymMatrix operator/(const G4ErrorSymMat
282 {
283 G4ErrorSymMatrix mret(mat1);
284 mret /= t;
285 return mret;
286 }
287
288 G4ErrorSymMatrix operator*(const G4ErrorSymMat
289 {
290 G4ErrorSymMatrix mret(mat1);
291 mret *= t;
292 return mret;
293 }
294
295 G4ErrorSymMatrix operator*(G4double t, const G
296 {
297 G4ErrorSymMatrix mret(mat1);
298 mret *= t;
299 return mret;
300 }
301
302 G4ErrorMatrix operator*(const G4ErrorMatrix& m
303 {
304 G4ErrorMatrix mret(mat1.num_row(), mat2.num_
305 CHK_DIM_1(mat1.num_col(), mat2.num_row(), *)
306 G4ErrorMatrixConstIter mit1, mit2, sp, snp;
307 G4double temp;
308 G4ErrorMatrixIter mir = mret.m.begin();
309 for(mit1 = mat1.m.begin();
310 mit1 < mat1.m.begin() + mat1.num_row() *
311 {
312 snp = mat2.m.begin();
313 for(int step = 1; step <= mat2.num_row();
314 {
315 mit2 = mit1;
316 sp = snp;
317 snp += step;
318 temp = 0;
319 while(sp < snp) // Loop checking, 06.08
320 {
321 temp += *(sp++) * (*(mit2++));
322 }
323 if(step < mat2.num_row())
324 { // only if we aren't on the last row
325 sp += step - 1;
326 for(int stept = step + 1; stept <= mat
327 {
328 temp += *sp * (*(mit2++));
329 if(stept < mat2.num_row())
330 sp += stept;
331 }
332 } // if(step
333 *(mir++) = temp;
334 } // for(step
335 } // for(mit1
336 return mret;
337 }
338
339 G4ErrorMatrix operator*(const G4ErrorSymMatrix
340 {
341 G4ErrorMatrix mret(mat1.num_row(), mat2.num_
342 CHK_DIM_1(mat1.num_col(), mat2.num_row(), *)
343 G4int step, stept;
344 G4ErrorMatrixConstIter mit1, mit2, sp, snp;
345 G4double temp;
346 G4ErrorMatrixIter mir = mret.m.begin();
347 for(step = 1, snp = mat1.m.begin(); step <=
348 {
349 for(mit1 = mat2.m.begin(); mit1 < mat2.m.b
350 {
351 mit2 = mit1;
352 sp = snp;
353 temp = 0;
354 while(sp < snp + step) // Loop checking
355 {
356 temp += *mit2 * (*(sp++));
357 mit2 += mat2.num_col();
358 }
359 sp += step - 1;
360 for(stept = step + 1; stept <= mat1.num_
361 {
362 temp += *mit2 * (*sp);
363 mit2 += mat2.num_col();
364 sp += stept;
365 }
366 *(mir++) = temp;
367 }
368 }
369 return mret;
370 }
371
372 G4ErrorMatrix operator*(const G4ErrorSymMatrix
373 const G4ErrorSymMatrix
374 {
375 G4ErrorMatrix mret(mat1.num_row(), mat1.num_
376 CHK_DIM_1(mat1.num_col(), mat2.num_row(), *)
377 G4int step1, stept1, step2, stept2;
378 G4ErrorMatrixConstIter snp1, sp1, snp2, sp2;
379 G4double temp;
380 G4ErrorMatrixIter mr = mret.m.begin();
381 for(step1 = 1, snp1 = mat1.m.begin(); step1
382 snp1 += step1++)
383 {
384 for(step2 = 1, snp2 = mat2.m.begin(); step
385 {
386 sp1 = snp1;
387 sp2 = snp2;
388 snp2 += step2;
389 temp = 0;
390 if(step1 < step2)
391 {
392 while(sp1 < snp1 + step1) // Loop che
393 {
394 temp += (*(sp1++)) * (*(sp2++));
395 }
396 sp1 += step1 - 1;
397 for(stept1 = step1 + 1; stept1 != step
398 {
399 temp += (*sp1) * (*(sp2++));
400 }
401 sp2 += step2 - 1;
402 for(stept2 = ++step2; stept2 <= mat2.n
403 sp1 += stept1++, sp2 += stept2++)
404 {
405 temp += (*sp1) * (*sp2);
406 }
407 }
408 else
409 {
410 while(sp2 < snp2) // Loop checking, 0
411 {
412 temp += (*(sp1++)) * (*(sp2++));
413 }
414 sp2 += step2 - 1;
415 for(stept2 = ++step2; stept2 != step1
416 {
417 temp += (*(sp1++)) * (*sp2);
418 }
419 sp1 += step1 - 1;
420 for(stept1 = step1 + 1; stept1 <= mat1
421 sp1 += stept1++, sp2 += stept2++)
422 {
423 temp += (*sp1) * (*sp2);
424 }
425 }
426 *(mr++) = temp;
427 }
428 }
429 return mret;
430 }
431
432 /* -------------------------------------------
433 This section contains the assignment and in
434 -------------------------------------------
435
436 G4ErrorMatrix& G4ErrorMatrix::operator+=(const
437 {
438 CHK_DIM_2(num_row(), mat2.num_row(), num_col
439 G4int n = num_col();
440 G4ErrorMatrixConstIter sjk = mat2.m.begin();
441 G4ErrorMatrixIter m1j = m.begin();
442 G4ErrorMatrixIter mj = m.begin();
443 // j >= k
444 for(G4int j = 1; j <= num_row(); j++)
445 {
446 G4ErrorMatrixIter mjk = mj;
447 G4ErrorMatrixIter mkj = m1j;
448 for(G4int k = 1; k <= j; k++)
449 {
450 *(mjk++) += *sjk;
451 if(j != k)
452 *mkj += *sjk;
453 sjk++;
454 mkj += n;
455 }
456 mj += n;
457 m1j++;
458 }
459 return (*this);
460 }
461
462 G4ErrorSymMatrix& G4ErrorSymMatrix::operator+=
463 {
464 CHK_DIM_2(num_row(), mat2.num_row(), num_col
465 SIMPLE_BOP(+=)
466 return (*this);
467 }
468
469 G4ErrorMatrix& G4ErrorMatrix::operator-=(const
470 {
471 CHK_DIM_2(num_row(), mat2.num_row(), num_col
472 G4int n = num_col();
473 G4ErrorMatrixConstIter sjk = mat2.m.begin();
474 G4ErrorMatrixIter m1j = m.begin();
475 G4ErrorMatrixIter mj = m.begin();
476 // j >= k
477 for(G4int j = 1; j <= num_row(); j++)
478 {
479 G4ErrorMatrixIter mjk = mj;
480 G4ErrorMatrixIter mkj = m1j;
481 for(G4int k = 1; k <= j; k++)
482 {
483 *(mjk++) -= *sjk;
484 if(j != k)
485 *mkj -= *sjk;
486 sjk++;
487 mkj += n;
488 }
489 mj += n;
490 m1j++;
491 }
492 return (*this);
493 }
494
495 G4ErrorSymMatrix& G4ErrorSymMatrix::operator-=
496 {
497 CHK_DIM_2(num_row(), mat2.num_row(), num_col
498 SIMPLE_BOP(-=)
499 return (*this);
500 }
501
502 G4ErrorSymMatrix& G4ErrorSymMatrix::operator/=
503 {
504 SIMPLE_UOP(/=)
505 return (*this);
506 }
507
508 G4ErrorSymMatrix& G4ErrorSymMatrix::operator*=
509 {
510 SIMPLE_UOP(*=)
511 return (*this);
512 }
513
514 G4ErrorMatrix& G4ErrorMatrix::operator=(const
515 {
516 if(mat1.nrow * mat1.nrow != size)
517 {
518 size = mat1.nrow * mat1.nrow;
519 m.resize(size);
520 }
521 nrow = mat1.nrow;
522 ncol = mat1.nrow;
523 G4int n = ncol;
524 G4ErrorMatrixConstIter sjk = mat1.m.begin();
525 G4ErrorMatrixIter m1j = m.begin();
526 G4ErrorMatrixIter mj = m.begin();
527 // j >= k
528 for(G4int j = 1; j <= num_row(); j++)
529 {
530 G4ErrorMatrixIter mjk = mj;
531 G4ErrorMatrixIter mkj = m1j;
532 for(G4int k = 1; k <= j; k++)
533 {
534 *(mjk++) = *sjk;
535 if(j != k)
536 *mkj = *sjk;
537 sjk++;
538 mkj += n;
539 }
540 mj += n;
541 m1j++;
542 }
543 return (*this);
544 }
545
546 G4ErrorSymMatrix& G4ErrorSymMatrix::operator=(
547 {
548 if(&mat1 == this)
549 {
550 return *this;
551 }
552 if(mat1.nrow != nrow)
553 {
554 nrow = mat1.nrow;
555 size = mat1.size;
556 m.resize(size);
557 }
558 m = mat1.m;
559 return (*this);
560 }
561
562 // Print the Matrix.
563
564 std::ostream& operator<<(std::ostream& os, con
565 {
566 os << G4endl;
567
568 // Fixed format needs 3 extra characters for
569 // while scientific needs 7
570
571 std::size_t width;
572 if(os.flags() & std::ios::fixed)
573 {
574 width = os.precision() + 3;
575 }
576 else
577 {
578 width = os.precision() + 7;
579 }
580 for(G4int irow = 1; irow <= q.num_row(); ++i
581 {
582 for(G4int icol = 1; icol <= q.num_col(); +
583 {
584 os.width(width);
585 os << q(irow, icol) << " ";
586 }
587 os << G4endl;
588 }
589 return os;
590 }
591
592 G4ErrorSymMatrix G4ErrorSymMatrix::apply(G4dou
593
594 {
595 G4ErrorSymMatrix mret(num_row());
596 G4ErrorMatrixConstIter a = m.cbegin();
597 G4ErrorMatrixIter b = mret.m.begin();
598 for(G4int ir = 1; ir <= num_row(); ++ir)
599 {
600 for(G4int ic = 1; ic <= ir; ++ic)
601 {
602 *(b++) = (*f)(*(a++), ir, ic);
603 }
604 }
605 return mret;
606 }
607
608 void G4ErrorSymMatrix::assign(const G4ErrorMat
609 {
610 if(mat1.nrow != nrow)
611 {
612 nrow = mat1.nrow;
613 size = nrow * (nrow + 1) / 2;
614 m.resize(size);
615 }
616 G4ErrorMatrixConstIter a = mat1.m.begin();
617 G4ErrorMatrixIter b = m.begin();
618 for(G4int r = 1; r <= nrow; r++)
619 {
620 G4ErrorMatrixConstIter d = a;
621 for(G4int c = 1; c <= r; c++)
622 {
623 *(b++) = *(d++);
624 }
625 a += nrow;
626 }
627 }
628
629 G4ErrorSymMatrix G4ErrorSymMatrix::similarity(
630 {
631 G4ErrorSymMatrix mret(mat1.num_row());
632 G4ErrorMatrix temp = mat1 * (*this);
633
634 // If mat1*(*this) has correct dimensions, t
635 // multiplication. So there is no need to ch
636
637 G4int n = mat1.num_col();
638 G4ErrorMatrixIter mr = mret.m.begin();
639 G4ErrorMatrixIter tempr1 = temp.m.begin();
640 for(G4int r = 1; r <= mret.num_row(); r++)
641 {
642 G4ErrorMatrixConstIter m1c1 = mat1.m.begin
643 for(G4int c = 1; c <= r; c++)
644 {
645 G4double tmp = 0.0;
646 G4ErrorMatrixIter tempri = tempr1;
647 G4ErrorMatrixConstIter m1ci = m1c1;
648 for(G4int i = 1; i <= mat1.num_col(); i+
649 {
650 tmp += (*(tempri++)) * (*(m1ci++));
651 }
652 *(mr++) = tmp;
653 m1c1 += n;
654 }
655 tempr1 += n;
656 }
657 return mret;
658 }
659
660 G4ErrorSymMatrix G4ErrorSymMatrix::similarity(
661 const G4ErrorSymMatrix& mat1) const
662 {
663 G4ErrorSymMatrix mret(mat1.num_row());
664 G4ErrorMatrix temp = mat1 * (*this);
665 G4int n = mat1.num_col();
666 G4ErrorMatrixIter mr = mret.m.begin();
667 G4ErrorMatrixIter tempr1 = temp.m.begin();
668 for(G4int r = 1; r <= mret.num_row(); r++)
669 {
670 G4ErrorMatrixConstIter m1c1 = mat1.m.begin
671 G4int c;
672 for(c = 1; c <= r; c++)
673 {
674 G4double tmp = 0.0;
675 G4ErrorMatrixIter tempri = tempr1;
676 G4ErrorMatrixConstIter m1ci = m1c1;
677 G4int i;
678 for(i = 1; i < c; i++)
679 {
680 tmp += (*(tempri++)) * (*(m1ci++));
681 }
682 for(i = c; i <= mat1.num_col(); i++)
683 {
684 tmp += (*(tempri++)) * (*(m1ci));
685 m1ci += i;
686 }
687 *(mr++) = tmp;
688 m1c1 += c;
689 }
690 tempr1 += n;
691 }
692 return mret;
693 }
694
695 G4ErrorSymMatrix G4ErrorSymMatrix::similarityT
696 {
697 G4ErrorSymMatrix mret(mat1.num_col());
698 G4ErrorMatrix temp = (*this) * mat1;
699 G4int n = mat1.num_col();
700 G4ErrorMatrixIter mrc = mret.m.begin();
701 G4ErrorMatrixIter temp1r = temp.m.begin();
702 for(G4int r = 1; r <= mret.num_row(); r++)
703 {
704 G4ErrorMatrixConstIter m11c = mat1.m.begin
705 for(G4int c = 1; c <= r; c++)
706 {
707 G4double tmp = 0.0;
708 G4ErrorMatrixIter tempir = temp1r;
709 G4ErrorMatrixConstIter m1ic = m11c;
710 for(G4int i = 1; i <= mat1.num_row(); i+
711 {
712 tmp += (*(tempir)) * (*(m1ic));
713 tempir += n;
714 m1ic += n;
715 }
716 *(mrc++) = tmp;
717 m11c++;
718 }
719 temp1r++;
720 }
721 return mret;
722 }
723
724 void G4ErrorSymMatrix::invert(G4int& ifail)
725 {
726 ifail = 0;
727
728 switch(nrow)
729 {
730 case 3:
731 {
732 G4double det, temp;
733 G4double t1, t2, t3;
734 G4double c11, c12, c13, c22, c23, c33;
735 c11 = (*(m.begin() + 2)) * (*(m.begin()
736 (*(m.begin() + 4)) * (*(m.begin()
737 c12 = (*(m.begin() + 4)) * (*(m.begin()
738 (*(m.begin() + 1)) * (*(m.begin()
739 c13 = (*(m.begin() + 1)) * (*(m.begin()
740 (*(m.begin() + 2)) * (*(m.begin()
741 c22 = (*(m.begin() + 5)) * (*m.begin())
742 (*(m.begin() + 3)) * (*(m.begin()
743 c23 = (*(m.begin() + 3)) * (*(m.begin()
744 (*(m.begin() + 4)) * (*m.begin());
745 c33 = (*m.begin()) * (*(m.begin() + 2))
746 (*(m.begin() + 1)) * (*(m.begin()
747 t1 = std::fabs(*m.begin());
748 t2 = std::fabs(*(m.begin() + 1));
749 t3 = std::fabs(*(m.begin() + 3));
750 if(t1 >= t2)
751 {
752 if(t3 >= t1)
753 {
754 temp = *(m.begin() + 3);
755 det = c23 * c12 - c22 * c13;
756 }
757 else
758 {
759 temp = *m.begin();
760 det = c22 * c33 - c23 * c23;
761 }
762 }
763 else if(t3 >= t2)
764 {
765 temp = *(m.begin() + 3);
766 det = c23 * c12 - c22 * c13;
767 }
768 else
769 {
770 temp = *(m.begin() + 1);
771 det = c13 * c23 - c12 * c33;
772 }
773 if(det == 0)
774 {
775 ifail = 1;
776 return;
777 }
778 {
779 G4double ss = temp / det;
780 G4ErrorMatrixIter mq = m.begin();
781 *(mq++) = ss * c11;
782 *(mq++) = ss * c12;
783 *(mq++) = ss * c22;
784 *(mq++) = ss * c13;
785 *(mq++) = ss * c23;
786 *(mq) = ss * c33;
787 }
788 }
789 break;
790 case 2:
791 {
792 G4double det, temp, ss;
793 det = (*m.begin()) * (*(m.begin() + 2))
794 (*(m.begin() + 1)) * (*(m.begin()
795 if(det == 0)
796 {
797 ifail = 1;
798 return;
799 }
800 ss = 1.0 / det;
801 *(m.begin() + 1) *= -ss;
802 temp = ss * (*(m.begin() + 2
803 *(m.begin() + 2) = ss * (*m.begin());
804 *m.begin() = temp;
805 break;
806 }
807 case 1:
808 {
809 if((*m.begin()) == 0)
810 {
811 ifail = 1;
812 return;
813 }
814 *m.begin() = 1.0 / (*m.begin());
815 break;
816 }
817 case 5:
818 {
819 invert5(ifail);
820 return;
821 }
822 case 6:
823 {
824 invert6(ifail);
825 return;
826 }
827 case 4:
828 {
829 invert4(ifail);
830 return;
831 }
832 default:
833 {
834 invertBunchKaufman(ifail);
835 return;
836 }
837 }
838 return; // inversion successful
839 }
840
841 G4double G4ErrorSymMatrix::determinant() const
842 {
843 static const G4int max_array = 20;
844
845 // ir must point to an array which is ***1 l
846
847 static std::vector<G4int> ir_vec(max_array +
848 if(ir_vec.size() <= static_cast<unsigned int
849 {
850 ir_vec.resize(nrow + 1);
851 }
852 G4int* ir = &ir_vec[0];
853
854 G4double det;
855 G4ErrorMatrix mt(*this);
856 G4int i = mt.dfact_matrix(det, ir);
857 if(i == 0)
858 {
859 return det;
860 }
861 return 0.0;
862 }
863
864 G4double G4ErrorSymMatrix::trace() const
865 {
866 G4double t = 0.0;
867 for(G4int i = 0; i < nrow; i++)
868 {
869 t += *(m.begin() + (i + 3) * i / 2);
870 }
871 return t;
872 }
873
874 void G4ErrorSymMatrix::invertBunchKaufman(G4in
875 {
876 // Bunch-Kaufman diagonal pivoting method
877 // It is decribed in J.R. Bunch, L. Kaufman
878 // "Some Stable Methods for Calculating Iner
879 // Linear Systems", Math. Comp. 31, p. 162-1
880 // Charles F. van Loan, "Matrix Computations
881 // has a bug.) and implemented in "lapack"
882 // Mario Stanke, 09/97
883
884 G4int i, j, k, ss;
885 G4int pivrow;
886
887 // Establish the two working-space arrays ne
888 // used as pointers to arrays of doubles and
889 // of length nrow. We do not want to reallo
890 // unless the size needs to grow. We do not
891 // by having a new without a delete that is
892
893 static const G4int max_array
894 static G4ThreadLocal std::vector<G4double>*
895 if(!xvec)
896 xvec = new std::vector<G4double>(max_array
897 static G4ThreadLocal std::vector<G4int>* piv
898 if(!pivv)
899 pivv = new std::vector<G4int>(max_array);
900 typedef std::vector<G4int>::iterator pivIter
901 if(xvec->size() < static_cast<unsigned int>(
902 xvec->resize(nrow);
903 if(pivv->size() < static_cast<unsigned int>(
904 pivv->resize(nrow);
905 // Note - resize should do nothing if the s
906 // but on VC++ there are indications
907 // Note - the data elements in a vector are
908 // so x[i] and piv[i] are optimally f
909 G4ErrorMatrixIter x = xvec->begin();
910 // x[i] is used as helper storage, needs to
911 pivIter piv = pivv->begin();
912 // piv[i] is used to store details of exchan
913
914 G4double temp1, temp2;
915 G4ErrorMatrixIter ip, mjj, iq;
916 G4double lambda, sigma;
917 const G4double alpha = .6404; // = (1+sqr
918 const G4double epsilon = 32 * DBL_EPSILON;
919 // whenever a sum of two doubles is below or
920 // it is set to zero.
921 // this constant could be set to zero but th
922 // doesn't neccessarily detect that a matrix
923
924 for(i = 0; i < nrow; i++)
925 {
926 piv[i] = i + 1;
927 }
928
929 ifail = 0;
930
931 // compute the factorization P*A*P^T = L * D
932 // L is unit lower triangular, D is direct s
933 // L and D^-1 are stored in A = *this, P is
934
935 for(j = 1; j < nrow; j += ss) // main loop
936 {
937 mjj = m.begin() + j * (j - 1) / 2 + j -
938 lambda = 0; // compute lambda = max of A(
939 pivrow = j + 1;
940 ip = m.begin() + (j + 1) * j / 2 + j -
941 for(i = j + 1; i <= nrow; ip += i++)
942 {
943 if(std::fabs(*ip) > lambda)
944 {
945 lambda = std::fabs(*ip);
946 pivrow = i;
947 }
948 }
949 if(lambda == 0)
950 {
951 if(*mjj == 0)
952 {
953 ifail = 1;
954 return;
955 }
956 ss = 1;
957 *mjj = 1. / *mjj;
958 }
959 else
960 {
961 if(std::fabs(*mjj) >= lambda * alpha)
962 {
963 ss = 1;
964 pivrow = j;
965 }
966 else
967 {
968 sigma = 0; // compute sigma = max A(p
969 ip = m.begin() + pivrow * (pivrow -
970 for(k = j; k < pivrow; k++)
971 {
972 if(std::fabs(*ip) > sigma)
973 sigma = std::fabs(*ip);
974 ip++;
975 }
976 if(sigma * std::fabs(*mjj) >= alpha *
977 {
978 ss = 1;
979 pivrow = j;
980 }
981 else if(std::fabs(*(m.begin() + pivrow
982 1)) >= alpha * sig
983 {
984 ss = 1;
985 }
986 else
987 {
988 ss = 2;
989 }
990 }
991 if(pivrow == j) // no permutation necce
992 {
993 piv[j - 1] = pivrow;
994 if(*mjj == 0)
995 {
996 ifail = 1;
997 return;
998 }
999 temp2 = *mjj = 1. / *mjj; // invert D
1000
1001 // update A(j+1:n, j+1,n)
1002 for(i = j + 1; i <= nrow; i++)
1003 {
1004 temp1 = *(m.begin() + i * (i - 1) /
1005 ip = m.begin() + i * (i - 1) / 2
1006 for(k = j + 1; k <= i; k++)
1007 {
1008 *ip -= temp1 * *(m.begin() + k *
1009 if(std::fabs(*ip) <= epsilon)
1010 {
1011 *ip = 0;
1012 }
1013 ip++;
1014 }
1015 }
1016 // update L
1017 ip = m.begin() + (j + 1) * j / 2 + j
1018 for(i = j + 1; i <= nrow; ip += i++)
1019 {
1020 *ip *= temp2;
1021 }
1022 }
1023 else if(ss == 1) // 1x1 pivot
1024 {
1025 piv[j - 1] = pivrow;
1026
1027 // interchange rows and columns j and
1028 // submatrix (j:n,j:n)
1029 ip = m.begin() + pivrow * (pivrow - 1
1030 for(i = j + 1; i < pivrow; i++, ip++)
1031 {
1032 temp1 = *(m.begin() + i * (i - 1) /
1033 *(m.begin() + i * (i - 1) / 2 + j -
1034 *ip
1035 }
1036 temp1 = *mjj;
1037 *mjj = *(m.begin() + pivrow * (pivro
1038 *(m.begin() + pivrow * (pivrow - 1) /
1039 ip = m.begin() + (pivrow + 1) * pivro
1040 iq = ip + pivrow - j;
1041 for(i = pivrow + 1; i <= nrow; ip +=
1042 {
1043 temp1 = *iq;
1044 *iq = *ip;
1045 *ip = temp1;
1046 }
1047
1048 if(*mjj == 0)
1049 {
1050 ifail = 1;
1051 return;
1052 }
1053 temp2 = *mjj = 1. / *mjj; // invert
1054
1055 // update A(j+1:n, j+1:n)
1056 for(i = j + 1; i <= nrow; i++)
1057 {
1058 temp1 = *(m.begin() + i * (i - 1) /
1059 ip = m.begin() + i * (i - 1) / 2
1060 for(k = j + 1; k <= i; k++)
1061 {
1062 *ip -= temp1 * *(m.begin() + k *
1063 if(std::fabs(*ip) <= epsilon)
1064 {
1065 *ip = 0;
1066 }
1067 ip++;
1068 }
1069 }
1070 // update L
1071 ip = m.begin() + (j + 1) * j / 2 + j
1072 for(i = j + 1; i <= nrow; ip += i++)
1073 {
1074 *ip *= temp2;
1075 }
1076 }
1077 else // ss=2, ie use a 2x2 pivot
1078 {
1079 piv[j - 1] = -pivrow;
1080 piv[j] = 0; // that means this i
1081
1082 if(j + 1 != pivrow)
1083 {
1084 // interchange rows and columns j+1
1085 // submatrix (j:n,j:n)
1086 ip = m.begin() + pivrow * (pivrow -
1087 for(i = j + 2; i < pivrow; i++, ip+
1088 {
1089 temp1 = *(m.begin() + i * (i - 1)
1090 *(m.begin() + i * (i - 1) / 2 + j
1091 *ip
1092 }
1093 temp1 = *(mjj + j + 1);
1094 *(mjj + j + 1) =
1095 *(m.begin() + pivrow * (pivrow -
1096 *(m.begin() + pivrow * (pivrow - 1)
1097 temp1
1098 *(mjj + j) = *(m.begin() + pivrow *
1099 *(m.begin() + pivrow * (pivrow - 1)
1100 ip = m.begin() + (pivrow + 1) * piv
1101 iq = ip + pivrow - (j + 1);
1102 for(i = pivrow + 1; i <= nrow; ip +
1103 {
1104 temp1 = *iq;
1105 *iq = *ip;
1106 *ip = temp1;
1107 }
1108 }
1109 // invert D(j:j+1,j:j+1)
1110 temp2 = *mjj * *(mjj + j + 1) - *(mjj
1111 if(temp2 == 0)
1112 {
1113 G4Exception("G4ErrorSymMatrix::bunc
1114 "GEANT4e-Notification",
1115 "Error in pivot choice!
1116 }
1117 temp2 = 1. / temp2;
1118
1119 // this quotient is guaranteed to exi
1120 // of the pivot
1121
1122 temp1 = *mjj;
1123 *mjj = *(mjj + j + 1) * tem
1124 *(mjj + j + 1) = temp1 * temp2;
1125 *(mjj + j) = -*(mjj + j) * temp2;
1126
1127 if(j < nrow - 1) // otherwise do not
1128 {
1129 // update A(j+2:n, j+2:n)
1130 for(i = j + 2; i <= nrow; i++)
1131 {
1132 ip = m.begin() + i * (i - 1) /
1133 temp1 = *ip * *mjj + *(ip + 1) *
1134 if(std::fabs(temp1) <= epsilon)
1135 {
1136 temp1 = 0;
1137 }
1138 temp2 = *ip * *(mjj + j) + *(ip +
1139 if(std::fabs(temp2) <= epsilon)
1140 {
1141 temp2 = 0;
1142 }
1143 for(k = j + 2; k <= i; k++)
1144 {
1145 ip = m.begin() + i * (i - 1) /
1146 iq = m.begin() + k * (k - 1) /
1147 *ip -= temp1 * *iq + temp2 * *(
1148 if(std::fabs(*ip) <= epsilon)
1149 {
1150 *ip = 0;
1151 }
1152 }
1153 }
1154 // update L
1155 for(i = j + 2; i <= nrow; i++)
1156 {
1157 ip = m.begin() + i * (i - 1) /
1158 temp1 = *ip * *mjj + *(ip + 1) *
1159 if(std::fabs(temp1) <= epsilon)
1160 {
1161 temp1 = 0;
1162 }
1163 *(ip + 1) = *ip * *(mjj + j) + *(
1164 if(std::fabs(*(ip + 1)) <= epsilo
1165 {
1166 *(ip + 1) = 0;
1167 }
1168 *ip = temp1;
1169 }
1170 }
1171 }
1172 }
1173 } // end of main loop over columns
1174
1175 if(j == nrow) // the the last pivot is 1x1
1176 {
1177 mjj = m.begin() + j * (j - 1) / 2 + j - 1
1178 if(*mjj == 0)
1179 {
1180 ifail = 1;
1181 return;
1182 }
1183 else
1184 {
1185 *mjj = 1. / *mjj;
1186 }
1187 } // end of last pivot code
1188
1189 // computing the inverse from the factoriza
1190
1191 for(j = nrow; j >= 1; j -= ss) // loop ove
1192 {
1193 mjj = m.begin() + j * (j - 1) / 2 + j - 1
1194 if(piv[j - 1] > 0) // 1x1 pivot, compute
1195 {
1196 ss = 1;
1197 if(j < nrow)
1198 {
1199 ip = m.begin() + (j + 1) * j / 2 + j
1200 for(i = 0; i < nrow - j; ip += 1 + j
1201 {
1202 x[i] = *ip;
1203 }
1204 for(i = j + 1; i <= nrow; i++)
1205 {
1206 temp2 = 0;
1207 ip = m.begin() + i * (i - 1) / 2
1208 for(k = 0; k <= i - j - 1; k++)
1209 {
1210 temp2 += *ip++ * x[k];
1211 }
1212 for(ip += i - 1; k < nrow - j; ip +
1213 {
1214 temp2 += *ip * x[k];
1215 }
1216 *(m.begin() + i * (i - 1) / 2 + j -
1217 }
1218 temp2 = 0;
1219 ip = m.begin() + (j + 1) * j / 2 +
1220 for(k = 0; k < nrow - j; ip += 1 + j
1221 {
1222 temp2 += x[k] * *ip;
1223 }
1224 *mjj -= temp2;
1225 }
1226 }
1227 else // 2x2 pivot, compute columns j and
1228 {
1229 if(piv[j - 1] != 0)
1230 {
1231 std::ostringstream message;
1232 message << "Error in pivot: " << piv[
1233 G4Exception("G4ErrorSymMatrix::invert
1234 "GEANT4e-Notification", J
1235 }
1236 ss = 2;
1237 if(j < nrow)
1238 {
1239 ip = m.begin() + (j + 1) * j / 2 + j
1240 for(i = 0; i < nrow - j; ip += 1 + j
1241 {
1242 x[i] = *ip;
1243 }
1244 for(i = j + 1; i <= nrow; i++)
1245 {
1246 temp2 = 0;
1247 ip = m.begin() + i * (i - 1) / 2
1248 for(k = 0; k <= i - j - 1; k++)
1249 {
1250 temp2 += *ip++ * x[k];
1251 }
1252 for(ip += i - 1; k < nrow - j; ip +
1253 {
1254 temp2 += *ip * x[k];
1255 }
1256 *(m.begin() + i * (i - 1) / 2 + j -
1257 }
1258 temp2 = 0;
1259 ip = m.begin() + (j + 1) * j / 2 +
1260 for(k = 0; k < nrow - j; ip += 1 + j
1261 {
1262 temp2 += x[k] * *ip;
1263 }
1264 *mjj -= temp2;
1265 temp2 = 0;
1266 ip = m.begin() + (j + 1) * j / 2 +
1267 for(i = j + 1; i <= nrow; ip += i++)
1268 {
1269 temp2 += *ip * *(ip + 1);
1270 }
1271 *(mjj - 1) -= temp2;
1272 ip = m.begin() + (j + 1) * j / 2 + j
1273 for(i = 0; i < nrow - j; ip += 1 + j
1274 {
1275 x[i] = *ip;
1276 }
1277 for(i = j + 1; i <= nrow; i++)
1278 {
1279 temp2 = 0;
1280 ip = m.begin() + i * (i - 1) / 2
1281 for(k = 0; k <= i - j - 1; k++)
1282 {
1283 temp2 += *ip++ * x[k];
1284 }
1285 for(ip += i - 1; k < nrow - j; ip +
1286 {
1287 temp2 += *ip * x[k];
1288 }
1289 *(m.begin() + i * (i - 1) / 2 + j -
1290 }
1291 temp2 = 0;
1292 ip = m.begin() + (j + 1) * j / 2 +
1293 for(k = 0; k < nrow - j; ip += 1 + j
1294 {
1295 temp2 += x[k] * *ip;
1296 }
1297 *(mjj - j) -= temp2;
1298 }
1299 }
1300
1301 // interchange rows and columns j and piv
1302 // or rows and columns j and -piv[j-2]
1303
1304 pivrow = (piv[j - 1] == 0) ? -piv[j - 2]
1305 ip = m.begin() + pivrow * (pivrow - 1
1306 for(i = j + 1; i < pivrow; i++, ip++)
1307 {
1308 temp1 = *(m.begin() + i * (i - 1) / 2 +
1309 *(m.begin() + i * (i - 1) / 2 + j - 1)
1310 *ip
1311 }
1312 temp1 = *mjj;
1313 *mjj = *(m.begin() + pivrow * (pivrow -
1314 *(m.begin() + pivrow * (pivrow - 1) / 2 +
1315 if(ss == 2)
1316 {
1317 temp1 = *(mjj - 1);
1318 *(mjj - 1) = *(m.begin() + pivrow * (pi
1319 *(m.begin() + pivrow * (pivrow - 1) / 2
1320 }
1321
1322 ip = m.begin() + (pivrow + 1) * pivrow /
1323 iq = ip + pivrow - j;
1324 for(i = pivrow + 1; i <= nrow; ip += i, i
1325 {
1326 temp1 = *iq;
1327 *iq = *ip;
1328 *ip = temp1;
1329 }
1330 } // end of loop over columns (in computin
1331
1332 return; // inversion successful
1333 }
1334
1335 G4ThreadLocal G4double G4ErrorSymMatrix::posD
1336 G4ThreadLocal G4double G4ErrorSymMatrix::posD
1337 G4ThreadLocal G4double G4ErrorSymMatrix::adju
1338 G4ThreadLocal G4double G4ErrorSymMatrix::adju
1339 const G4double G4ErrorSymMatrix::CHOLESKY_THR
1340 const G4double G4ErrorSymMatrix::CHOLESKY_THR
1341 const G4double G4ErrorSymMatrix::CHOLESKY_CRE
1342 const G4double G4ErrorSymMatrix::CHOLESKY_CRE
1343
1344 // Aij are indices for a 6x6 symmetric matrix
1345 // The indices for 5x5 or 4x4 symmetric m
1346 // ignoring all combinations with an inde
1347
1348 #define A00 0
1349 #define A01 1
1350 #define A02 3
1351 #define A03 6
1352 #define A04 10
1353 #define A05 15
1354
1355 #define A10 1
1356 #define A11 2
1357 #define A12 4
1358 #define A13 7
1359 #define A14 11
1360 #define A15 16
1361
1362 #define A20 3
1363 #define A21 4
1364 #define A22 5
1365 #define A23 8
1366 #define A24 12
1367 #define A25 17
1368
1369 #define A30 6
1370 #define A31 7
1371 #define A32 8
1372 #define A33 9
1373 #define A34 13
1374 #define A35 18
1375
1376 #define A40 10
1377 #define A41 11
1378 #define A42 12
1379 #define A43 13
1380 #define A44 14
1381 #define A45 19
1382
1383 #define A50 15
1384 #define A51 16
1385 #define A52 17
1386 #define A53 18
1387 #define A54 19
1388 #define A55 20
1389
1390 void G4ErrorSymMatrix::invert5(G4int& ifail)
1391 {
1392 if(posDefFraction5x5 >= CHOLESKY_THRESHOLD_
1393 {
1394 invertCholesky5(ifail);
1395 posDefFraction5x5 = .9 * posDefFraction5x
1396 if(ifail != 0) // Cholesky failed -- inv
1397 {
1398 invertHaywood5(ifail);
1399 }
1400 }
1401 else
1402 {
1403 if(posDefFraction5x5 + adjustment5x5 >= C
1404 {
1405 invertCholesky5(ifail);
1406 posDefFraction5x5 = .9 * posDefFraction
1407 if(ifail != 0) // Cholesky failed -- i
1408 {
1409 invertHaywood5(ifail);
1410 adjustment5x5 = 0;
1411 }
1412 }
1413 else
1414 {
1415 invertHaywood5(ifail);
1416 adjustment5x5 += CHOLESKY_CREEP_5x5;
1417 }
1418 }
1419 return;
1420 }
1421
1422 void G4ErrorSymMatrix::invert6(G4int& ifail)
1423 {
1424 if(posDefFraction6x6 >= CHOLESKY_THRESHOLD_
1425 {
1426 invertCholesky6(ifail);
1427 posDefFraction6x6 = .9 * posDefFraction6x
1428 if(ifail != 0) // Cholesky failed -- inv
1429 {
1430 invertHaywood6(ifail);
1431 }
1432 }
1433 else
1434 {
1435 if(posDefFraction6x6 + adjustment6x6 >= C
1436 {
1437 invertCholesky6(ifail);
1438 posDefFraction6x6 = .9 * posDefFraction
1439 if(ifail != 0) // Cholesky failed -- i
1440 {
1441 invertHaywood6(ifail);
1442 adjustment6x6 = 0;
1443 }
1444 }
1445 else
1446 {
1447 invertHaywood6(ifail);
1448 adjustment6x6 += CHOLESKY_CREEP_6x6;
1449 }
1450 }
1451 return;
1452 }
1453
1454 void G4ErrorSymMatrix::invertHaywood5(G4int&
1455 {
1456 ifail = 0;
1457
1458 // Find all NECESSARY 2x2 dets: (25 of the
1459
1460 G4double Det2_23_01 = m[A20] * m[A31] - m[A
1461 G4double Det2_23_02 = m[A20] * m[A32] - m[A
1462 G4double Det2_23_03 = m[A20] * m[A33] - m[A
1463 G4double Det2_23_12 = m[A21] * m[A32] - m[A
1464 G4double Det2_23_13 = m[A21] * m[A33] - m[A
1465 G4double Det2_23_23 = m[A22] * m[A33] - m[A
1466 G4double Det2_24_01 = m[A20] * m[A41] - m[A
1467 G4double Det2_24_02 = m[A20] * m[A42] - m[A
1468 G4double Det2_24_03 = m[A20] * m[A43] - m[A
1469 G4double Det2_24_04 = m[A20] * m[A44] - m[A
1470 G4double Det2_24_12 = m[A21] * m[A42] - m[A
1471 G4double Det2_24_13 = m[A21] * m[A43] - m[A
1472 G4double Det2_24_14 = m[A21] * m[A44] - m[A
1473 G4double Det2_24_23 = m[A22] * m[A43] - m[A
1474 G4double Det2_24_24 = m[A22] * m[A44] - m[A
1475 G4double Det2_34_01 = m[A30] * m[A41] - m[A
1476 G4double Det2_34_02 = m[A30] * m[A42] - m[A
1477 G4double Det2_34_03 = m[A30] * m[A43] - m[A
1478 G4double Det2_34_04 = m[A30] * m[A44] - m[A
1479 G4double Det2_34_12 = m[A31] * m[A42] - m[A
1480 G4double Det2_34_13 = m[A31] * m[A43] - m[A
1481 G4double Det2_34_14 = m[A31] * m[A44] - m[A
1482 G4double Det2_34_23 = m[A32] * m[A43] - m[A
1483 G4double Det2_34_24 = m[A32] * m[A44] - m[A
1484 G4double Det2_34_34 = m[A33] * m[A44] - m[A
1485
1486 // Find all NECESSARY 3x3 dets: (30 of th
1487
1488 G4double Det3_123_012 =
1489 m[A10] * Det2_23_12 - m[A11] * Det2_23_02
1490 G4double Det3_123_013 =
1491 m[A10] * Det2_23_13 - m[A11] * Det2_23_03
1492 G4double Det3_123_023 =
1493 m[A10] * Det2_23_23 - m[A12] * Det2_23_03
1494 G4double Det3_123_123 =
1495 m[A11] * Det2_23_23 - m[A12] * Det2_23_13
1496 G4double Det3_124_012 =
1497 m[A10] * Det2_24_12 - m[A11] * Det2_24_02
1498 G4double Det3_124_013 =
1499 m[A10] * Det2_24_13 - m[A11] * Det2_24_03
1500 G4double Det3_124_014 =
1501 m[A10] * Det2_24_14 - m[A11] * Det2_24_04
1502 G4double Det3_124_023 =
1503 m[A10] * Det2_24_23 - m[A12] * Det2_24_03
1504 G4double Det3_124_024 =
1505 m[A10] * Det2_24_24 - m[A12] * Det2_24_04
1506 G4double Det3_124_123 =
1507 m[A11] * Det2_24_23 - m[A12] * Det2_24_13
1508 G4double Det3_124_124 =
1509 m[A11] * Det2_24_24 - m[A12] * Det2_24_14
1510 G4double Det3_134_012 =
1511 m[A10] * Det2_34_12 - m[A11] * Det2_34_02
1512 G4double Det3_134_013 =
1513 m[A10] * Det2_34_13 - m[A11] * Det2_34_03
1514 G4double Det3_134_014 =
1515 m[A10] * Det2_34_14 - m[A11] * Det2_34_04
1516 G4double Det3_134_023 =
1517 m[A10] * Det2_34_23 - m[A12] * Det2_34_03
1518 G4double Det3_134_024 =
1519 m[A10] * Det2_34_24 - m[A12] * Det2_34_04
1520 G4double Det3_134_034 =
1521 m[A10] * Det2_34_34 - m[A13] * Det2_34_04
1522 G4double Det3_134_123 =
1523 m[A11] * Det2_34_23 - m[A12] * Det2_34_13
1524 G4double Det3_134_124 =
1525 m[A11] * Det2_34_24 - m[A12] * Det2_34_14
1526 G4double Det3_134_134 =
1527 m[A11] * Det2_34_34 - m[A13] * Det2_34_14
1528 G4double Det3_234_012 =
1529 m[A20] * Det2_34_12 - m[A21] * Det2_34_02
1530 G4double Det3_234_013 =
1531 m[A20] * Det2_34_13 - m[A21] * Det2_34_03
1532 G4double Det3_234_014 =
1533 m[A20] * Det2_34_14 - m[A21] * Det2_34_04
1534 G4double Det3_234_023 =
1535 m[A20] * Det2_34_23 - m[A22] * Det2_34_03
1536 G4double Det3_234_024 =
1537 m[A20] * Det2_34_24 - m[A22] * Det2_34_04
1538 G4double Det3_234_034 =
1539 m[A20] * Det2_34_34 - m[A23] * Det2_34_04
1540 G4double Det3_234_123 =
1541 m[A21] * Det2_34_23 - m[A22] * Det2_34_13
1542 G4double Det3_234_124 =
1543 m[A21] * Det2_34_24 - m[A22] * Det2_34_14
1544 G4double Det3_234_134 =
1545 m[A21] * Det2_34_34 - m[A23] * Det2_34_14
1546 G4double Det3_234_234 =
1547 m[A22] * Det2_34_34 - m[A23] * Det2_34_24
1548
1549 // Find all NECESSARY 4x4 dets: (15 of th
1550
1551 G4double Det4_0123_0123 = m[A00] * Det3_123
1552 m[A02] * Det3_123
1553 G4double Det4_0124_0123 = m[A00] * Det3_124
1554 m[A02] * Det3_124
1555 G4double Det4_0124_0124 = m[A00] * Det3_124
1556 m[A02] * Det3_124
1557 G4double Det4_0134_0123 = m[A00] * Det3_134
1558 m[A02] * Det3_134
1559 G4double Det4_0134_0124 = m[A00] * Det3_134
1560 m[A02] * Det3_134
1561 G4double Det4_0134_0134 = m[A00] * Det3_134
1562 m[A03] * Det3_134
1563 G4double Det4_0234_0123 = m[A00] * Det3_234
1564 m[A02] * Det3_234
1565 G4double Det4_0234_0124 = m[A00] * Det3_234
1566 m[A02] * Det3_234
1567 G4double Det4_0234_0134 = m[A00] * Det3_234
1568 m[A03] * Det3_234
1569 G4double Det4_0234_0234 = m[A00] * Det3_234
1570 m[A03] * Det3_234
1571 G4double Det4_1234_0123 = m[A10] * Det3_234
1572 m[A12] * Det3_234
1573 G4double Det4_1234_0124 = m[A10] * Det3_234
1574 m[A12] * Det3_234
1575 G4double Det4_1234_0134 = m[A10] * Det3_234
1576 m[A13] * Det3_234
1577 G4double Det4_1234_0234 = m[A10] * Det3_234
1578 m[A13] * Det3_234
1579 G4double Det4_1234_1234 = m[A11] * Det3_234
1580 m[A13] * Det3_234
1581
1582 // Find the 5x5 det:
1583
1584 G4double det = m[A00] * Det4_1234_1234 - m[
1585 m[A02] * Det4_1234_0134 - m[
1586 m[A04] * Det4_1234_0123;
1587
1588 if(det == 0)
1589 {
1590 ifail = 1;
1591 return;
1592 }
1593
1594 G4double oneOverDet = 1.0 / det;
1595 G4double mn1OverDet = -oneOverDet;
1596
1597 m[A00] = Det4_1234_1234 * oneOverDet;
1598 m[A01] = Det4_1234_0234 * mn1OverDet;
1599 m[A02] = Det4_1234_0134 * oneOverDet;
1600 m[A03] = Det4_1234_0124 * mn1OverDet;
1601 m[A04] = Det4_1234_0123 * oneOverDet;
1602
1603 m[A11] = Det4_0234_0234 * oneOverDet;
1604 m[A12] = Det4_0234_0134 * mn1OverDet;
1605 m[A13] = Det4_0234_0124 * oneOverDet;
1606 m[A14] = Det4_0234_0123 * mn1OverDet;
1607
1608 m[A22] = Det4_0134_0134 * oneOverDet;
1609 m[A23] = Det4_0134_0124 * mn1OverDet;
1610 m[A24] = Det4_0134_0123 * oneOverDet;
1611
1612 m[A33] = Det4_0124_0124 * oneOverDet;
1613 m[A34] = Det4_0124_0123 * mn1OverDet;
1614
1615 m[A44] = Det4_0123_0123 * oneOverDet;
1616
1617 return;
1618 }
1619
1620 void G4ErrorSymMatrix::invertHaywood6(G4int&
1621 {
1622 ifail = 0;
1623
1624 // Find all NECESSARY 2x2 dets: (39 of the
1625
1626 G4double Det2_34_01 = m[A30] * m[A41] - m[A
1627 G4double Det2_34_02 = m[A30] * m[A42] - m[A
1628 G4double Det2_34_03 = m[A30] * m[A43] - m[A
1629 G4double Det2_34_04 = m[A30] * m[A44] - m[A
1630 G4double Det2_34_12 = m[A31] * m[A42] - m[A
1631 G4double Det2_34_13 = m[A31] * m[A43] - m[A
1632 G4double Det2_34_14 = m[A31] * m[A44] - m[A
1633 G4double Det2_34_23 = m[A32] * m[A43] - m[A
1634 G4double Det2_34_24 = m[A32] * m[A44] - m[A
1635 G4double Det2_34_34 = m[A33] * m[A44] - m[A
1636 G4double Det2_35_01 = m[A30] * m[A51] - m[A
1637 G4double Det2_35_02 = m[A30] * m[A52] - m[A
1638 G4double Det2_35_03 = m[A30] * m[A53] - m[A
1639 G4double Det2_35_04 = m[A30] * m[A54] - m[A
1640 G4double Det2_35_05 = m[A30] * m[A55] - m[A
1641 G4double Det2_35_12 = m[A31] * m[A52] - m[A
1642 G4double Det2_35_13 = m[A31] * m[A53] - m[A
1643 G4double Det2_35_14 = m[A31] * m[A54] - m[A
1644 G4double Det2_35_15 = m[A31] * m[A55] - m[A
1645 G4double Det2_35_23 = m[A32] * m[A53] - m[A
1646 G4double Det2_35_24 = m[A32] * m[A54] - m[A
1647 G4double Det2_35_25 = m[A32] * m[A55] - m[A
1648 G4double Det2_35_34 = m[A33] * m[A54] - m[A
1649 G4double Det2_35_35 = m[A33] * m[A55] - m[A
1650 G4double Det2_45_01 = m[A40] * m[A51] - m[A
1651 G4double Det2_45_02 = m[A40] * m[A52] - m[A
1652 G4double Det2_45_03 = m[A40] * m[A53] - m[A
1653 G4double Det2_45_04 = m[A40] * m[A54] - m[A
1654 G4double Det2_45_05 = m[A40] * m[A55] - m[A
1655 G4double Det2_45_12 = m[A41] * m[A52] - m[A
1656 G4double Det2_45_13 = m[A41] * m[A53] - m[A
1657 G4double Det2_45_14 = m[A41] * m[A54] - m[A
1658 G4double Det2_45_15 = m[A41] * m[A55] - m[A
1659 G4double Det2_45_23 = m[A42] * m[A53] - m[A
1660 G4double Det2_45_24 = m[A42] * m[A54] - m[A
1661 G4double Det2_45_25 = m[A42] * m[A55] - m[A
1662 G4double Det2_45_34 = m[A43] * m[A54] - m[A
1663 G4double Det2_45_35 = m[A43] * m[A55] - m[A
1664 G4double Det2_45_45 = m[A44] * m[A55] - m[A
1665
1666 // Find all NECESSARY 3x3 dets: (65 of the
1667
1668 G4double Det3_234_012 =
1669 m[A20] * Det2_34_12 - m[A21] * Det2_34_02
1670 G4double Det3_234_013 =
1671 m[A20] * Det2_34_13 - m[A21] * Det2_34_03
1672 G4double Det3_234_014 =
1673 m[A20] * Det2_34_14 - m[A21] * Det2_34_04
1674 G4double Det3_234_023 =
1675 m[A20] * Det2_34_23 - m[A22] * Det2_34_03
1676 G4double Det3_234_024 =
1677 m[A20] * Det2_34_24 - m[A22] * Det2_34_04
1678 G4double Det3_234_034 =
1679 m[A20] * Det2_34_34 - m[A23] * Det2_34_04
1680 G4double Det3_234_123 =
1681 m[A21] * Det2_34_23 - m[A22] * Det2_34_13
1682 G4double Det3_234_124 =
1683 m[A21] * Det2_34_24 - m[A22] * Det2_34_14
1684 G4double Det3_234_134 =
1685 m[A21] * Det2_34_34 - m[A23] * Det2_34_14
1686 G4double Det3_234_234 =
1687 m[A22] * Det2_34_34 - m[A23] * Det2_34_24
1688 G4double Det3_235_012 =
1689 m[A20] * Det2_35_12 - m[A21] * Det2_35_02
1690 G4double Det3_235_013 =
1691 m[A20] * Det2_35_13 - m[A21] * Det2_35_03
1692 G4double Det3_235_014 =
1693 m[A20] * Det2_35_14 - m[A21] * Det2_35_04
1694 G4double Det3_235_015 =
1695 m[A20] * Det2_35_15 - m[A21] * Det2_35_05
1696 G4double Det3_235_023 =
1697 m[A20] * Det2_35_23 - m[A22] * Det2_35_03
1698 G4double Det3_235_024 =
1699 m[A20] * Det2_35_24 - m[A22] * Det2_35_04
1700 G4double Det3_235_025 =
1701 m[A20] * Det2_35_25 - m[A22] * Det2_35_05
1702 G4double Det3_235_034 =
1703 m[A20] * Det2_35_34 - m[A23] * Det2_35_04
1704 G4double Det3_235_035 =
1705 m[A20] * Det2_35_35 - m[A23] * Det2_35_05
1706 G4double Det3_235_123 =
1707 m[A21] * Det2_35_23 - m[A22] * Det2_35_13
1708 G4double Det3_235_124 =
1709 m[A21] * Det2_35_24 - m[A22] * Det2_35_14
1710 G4double Det3_235_125 =
1711 m[A21] * Det2_35_25 - m[A22] * Det2_35_15
1712 G4double Det3_235_134 =
1713 m[A21] * Det2_35_34 - m[A23] * Det2_35_14
1714 G4double Det3_235_135 =
1715 m[A21] * Det2_35_35 - m[A23] * Det2_35_15
1716 G4double Det3_235_234 =
1717 m[A22] * Det2_35_34 - m[A23] * Det2_35_24
1718 G4double Det3_235_235 =
1719 m[A22] * Det2_35_35 - m[A23] * Det2_35_25
1720 G4double Det3_245_012 =
1721 m[A20] * Det2_45_12 - m[A21] * Det2_45_02
1722 G4double Det3_245_013 =
1723 m[A20] * Det2_45_13 - m[A21] * Det2_45_03
1724 G4double Det3_245_014 =
1725 m[A20] * Det2_45_14 - m[A21] * Det2_45_04
1726 G4double Det3_245_015 =
1727 m[A20] * Det2_45_15 - m[A21] * Det2_45_05
1728 G4double Det3_245_023 =
1729 m[A20] * Det2_45_23 - m[A22] * Det2_45_03
1730 G4double Det3_245_024 =
1731 m[A20] * Det2_45_24 - m[A22] * Det2_45_04
1732 G4double Det3_245_025 =
1733 m[A20] * Det2_45_25 - m[A22] * Det2_45_05
1734 G4double Det3_245_034 =
1735 m[A20] * Det2_45_34 - m[A23] * Det2_45_04
1736 G4double Det3_245_035 =
1737 m[A20] * Det2_45_35 - m[A23] * Det2_45_05
1738 G4double Det3_245_045 =
1739 m[A20] * Det2_45_45 - m[A24] * Det2_45_05
1740 G4double Det3_245_123 =
1741 m[A21] * Det2_45_23 - m[A22] * Det2_45_13
1742 G4double Det3_245_124 =
1743 m[A21] * Det2_45_24 - m[A22] * Det2_45_14
1744 G4double Det3_245_125 =
1745 m[A21] * Det2_45_25 - m[A22] * Det2_45_15
1746 G4double Det3_245_134 =
1747 m[A21] * Det2_45_34 - m[A23] * Det2_45_14
1748 G4double Det3_245_135 =
1749 m[A21] * Det2_45_35 - m[A23] * Det2_45_15
1750 G4double Det3_245_145 =
1751 m[A21] * Det2_45_45 - m[A24] * Det2_45_15
1752 G4double Det3_245_234 =
1753 m[A22] * Det2_45_34 - m[A23] * Det2_45_24
1754 G4double Det3_245_235 =
1755 m[A22] * Det2_45_35 - m[A23] * Det2_45_25
1756 G4double Det3_245_245 =
1757 m[A22] * Det2_45_45 - m[A24] * Det2_45_25
1758 G4double Det3_345_012 =
1759 m[A30] * Det2_45_12 - m[A31] * Det2_45_02
1760 G4double Det3_345_013 =
1761 m[A30] * Det2_45_13 - m[A31] * Det2_45_03
1762 G4double Det3_345_014 =
1763 m[A30] * Det2_45_14 - m[A31] * Det2_45_04
1764 G4double Det3_345_015 =
1765 m[A30] * Det2_45_15 - m[A31] * Det2_45_05
1766 G4double Det3_345_023 =
1767 m[A30] * Det2_45_23 - m[A32] * Det2_45_03
1768 G4double Det3_345_024 =
1769 m[A30] * Det2_45_24 - m[A32] * Det2_45_04
1770 G4double Det3_345_025 =
1771 m[A30] * Det2_45_25 - m[A32] * Det2_45_05
1772 G4double Det3_345_034 =
1773 m[A30] * Det2_45_34 - m[A33] * Det2_45_04
1774 G4double Det3_345_035 =
1775 m[A30] * Det2_45_35 - m[A33] * Det2_45_05
1776 G4double Det3_345_045 =
1777 m[A30] * Det2_45_45 - m[A34] * Det2_45_05
1778 G4double Det3_345_123 =
1779 m[A31] * Det2_45_23 - m[A32] * Det2_45_13
1780 G4double Det3_345_124 =
1781 m[A31] * Det2_45_24 - m[A32] * Det2_45_14
1782 G4double Det3_345_125 =
1783 m[A31] * Det2_45_25 - m[A32] * Det2_45_15
1784 G4double Det3_345_134 =
1785 m[A31] * Det2_45_34 - m[A33] * Det2_45_14
1786 G4double Det3_345_135 =
1787 m[A31] * Det2_45_35 - m[A33] * Det2_45_15
1788 G4double Det3_345_145 =
1789 m[A31] * Det2_45_45 - m[A34] * Det2_45_15
1790 G4double Det3_345_234 =
1791 m[A32] * Det2_45_34 - m[A33] * Det2_45_24
1792 G4double Det3_345_235 =
1793 m[A32] * Det2_45_35 - m[A33] * Det2_45_25
1794 G4double Det3_345_245 =
1795 m[A32] * Det2_45_45 - m[A34] * Det2_45_25
1796 G4double Det3_345_345 =
1797 m[A33] * Det2_45_45 - m[A34] * Det2_45_35
1798
1799 // Find all NECESSARY 4x4 dets: (55 of the
1800
1801 G4double Det4_1234_0123 = m[A10] * Det3_234
1802 m[A12] * Det3_234
1803 G4double Det4_1234_0124 = m[A10] * Det3_234
1804 m[A12] * Det3_234
1805 G4double Det4_1234_0134 = m[A10] * Det3_234
1806 m[A13] * Det3_234
1807 G4double Det4_1234_0234 = m[A10] * Det3_234
1808 m[A13] * Det3_234
1809 G4double Det4_1234_1234 = m[A11] * Det3_234
1810 m[A13] * Det3_234
1811 G4double Det4_1235_0123 = m[A10] * Det3_235
1812 m[A12] * Det3_235
1813 G4double Det4_1235_0124 = m[A10] * Det3_235
1814 m[A12] * Det3_235
1815 G4double Det4_1235_0125 = m[A10] * Det3_235
1816 m[A12] * Det3_235
1817 G4double Det4_1235_0134 = m[A10] * Det3_235
1818 m[A13] * Det3_235
1819 G4double Det4_1235_0135 = m[A10] * Det3_235
1820 m[A13] * Det3_235
1821 G4double Det4_1235_0234 = m[A10] * Det3_235
1822 m[A13] * Det3_235
1823 G4double Det4_1235_0235 = m[A10] * Det3_235
1824 m[A13] * Det3_235
1825 G4double Det4_1235_1234 = m[A11] * Det3_235
1826 m[A13] * Det3_235
1827 G4double Det4_1235_1235 = m[A11] * Det3_235
1828 m[A13] * Det3_235
1829 G4double Det4_1245_0123 = m[A10] * Det3_245
1830 m[A12] * Det3_245
1831 G4double Det4_1245_0124 = m[A10] * Det3_245
1832 m[A12] * Det3_245
1833 G4double Det4_1245_0125 = m[A10] * Det3_245
1834 m[A12] * Det3_245
1835 G4double Det4_1245_0134 = m[A10] * Det3_245
1836 m[A13] * Det3_245
1837 G4double Det4_1245_0135 = m[A10] * Det3_245
1838 m[A13] * Det3_245
1839 G4double Det4_1245_0145 = m[A10] * Det3_245
1840 m[A14] * Det3_245
1841 G4double Det4_1245_0234 = m[A10] * Det3_245
1842 m[A13] * Det3_245
1843 G4double Det4_1245_0235 = m[A10] * Det3_245
1844 m[A13] * Det3_245
1845 G4double Det4_1245_0245 = m[A10] * Det3_245
1846 m[A14] * Det3_245
1847 G4double Det4_1245_1234 = m[A11] * Det3_245
1848 m[A13] * Det3_245
1849 G4double Det4_1245_1235 = m[A11] * Det3_245
1850 m[A13] * Det3_245
1851 G4double Det4_1245_1245 = m[A11] * Det3_245
1852 m[A14] * Det3_245
1853 G4double Det4_1345_0123 = m[A10] * Det3_345
1854 m[A12] * Det3_345
1855 G4double Det4_1345_0124 = m[A10] * Det3_345
1856 m[A12] * Det3_345
1857 G4double Det4_1345_0125 = m[A10] * Det3_345
1858 m[A12] * Det3_345
1859 G4double Det4_1345_0134 = m[A10] * Det3_345
1860 m[A13] * Det3_345
1861 G4double Det4_1345_0135 = m[A10] * Det3_345
1862 m[A13] * Det3_345
1863 G4double Det4_1345_0145 = m[A10] * Det3_345
1864 m[A14] * Det3_345
1865 G4double Det4_1345_0234 = m[A10] * Det3_345
1866 m[A13] * Det3_345
1867 G4double Det4_1345_0235 = m[A10] * Det3_345
1868 m[A13] * Det3_345
1869 G4double Det4_1345_0245 = m[A10] * Det3_345
1870 m[A14] * Det3_345
1871 G4double Det4_1345_0345 = m[A10] * Det3_345
1872 m[A14] * Det3_345
1873 G4double Det4_1345_1234 = m[A11] * Det3_345
1874 m[A13] * Det3_345
1875 G4double Det4_1345_1235 = m[A11] * Det3_345
1876 m[A13] * Det3_345
1877 G4double Det4_1345_1245 = m[A11] * Det3_345
1878 m[A14] * Det3_345
1879 G4double Det4_1345_1345 = m[A11] * Det3_345
1880 m[A14] * Det3_345
1881 G4double Det4_2345_0123 = m[A20] * Det3_345
1882 m[A22] * Det3_345
1883 G4double Det4_2345_0124 = m[A20] * Det3_345
1884 m[A22] * Det3_345
1885 G4double Det4_2345_0125 = m[A20] * Det3_345
1886 m[A22] * Det3_345
1887 G4double Det4_2345_0134 = m[A20] * Det3_345
1888 m[A23] * Det3_345
1889 G4double Det4_2345_0135 = m[A20] * Det3_345
1890 m[A23] * Det3_345
1891 G4double Det4_2345_0145 = m[A20] * Det3_345
1892 m[A24] * Det3_345
1893 G4double Det4_2345_0234 = m[A20] * Det3_345
1894 m[A23] * Det3_345
1895 G4double Det4_2345_0235 = m[A20] * Det3_345
1896 m[A23] * Det3_345
1897 G4double Det4_2345_0245 = m[A20] * Det3_345
1898 m[A24] * Det3_345
1899 G4double Det4_2345_0345 = m[A20] * Det3_345
1900 m[A24] * Det3_345
1901 G4double Det4_2345_1234 = m[A21] * Det3_345
1902 m[A23] * Det3_345
1903 G4double Det4_2345_1235 = m[A21] * Det3_345
1904 m[A23] * Det3_345
1905 G4double Det4_2345_1245 = m[A21] * Det3_345
1906 m[A24] * Det3_345
1907 G4double Det4_2345_1345 = m[A21] * Det3_345
1908 m[A24] * Det3_345
1909 G4double Det4_2345_2345 = m[A22] * Det3_345
1910 m[A24] * Det3_345
1911
1912 // Find all NECESSARY 5x5 dets: (19 of the
1913
1914 G4double Det5_01234_01234 =
1915 m[A00] * Det4_1234_1234 - m[A01] * Det4_1
1916 m[A02] * Det4_1234_0134 - m[A03] * Det4_1
1917 G4double Det5_01235_01234 =
1918 m[A00] * Det4_1235_1234 - m[A01] * Det4_1
1919 m[A02] * Det4_1235_0134 - m[A03] * Det4_1
1920 G4double Det5_01235_01235 =
1921 m[A00] * Det4_1235_1235 - m[A01] * Det4_1
1922 m[A02] * Det4_1235_0135 - m[A03] * Det4_1
1923 G4double Det5_01245_01234 =
1924 m[A00] * Det4_1245_1234 - m[A01] * Det4_1
1925 m[A02] * Det4_1245_0134 - m[A03] * Det4_1
1926 G4double Det5_01245_01235 =
1927 m[A00] * Det4_1245_1235 - m[A01] * Det4_1
1928 m[A02] * Det4_1245_0135 - m[A03] * Det4_1
1929 G4double Det5_01245_01245 =
1930 m[A00] * Det4_1245_1245 - m[A01] * Det4_1
1931 m[A02] * Det4_1245_0145 - m[A04] * Det4_1
1932 G4double Det5_01345_01234 =
1933 m[A00] * Det4_1345_1234 - m[A01] * Det4_1
1934 m[A02] * Det4_1345_0134 - m[A03] * Det4_1
1935 G4double Det5_01345_01235 =
1936 m[A00] * Det4_1345_1235 - m[A01] * Det4_1
1937 m[A02] * Det4_1345_0135 - m[A03] * Det4_1
1938 G4double Det5_01345_01245 =
1939 m[A00] * Det4_1345_1245 - m[A01] * Det4_1
1940 m[A02] * Det4_1345_0145 - m[A04] * Det4_1
1941 G4double Det5_01345_01345 =
1942 m[A00] * Det4_1345_1345 - m[A01] * Det4_1
1943 m[A03] * Det4_1345_0145 - m[A04] * Det4_1
1944 G4double Det5_02345_01234 =
1945 m[A00] * Det4_2345_1234 - m[A01] * Det4_2
1946 m[A02] * Det4_2345_0134 - m[A03] * Det4_2
1947 G4double Det5_02345_01235 =
1948 m[A00] * Det4_2345_1235 - m[A01] * Det4_2
1949 m[A02] * Det4_2345_0135 - m[A03] * Det4_2
1950 G4double Det5_02345_01245 =
1951 m[A00] * Det4_2345_1245 - m[A01] * Det4_2
1952 m[A02] * Det4_2345_0145 - m[A04] * Det4_2
1953 G4double Det5_02345_01345 =
1954 m[A00] * Det4_2345_1345 - m[A01] * Det4_2
1955 m[A03] * Det4_2345_0145 - m[A04] * Det4_2
1956 G4double Det5_02345_02345 =
1957 m[A00] * Det4_2345_2345 - m[A02] * Det4_2
1958 m[A03] * Det4_2345_0245 - m[A04] * Det4_2
1959 G4double Det5_12345_01234 =
1960 m[A10] * Det4_2345_1234 - m[A11] * Det4_2
1961 m[A12] * Det4_2345_0134 - m[A13] * Det4_2
1962 G4double Det5_12345_01235 =
1963 m[A10] * Det4_2345_1235 - m[A11] * Det4_2
1964 m[A12] * Det4_2345_0135 - m[A13] * Det4_2
1965 G4double Det5_12345_01245 =
1966 m[A10] * Det4_2345_1245 - m[A11] * Det4_2
1967 m[A12] * Det4_2345_0145 - m[A14] * Det4_2
1968 G4double Det5_12345_01345 =
1969 m[A10] * Det4_2345_1345 - m[A11] * Det4_2
1970 m[A13] * Det4_2345_0145 - m[A14] * Det4_2
1971 G4double Det5_12345_02345 =
1972 m[A10] * Det4_2345_2345 - m[A12] * Det4_2
1973 m[A13] * Det4_2345_0245 - m[A14] * Det4_2
1974 G4double Det5_12345_12345 =
1975 m[A11] * Det4_2345_2345 - m[A12] * Det4_2
1976 m[A13] * Det4_2345_1245 - m[A14] * Det4_2
1977
1978 // Find the determinant
1979
1980 G4double det = m[A00] * Det5_12345_12345 -
1981 m[A02] * Det5_12345_01345 -
1982 m[A04] * Det5_12345_01235 -
1983
1984 if(det == 0)
1985 {
1986 ifail = 1;
1987 return;
1988 }
1989
1990 G4double oneOverDet = 1.0 / det;
1991 G4double mn1OverDet = -oneOverDet;
1992
1993 m[A00] = Det5_12345_12345 * oneOverDet;
1994 m[A01] = Det5_12345_02345 * mn1OverDet;
1995 m[A02] = Det5_12345_01345 * oneOverDet;
1996 m[A03] = Det5_12345_01245 * mn1OverDet;
1997 m[A04] = Det5_12345_01235 * oneOverDet;
1998 m[A05] = Det5_12345_01234 * mn1OverDet;
1999
2000 m[A11] = Det5_02345_02345 * oneOverDet;
2001 m[A12] = Det5_02345_01345 * mn1OverDet;
2002 m[A13] = Det5_02345_01245 * oneOverDet;
2003 m[A14] = Det5_02345_01235 * mn1OverDet;
2004 m[A15] = Det5_02345_01234 * oneOverDet;
2005
2006 m[A22] = Det5_01345_01345 * oneOverDet;
2007 m[A23] = Det5_01345_01245 * mn1OverDet;
2008 m[A24] = Det5_01345_01235 * oneOverDet;
2009 m[A25] = Det5_01345_01234 * mn1OverDet;
2010
2011 m[A33] = Det5_01245_01245 * oneOverDet;
2012 m[A34] = Det5_01245_01235 * mn1OverDet;
2013 m[A35] = Det5_01245_01234 * oneOverDet;
2014
2015 m[A44] = Det5_01235_01235 * oneOverDet;
2016 m[A45] = Det5_01235_01234 * mn1OverDet;
2017
2018 m[A55] = Det5_01234_01234 * oneOverDet;
2019
2020 return;
2021 }
2022
2023 void G4ErrorSymMatrix::invertCholesky5(G4int&
2024 {
2025 // Invert by
2026 //
2027 // a) decomposing M = G*G^T with G lower tr
2028 // (if M is not positive definite this w
2029 // b) inverting G to form H
2030 // c) multiplying H^T * H to get M^-1.
2031 //
2032 // If the matrix is pos. def. it is inverte
2033 // If the matrix is not pos. def. it remain
2034
2035 G4double h10; // below-diagonal elements o
2036 G4double h20, h21;
2037 G4double h30, h31, h32;
2038 G4double h40, h41, h42, h43;
2039
2040 G4double h00, h11, h22, h33, h44; // 1/dia
2041 // diago
2042
2043 G4double g10; // below-diagonal elements o
2044 G4double g20, g21;
2045 G4double g30, g31, g32;
2046 G4double g40, g41, g42, g43;
2047
2048 ifail = 1; // We start by assuing we won't
2049
2050 // Form G -- compute diagonal members of H
2051 //-------
2052
2053 // Scale first column by 1/sqrt(A00)
2054
2055 h00 = m[A00];
2056 if(h00 <= 0)
2057 {
2058 return;
2059 }
2060 h00 = 1.0 / std::sqrt(h00);
2061
2062 g10 = m[A10] * h00;
2063 g20 = m[A20] * h00;
2064 g30 = m[A30] * h00;
2065 g40 = m[A40] * h00;
2066
2067 // Form G11 (actually, just h11)
2068
2069 h11 = m[A11] - (g10 * g10);
2070 if(h11 <= 0)
2071 {
2072 return;
2073 }
2074 h11 = 1.0 / std::sqrt(h11);
2075
2076 // Subtract inter-column column dot product
2077 // scale to get column 1 of G
2078
2079 g21 = (m[A21] - (g10 * g20)) * h11;
2080 g31 = (m[A31] - (g10 * g30)) * h11;
2081 g41 = (m[A41] - (g10 * g40)) * h11;
2082
2083 // Form G22 (actually, just h22)
2084
2085 h22 = m[A22] - (g20 * g20) - (g21 * g21);
2086 if(h22 <= 0)
2087 {
2088 return;
2089 }
2090 h22 = 1.0 / std::sqrt(h22);
2091
2092 // Subtract inter-column column dot product
2093 // scale to get column 2 of G
2094
2095 g32 = (m[A32] - (g20 * g30) - (g21 * g31))
2096 g42 = (m[A42] - (g20 * g40) - (g21 * g41))
2097
2098 // Form G33 (actually, just h33)
2099
2100 h33 = m[A33] - (g30 * g30) - (g31 * g31) -
2101 if(h33 <= 0)
2102 {
2103 return;
2104 }
2105 h33 = 1.0 / std::sqrt(h33);
2106
2107 // Subtract inter-column column dot product
2108
2109 g43 = (m[A43] - (g30 * g40) - (g31 * g41) -
2110
2111 // Finally form h44 - if this is possible i
2112
2113 h44 = m[A44] - (g40 * g40) - (g41 * g41) -
2114 if(h44 <= 0)
2115 {
2116 return;
2117 }
2118 h44 = 1.0 / std::sqrt(h44);
2119
2120 // Form H = 1/G -- diagonal members of H ar
2121 //-------------
2122
2123 // The order here is dictated by speed cons
2124
2125 h43 = -h33 * g43 * h44;
2126 h32 = -h22 * g32 * h33;
2127 h42 = -h22 * (g32 * h43 + g42 * h44);
2128 h21 = -h11 * g21 * h22;
2129 h31 = -h11 * (g21 * h32 + g31 * h33);
2130 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 *
2131 h10 = -h00 * g10 * h11;
2132 h20 = -h00 * (g10 * h21 + g20 * h22);
2133 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 *
2134 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 *
2135
2136 // Change this to its inverse = H^T*H
2137 //------------------------------------
2138
2139 m[A00] = h00 * h00 + h10 * h10 + h20 * h20
2140 m[A01] = h10 * h11 + h20 * h21 + h30 * h31
2141 m[A11] = h11 * h11 + h21 * h21 + h31 * h31
2142 m[A02] = h20 * h22 + h30 * h32 + h40 * h42;
2143 m[A12] = h21 * h22 + h31 * h32 + h41 * h42;
2144 m[A22] = h22 * h22 + h32 * h32 + h42 * h42;
2145 m[A03] = h30 * h33 + h40 * h43;
2146 m[A13] = h31 * h33 + h41 * h43;
2147 m[A23] = h32 * h33 + h42 * h43;
2148 m[A33] = h33 * h33 + h43 * h43;
2149 m[A04] = h40 * h44;
2150 m[A14] = h41 * h44;
2151 m[A24] = h42 * h44;
2152 m[A34] = h43 * h44;
2153 m[A44] = h44 * h44;
2154
2155 ifail = 0;
2156 return;
2157 }
2158
2159 void G4ErrorSymMatrix::invertCholesky6(G4int&
2160 {
2161 // Invert by
2162 //
2163 // a) decomposing M = G*G^T with G lower tr
2164 // (if M is not positive definite this w
2165 // b) inverting G to form H
2166 // c) multiplying H^T * H to get M^-1.
2167 //
2168 // If the matrix is pos. def. it is inverte
2169 // If the matrix is not pos. def. it remain
2170
2171 G4double h10; // below-diagonal elements o
2172 G4double h20, h21;
2173 G4double h30, h31, h32;
2174 G4double h40, h41, h42, h43;
2175 G4double h50, h51, h52, h53, h54;
2176
2177 G4double h00, h11, h22, h33, h44, h55; //
2178 //
2179
2180 G4double g10; // below-diagonal elements o
2181 G4double g20, g21;
2182 G4double g30, g31, g32;
2183 G4double g40, g41, g42, g43;
2184 G4double g50, g51, g52, g53, g54;
2185
2186 ifail = 1; // We start by assuing we won't
2187
2188 // Form G -- compute diagonal members of H
2189 //-------
2190
2191 // Scale first column by 1/sqrt(A00)
2192
2193 h00 = m[A00];
2194 if(h00 <= 0)
2195 {
2196 return;
2197 }
2198 h00 = 1.0 / std::sqrt(h00);
2199
2200 g10 = m[A10] * h00;
2201 g20 = m[A20] * h00;
2202 g30 = m[A30] * h00;
2203 g40 = m[A40] * h00;
2204 g50 = m[A50] * h00;
2205
2206 // Form G11 (actually, just h11)
2207
2208 h11 = m[A11] - (g10 * g10);
2209 if(h11 <= 0)
2210 {
2211 return;
2212 }
2213 h11 = 1.0 / std::sqrt(h11);
2214
2215 // Subtract inter-column column dot product
2216 // scale to get column 1 of G
2217
2218 g21 = (m[A21] - (g10 * g20)) * h11;
2219 g31 = (m[A31] - (g10 * g30)) * h11;
2220 g41 = (m[A41] - (g10 * g40)) * h11;
2221 g51 = (m[A51] - (g10 * g50)) * h11;
2222
2223 // Form G22 (actually, just h22)
2224
2225 h22 = m[A22] - (g20 * g20) - (g21 * g21);
2226 if(h22 <= 0)
2227 {
2228 return;
2229 }
2230 h22 = 1.0 / std::sqrt(h22);
2231
2232 // Subtract inter-column column dot product
2233 // scale to get column 2 of G
2234
2235 g32 = (m[A32] - (g20 * g30) - (g21 * g31))
2236 g42 = (m[A42] - (g20 * g40) - (g21 * g41))
2237 g52 = (m[A52] - (g20 * g50) - (g21 * g51))
2238
2239 // Form G33 (actually, just h33)
2240
2241 h33 = m[A33] - (g30 * g30) - (g31 * g31) -
2242 if(h33 <= 0)
2243 {
2244 return;
2245 }
2246 h33 = 1.0 / std::sqrt(h33);
2247
2248 // Subtract inter-column column dot product
2249 // scale to get column 3 of G
2250
2251 g43 = (m[A43] - (g30 * g40) - (g31 * g41) -
2252 g53 = (m[A53] - (g30 * g50) - (g31 * g51) -
2253
2254 // Form G44 (actually, just h44)
2255
2256 h44 = m[A44] - (g40 * g40) - (g41 * g41) -
2257 if(h44 <= 0)
2258 {
2259 return;
2260 }
2261 h44 = 1.0 / std::sqrt(h44);
2262
2263 // Subtract inter-column column dot product
2264
2265 g54 = (m[A54] - (g40 * g50) - (g41 * g51) -
2266
2267 // Finally form h55 - if this is possible i
2268
2269 h55 = m[A55] - (g50 * g50) - (g51 * g51) -
2270 (g54 * g54);
2271 if(h55 <= 0)
2272 {
2273 return;
2274 }
2275 h55 = 1.0 / std::sqrt(h55);
2276
2277 // Form H = 1/G -- diagonal members of H ar
2278 //-------------
2279
2280 // The order here is dictated by speed cons
2281
2282 h54 = -h44 * g54 * h55;
2283 h43 = -h33 * g43 * h44;
2284 h53 = -h33 * (g43 * h54 + g53 * h55);
2285 h32 = -h22 * g32 * h33;
2286 h42 = -h22 * (g32 * h43 + g42 * h44);
2287 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 *
2288 h21 = -h11 * g21 * h22;
2289 h31 = -h11 * (g21 * h32 + g31 * h33);
2290 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 *
2291 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 *
2292 h10 = -h00 * g10 * h11;
2293 h20 = -h00 * (g10 * h21 + g20 * h22);
2294 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 *
2295 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 *
2296 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 *
2297
2298 // Change this to its inverse = H^T*H
2299 //------------------------------------
2300
2301 m[A00] =
2302 h00 * h00 + h10 * h10 + h20 * h20 + h30 *
2303 m[A01] = h10 * h11 + h20 * h21 + h30 * h31
2304 m[A11] = h11 * h11 + h21 * h21 + h31 * h31
2305 m[A02] = h20 * h22 + h30 * h32 + h40 * h42
2306 m[A12] = h21 * h22 + h31 * h32 + h41 * h42
2307 m[A22] = h22 * h22 + h32 * h32 + h42 * h42
2308 m[A03] = h30 * h33 + h40 * h43 + h50 * h53;
2309 m[A13] = h31 * h33 + h41 * h43 + h51 * h53;
2310 m[A23] = h32 * h33 + h42 * h43 + h52 * h53;
2311 m[A33] = h33 * h33 + h43 * h43 + h53 * h53;
2312 m[A04] = h40 * h44 + h50 * h54;
2313 m[A14] = h41 * h44 + h51 * h54;
2314 m[A24] = h42 * h44 + h52 * h54;
2315 m[A34] = h43 * h44 + h53 * h54;
2316 m[A44] = h44 * h44 + h54 * h54;
2317 m[A05] = h50 * h55;
2318 m[A15] = h51 * h55;
2319 m[A25] = h52 * h55;
2320 m[A35] = h53 * h55;
2321 m[A45] = h54 * h55;
2322 m[A55] = h55 * h55;
2323
2324 ifail = 0;
2325 return;
2326 }
2327
2328 void G4ErrorSymMatrix::invert4(G4int& ifail)
2329 {
2330 ifail = 0;
2331
2332 // Find all NECESSARY 2x2 dets: (14 of the
2333
2334 G4double Det2_12_01 = m[A10] * m[A21] - m[A
2335 G4double Det2_12_02 = m[A10] * m[A22] - m[A
2336 G4double Det2_12_12 = m[A11] * m[A22] - m[A
2337 G4double Det2_13_01 = m[A10] * m[A31] - m[A
2338 G4double Det2_13_02 = m[A10] * m[A32] - m[A
2339 G4double Det2_13_03 = m[A10] * m[A33] - m[A
2340 G4double Det2_13_12 = m[A11] * m[A32] - m[A
2341 G4double Det2_13_13 = m[A11] * m[A33] - m[A
2342 G4double Det2_23_01 = m[A20] * m[A31] - m[A
2343 G4double Det2_23_02 = m[A20] * m[A32] - m[A
2344 G4double Det2_23_03 = m[A20] * m[A33] - m[A
2345 G4double Det2_23_12 = m[A21] * m[A32] - m[A
2346 G4double Det2_23_13 = m[A21] * m[A33] - m[A
2347 G4double Det2_23_23 = m[A22] * m[A33] - m[A
2348
2349 // Find all NECESSARY 3x3 dets: (10 of th
2350
2351 G4double Det3_012_012 =
2352 m[A00] * Det2_12_12 - m[A01] * Det2_12_02
2353 G4double Det3_013_012 =
2354 m[A00] * Det2_13_12 - m[A01] * Det2_13_02
2355 G4double Det3_013_013 =
2356 m[A00] * Det2_13_13 - m[A01] * Det2_13_03
2357 G4double Det3_023_012 =
2358 m[A00] * Det2_23_12 - m[A01] * Det2_23_02
2359 G4double Det3_023_013 =
2360 m[A00] * Det2_23_13 - m[A01] * Det2_23_03
2361 G4double Det3_023_023 =
2362 m[A00] * Det2_23_23 - m[A02] * Det2_23_03
2363 G4double Det3_123_012 =
2364 m[A10] * Det2_23_12 - m[A11] * Det2_23_02
2365 G4double Det3_123_013 =
2366 m[A10] * Det2_23_13 - m[A11] * Det2_23_03
2367 G4double Det3_123_023 =
2368 m[A10] * Det2_23_23 - m[A12] * Det2_23_03
2369 G4double Det3_123_123 =
2370 m[A11] * Det2_23_23 - m[A12] * Det2_23_13
2371
2372 // Find the 4x4 det:
2373
2374 G4double det = m[A00] * Det3_123_123 - m[A0
2375 m[A02] * Det3_123_013 - m[A0
2376
2377 if(det == 0)
2378 {
2379 ifail = 1;
2380 return;
2381 }
2382
2383 G4double oneOverDet = 1.0 / det;
2384 G4double mn1OverDet = -oneOverDet;
2385
2386 m[A00] = Det3_123_123 * oneOverDet;
2387 m[A01] = Det3_123_023 * mn1OverDet;
2388 m[A02] = Det3_123_013 * oneOverDet;
2389 m[A03] = Det3_123_012 * mn1OverDet;
2390
2391 m[A11] = Det3_023_023 * oneOverDet;
2392 m[A12] = Det3_023_013 * mn1OverDet;
2393 m[A13] = Det3_023_012 * oneOverDet;
2394
2395 m[A22] = Det3_013_013 * oneOverDet;
2396 m[A23] = Det3_013_012 * mn1OverDet;
2397
2398 m[A33] = Det3_012_012 * oneOverDet;
2399
2400 return;
2401 }
2402
2403 void G4ErrorSymMatrix::invertHaywood4(G4int&
2404 {
2405 invert4(ifail); // For the 4x4 case, the m
2406 // the Haywood method.
2407 }
2408