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