Geant4-11
Data Structures | Public Member Functions | Protected Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | Friends
G4ErrorSymMatrix Class Reference

#include <G4ErrorSymMatrix.hh>

Data Structures

class  G4ErrorSymMatrix_row
 
class  G4ErrorSymMatrix_row_const
 

Public Member Functions

G4ErrorSymMatrix apply (G4double(*f)(G4double, G4int, G4int)) const
 
void assign (const G4ErrorMatrix &m2)
 
void assign (const G4ErrorSymMatrix &m2)
 
G4double determinant () const
 
G4double & fast (G4int row, G4int col)
 
const G4double & fast (G4int row, G4int col) const
 
 G4ErrorSymMatrix ()
 
 G4ErrorSymMatrix (const G4ErrorSymMatrix &m1)
 
 G4ErrorSymMatrix (G4ErrorSymMatrix &&)=default
 
 G4ErrorSymMatrix (G4int p)
 
 G4ErrorSymMatrix (G4int p, G4int)
 
G4ErrorSymMatrix inverse (G4int &ifail) const
 
void invert (G4int &ifail)
 
void invertBunchKaufman (G4int &ifail)
 
void invertCholesky5 (G4int &ifail)
 
void invertCholesky6 (G4int &ifail)
 
void invertHaywood4 (G4int &ifail)
 
void invertHaywood5 (G4int &ifail)
 
void invertHaywood6 (G4int &ifail)
 
G4int num_col () const
 
G4int num_row () const
 
G4double & operator() (G4int row, G4int col)
 
const G4double & operator() (G4int row, G4int col) const
 
G4ErrorSymMatrix & operator*= (G4double t)
 
G4ErrorSymMatrix & operator+= (const G4ErrorSymMatrix &m2)
 
G4ErrorSymMatrix operator- () const
 
G4ErrorSymMatrix & operator-= (const G4ErrorSymMatrix &m2)
 
G4ErrorSymMatrix & operator/= (G4double t)
 
G4ErrorSymMatrix & operator= (const G4ErrorSymMatrix &m2)
 
G4ErrorSymMatrix & operator= (G4ErrorSymMatrix &&)=default
 
G4ErrorSymMatrix_row operator[] (G4int)
 
G4ErrorSymMatrix_row_const operator[] (G4int) const
 
G4ErrorSymMatrix similarity (const G4ErrorMatrix &m1) const
 
G4ErrorSymMatrix similarity (const G4ErrorSymMatrix &m1) const
 
G4ErrorSymMatrix similarityT (const G4ErrorMatrix &m1) const
 
G4ErrorSymMatrix sub (G4int min_row, G4int max_row)
 
G4ErrorSymMatrix sub (G4int min_row, G4int max_row) const
 
void sub (G4int row, const G4ErrorSymMatrix &m1)
 
G4ErrorSymMatrix T () const
 
G4double trace () const
 
virtual ~G4ErrorSymMatrix ()
 

Protected Member Functions

G4int num_size () const
 

Private Member Functions

void invert4 (G4int &ifail)
 
void invert5 (G4int &ifail)
 
void invert6 (G4int &ifail)
 

Private Attributes

std::vector< G4double > m
 
G4int nrow
 
G4int size
 

Static Private Attributes

static G4ThreadLocal G4double adjustment5x5 = 0.0
 
static G4ThreadLocal G4double adjustment6x6 = 0.0
 
static const G4double CHOLESKY_CREEP_5x5 = .005
 
static const G4double CHOLESKY_CREEP_6x6 = .002
 
static const G4double CHOLESKY_THRESHOLD_5x5 = .5
 
static const G4double CHOLESKY_THRESHOLD_6x6 = .2
 
static G4ThreadLocal G4double posDefFraction5x5 = 1.0
 
static G4ThreadLocal G4double posDefFraction6x6 = 1.0
 

Friends

G4double condition (const G4ErrorSymMatrix &m)
 
void diag_step (G4ErrorSymMatrix *t, G4ErrorMatrix *u, G4int begin, G4int end)
 
void diag_step (G4ErrorSymMatrix *t, G4int begin, G4int end)
 
G4ErrorMatrix diagonalize (G4ErrorSymMatrix *s)
 
class G4ErrorMatrix
 
class G4ErrorSymMatrix_row
 
class G4ErrorSymMatrix_row_const
 
void house_with_update2 (G4ErrorSymMatrix *a, G4ErrorMatrix *v, G4int row, G4int col)
 
G4ErrorMatrix operator* (const G4ErrorMatrix &m1, const G4ErrorSymMatrix &m2)
 
G4ErrorMatrix operator* (const G4ErrorSymMatrix &m1, const G4ErrorMatrix &m2)
 
G4ErrorMatrix operator* (const G4ErrorSymMatrix &m1, const G4ErrorSymMatrix &m2)
 
G4ErrorSymMatrix operator+ (const G4ErrorSymMatrix &m1, const G4ErrorSymMatrix &m2)
 
G4ErrorSymMatrix operator- (const G4ErrorSymMatrix &m1, const G4ErrorSymMatrix &m2)
 
void tridiagonal (G4ErrorSymMatrix *a, G4ErrorMatrix *hsm)
 

Detailed Description

Definition at line 43 of file G4ErrorSymMatrix.hh.

Constructor & Destructor Documentation

◆ G4ErrorSymMatrix() [1/5]

G4ErrorSymMatrix::G4ErrorSymMatrix ( )
inline

◆ G4ErrorSymMatrix() [2/5]

G4ErrorSymMatrix::G4ErrorSymMatrix ( G4int  p)
explicit

Definition at line 75 of file G4ErrorSymMatrix.cc.

76 : m(p * (p + 1) / 2)
77 , nrow(p)
78{
79 size = nrow * (nrow + 1) / 2;
80 m.assign(size, 0);
81}
G4ErrorSymMatrix::m
std::vector< G4double > m
Definition: G4ErrorSymMatrix.hh:219

References m, nrow, and size.

◆ G4ErrorSymMatrix() [3/5]

G4ErrorSymMatrix::G4ErrorSymMatrix ( G4int  p,
G4int  init 
)

Definition at line 83 of file G4ErrorSymMatrix.cc.

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: initialization must be 0 or 1.");
107 }
108}
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4int
int G4int
Definition: G4Types.hh:85
G4ErrorMatrix::error
static void error(const char *s)
Definition: G4ErrorMatrix.cc:855

References G4ErrorMatrix::error(), m, nrow, and size.

◆ G4ErrorSymMatrix() [4/5]

G4ErrorSymMatrix::G4ErrorSymMatrix ( const G4ErrorSymMatrix &  m1)

Definition at line 116 of file G4ErrorSymMatrix.cc.

117 : m(mat1.size)
118 , nrow(mat1.nrow)
119 , size(mat1.size)
120{
121 m = mat1.m;
122}

References m.

◆ G4ErrorSymMatrix() [5/5]

G4ErrorSymMatrix::G4ErrorSymMatrix ( G4ErrorSymMatrix &&  )
default

◆ ~G4ErrorSymMatrix()

G4ErrorSymMatrix::~G4ErrorSymMatrix ( )
virtual

Definition at line 114 of file G4ErrorSymMatrix.cc.

114{}

Member Function Documentation

◆ apply()

G4ErrorSymMatrix G4ErrorSymMatrix::apply ( G4double(*)(G4double, G4int, G4int)  f) const

Definition at line 592 of file G4ErrorSymMatrix.cc.

594{
595 G4ErrorSymMatrix mret(num_row());
596 G4ErrorMatrixConstIter a = m.begin();
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}
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
G4ErrorSymMatrix::num_row
G4int num_row() const

References m, and num_row().

◆ assign() [1/2]

void G4ErrorSymMatrix::assign ( const G4ErrorMatrix &  m2)

Definition at line 608 of file G4ErrorSymMatrix.cc.

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}

References G4ErrorMatrix::m, m, G4ErrorMatrix::nrow, nrow, and size.

◆ assign() [2/2]

void G4ErrorSymMatrix::assign ( const G4ErrorSymMatrix &  m2)

◆ determinant()

G4double G4ErrorSymMatrix::determinant ( ) const

Definition at line 841 of file G4ErrorSymMatrix.cc.

842{
843 static const G4int max_array = 20;
844
845 // ir must point to an array which is ***1 longer than*** nrow
846
847 static std::vector<G4int> ir_vec(max_array + 1);
848 if(ir_vec.size() <= static_cast<unsigned int>(nrow))
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}
G4double
double G4double
Definition: G4Types.hh:83

References G4ErrorMatrix::dfact_matrix(), and nrow.

◆ fast() [1/2]

G4double & G4ErrorSymMatrix::fast ( G4int  row,
G4int  col 
)

◆ fast() [2/2]

const G4double & G4ErrorSymMatrix::fast ( G4int  row,
G4int  col 
) const

◆ inverse()

G4ErrorSymMatrix G4ErrorSymMatrix::inverse ( G4int &  ifail) const
inline

◆ invert()

void G4ErrorSymMatrix::invert ( G4int &  ifail)

Definition at line 724 of file G4ErrorSymMatrix.cc.

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() + 5)) -
736 (*(m.begin() + 4)) * (*(m.begin() + 4));
737 c12 = (*(m.begin() + 4)) * (*(m.begin() + 3)) -
738 (*(m.begin() + 1)) * (*(m.begin() + 5));
739 c13 = (*(m.begin() + 1)) * (*(m.begin() + 4)) -
740 (*(m.begin() + 2)) * (*(m.begin() + 3));
741 c22 = (*(m.begin() + 5)) * (*m.begin()) -
742 (*(m.begin() + 3)) * (*(m.begin() + 3));
743 c23 = (*(m.begin() + 3)) * (*(m.begin() + 1)) -
744 (*(m.begin() + 4)) * (*m.begin());
745 c33 = (*m.begin()) * (*(m.begin() + 2)) -
746 (*(m.begin() + 1)) * (*(m.begin() + 1));
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() + 1));
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}
G4ErrorSymMatrix::invert6
void invert6(G4int &ifail)
Definition: G4ErrorSymMatrix.cc:1422
G4ErrorSymMatrix::invert4
void invert4(G4int &ifail)
Definition: G4ErrorSymMatrix.cc:2328
G4ErrorSymMatrix::invert5
void invert5(G4int &ifail)
Definition: G4ErrorSymMatrix.cc:1390
G4ErrorSymMatrix::invertBunchKaufman
void invertBunchKaufman(G4int &ifail)
Definition: G4ErrorSymMatrix.cc:874

References invert4(), invert5(), invert6(), invertBunchKaufman(), m, and nrow.

◆ invert4()

void G4ErrorSymMatrix::invert4 ( G4int &  ifail)
private

Definition at line 2328 of file G4ErrorSymMatrix.cc.

2329{
2330 ifail = 0;
2331
2332 // Find all NECESSARY 2x2 dets: (14 of them)
2333
2334 G4double Det2_12_01 = m[A10] * m[A21] - m[A11] * m[A20];
2335 G4double Det2_12_02 = m[A10] * m[A22] - m[A12] * m[A20];
2336 G4double Det2_12_12 = m[A11] * m[A22] - m[A12] * m[A21];
2337 G4double Det2_13_01 = m[A10] * m[A31] - m[A11] * m[A30];
2338 G4double Det2_13_02 = m[A10] * m[A32] - m[A12] * m[A30];
2339 G4double Det2_13_03 = m[A10] * m[A33] - m[A13] * m[A30];
2340 G4double Det2_13_12 = m[A11] * m[A32] - m[A12] * m[A31];
2341 G4double Det2_13_13 = m[A11] * m[A33] - m[A13] * m[A31];
2342 G4double Det2_23_01 = m[A20] * m[A31] - m[A21] * m[A30];
2343 G4double Det2_23_02 = m[A20] * m[A32] - m[A22] * m[A30];
2344 G4double Det2_23_03 = m[A20] * m[A33] - m[A23] * m[A30];
2345 G4double Det2_23_12 = m[A21] * m[A32] - m[A22] * m[A31];
2346 G4double Det2_23_13 = m[A21] * m[A33] - m[A23] * m[A31];
2347 G4double Det2_23_23 = m[A22] * m[A33] - m[A23] * m[A32];
2348
2349 // Find all NECESSARY 3x3 dets: (10 of them)
2350
2351 G4double Det3_012_012 =
2352 m[A00] * Det2_12_12 - m[A01] * Det2_12_02 + m[A02] * Det2_12_01;
2353 G4double Det3_013_012 =
2354 m[A00] * Det2_13_12 - m[A01] * Det2_13_02 + m[A02] * Det2_13_01;
2355 G4double Det3_013_013 =
2356 m[A00] * Det2_13_13 - m[A01] * Det2_13_03 + m[A03] * Det2_13_01;
2357 G4double Det3_023_012 =
2358 m[A00] * Det2_23_12 - m[A01] * Det2_23_02 + m[A02] * Det2_23_01;
2359 G4double Det3_023_013 =
2360 m[A00] * Det2_23_13 - m[A01] * Det2_23_03 + m[A03] * Det2_23_01;
2361 G4double Det3_023_023 =
2362 m[A00] * Det2_23_23 - m[A02] * Det2_23_03 + m[A03] * Det2_23_02;
2363 G4double Det3_123_012 =
2364 m[A10] * Det2_23_12 - m[A11] * Det2_23_02 + m[A12] * Det2_23_01;
2365 G4double Det3_123_013 =
2366 m[A10] * Det2_23_13 - m[A11] * Det2_23_03 + m[A13] * Det2_23_01;
2367 G4double Det3_123_023 =
2368 m[A10] * Det2_23_23 - m[A12] * Det2_23_03 + m[A13] * Det2_23_02;
2369 G4double Det3_123_123 =
2370 m[A11] * Det2_23_23 - m[A12] * Det2_23_13 + m[A13] * Det2_23_12;
2371
2372 // Find the 4x4 det:
2373
2374 G4double det = m[A00] * Det3_123_123 - m[A01] * Det3_123_023 +
2375 m[A02] * Det3_123_013 - m[A03] * Det3_123_012;
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}
A20
#define A20
Definition: G4ErrorSymMatrix.cc:1362
A01
#define A01
Definition: G4ErrorSymMatrix.cc:1349
A23
#define A23
Definition: G4ErrorSymMatrix.cc:1365
A13
#define A13
Definition: G4ErrorSymMatrix.cc:1358
A00
#define A00
Definition: G4ErrorSymMatrix.cc:1348
A02
#define A02
Definition: G4ErrorSymMatrix.cc:1350
A22
#define A22
Definition: G4ErrorSymMatrix.cc:1364
A30
#define A30
Definition: G4ErrorSymMatrix.cc:1369
A12
#define A12
Definition: G4ErrorSymMatrix.cc:1357
A03
#define A03
Definition: G4ErrorSymMatrix.cc:1351
A21
#define A21
Definition: G4ErrorSymMatrix.cc:1363
A11
#define A11
Definition: G4ErrorSymMatrix.cc:1356
A10
#define A10
Definition: G4ErrorSymMatrix.cc:1355
A32
#define A32
Definition: G4ErrorSymMatrix.cc:1371
A31
#define A31
Definition: G4ErrorSymMatrix.cc:1370
A33
#define A33
Definition: G4ErrorSymMatrix.cc:1372

References A00, A01, A02, A03, A10, A11, A12, A13, A20, A21, A22, A23, A30, A31, A32, A33, and m.

Referenced by invert(), and invertHaywood4().

◆ invert5()

void G4ErrorSymMatrix::invert5 ( G4int &  ifail)
private

Definition at line 1390 of file G4ErrorSymMatrix.cc.

1391{
1392 if(posDefFraction5x5 >= CHOLESKY_THRESHOLD_5x5)
1393 {
1394 invertCholesky5(ifail);
1395 posDefFraction5x5 = .9 * posDefFraction5x5 + .1 * (1 - ifail);
1396 if(ifail != 0) // Cholesky failed -- invert using Haywood
1397 {
1398 invertHaywood5(ifail);
1399 }
1400 }
1401 else
1402 {
1403 if(posDefFraction5x5 + adjustment5x5 >= CHOLESKY_THRESHOLD_5x5)
1404 {
1405 invertCholesky5(ifail);
1406 posDefFraction5x5 = .9 * posDefFraction5x5 + .1 * (1 - ifail);
1407 if(ifail != 0) // Cholesky failed -- invert using Haywood
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}
G4ErrorSymMatrix::invertHaywood5
void invertHaywood5(G4int &ifail)
Definition: G4ErrorSymMatrix.cc:1454
G4ErrorSymMatrix::adjustment5x5
static G4ThreadLocal G4double adjustment5x5
Definition: G4ErrorSymMatrix.hh:225
G4ErrorSymMatrix::CHOLESKY_CREEP_5x5
static const G4double CHOLESKY_CREEP_5x5
Definition: G4ErrorSymMatrix.hh:227
G4ErrorSymMatrix::invertCholesky5
void invertCholesky5(G4int &ifail)
Definition: G4ErrorSymMatrix.cc:2023
G4ErrorSymMatrix::posDefFraction5x5
static G4ThreadLocal G4double posDefFraction5x5
Definition: G4ErrorSymMatrix.hh:224
G4ErrorSymMatrix::CHOLESKY_THRESHOLD_5x5
static const G4double CHOLESKY_THRESHOLD_5x5
Definition: G4ErrorSymMatrix.hh:226

References adjustment5x5, CHOLESKY_CREEP_5x5, CHOLESKY_THRESHOLD_5x5, invertCholesky5(), invertHaywood5(), and posDefFraction5x5.

Referenced by invert().

◆ invert6()

void G4ErrorSymMatrix::invert6 ( G4int &  ifail)
private

Definition at line 1422 of file G4ErrorSymMatrix.cc.

1423{
1424 if(posDefFraction6x6 >= CHOLESKY_THRESHOLD_6x6)
1425 {
1426 invertCholesky6(ifail);
1427 posDefFraction6x6 = .9 * posDefFraction6x6 + .1 * (1 - ifail);
1428 if(ifail != 0) // Cholesky failed -- invert using Haywood
1429 {
1430 invertHaywood6(ifail);
1431 }
1432 }
1433 else
1434 {
1435 if(posDefFraction6x6 + adjustment6x6 >= CHOLESKY_THRESHOLD_6x6)
1436 {
1437 invertCholesky6(ifail);
1438 posDefFraction6x6 = .9 * posDefFraction6x6 + .1 * (1 - ifail);
1439 if(ifail != 0) // Cholesky failed -- invert using Haywood
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}
G4ErrorSymMatrix::invertCholesky6
void invertCholesky6(G4int &ifail)
Definition: G4ErrorSymMatrix.cc:2159
G4ErrorSymMatrix::invertHaywood6
void invertHaywood6(G4int &ifail)
Definition: G4ErrorSymMatrix.cc:1620
G4ErrorSymMatrix::posDefFraction6x6
static G4ThreadLocal G4double posDefFraction6x6
Definition: G4ErrorSymMatrix.hh:229
G4ErrorSymMatrix::adjustment6x6
static G4ThreadLocal G4double adjustment6x6
Definition: G4ErrorSymMatrix.hh:230
G4ErrorSymMatrix::CHOLESKY_CREEP_6x6
static const G4double CHOLESKY_CREEP_6x6
Definition: G4ErrorSymMatrix.hh:232
G4ErrorSymMatrix::CHOLESKY_THRESHOLD_6x6
static const G4double CHOLESKY_THRESHOLD_6x6
Definition: G4ErrorSymMatrix.hh:231

References adjustment6x6, CHOLESKY_CREEP_6x6, CHOLESKY_THRESHOLD_6x6, invertCholesky6(), invertHaywood6(), and posDefFraction6x6.

Referenced by invert().

◆ invertBunchKaufman()

void G4ErrorSymMatrix::invertBunchKaufman ( G4int &  ifail)

Definition at line 874 of file G4ErrorSymMatrix.cc.

875{
876 // Bunch-Kaufman diagonal pivoting method
877 // It is decribed in J.R. Bunch, L. Kaufman (1977).
878 // "Some Stable Methods for Calculating Inertia and Solving Symmetric
879 // Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub,
880 // Charles F. van Loan, "Matrix Computations" (the second edition
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 needed: x and piv are
888 // used as pointers to arrays of doubles and ints respectively, each
889 // of length nrow. We do not want to reallocate each time through
890 // unless the size needs to grow. We do not want to leak memory, even
891 // by having a new without a delete that is only done once.
892
893 static const G4int max_array = 25;
894 static G4ThreadLocal std::vector<G4double>* xvec = 0;
895 if(!xvec)
896 xvec = new std::vector<G4double>(max_array);
897 static G4ThreadLocal std::vector<G4int>* pivv = 0;
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>(nrow))
902 xvec->resize(nrow);
903 if(pivv->size() < static_cast<unsigned int>(nrow))
904 pivv->resize(nrow);
905 // Note - resize should do nothing if the size is already larger than nrow,
906 // but on VC++ there are indications that it does so we check.
907 // Note - the data elements in a vector are guaranteed to be contiguous,
908 // so x[i] and piv[i] are optimally fast.
909 G4ErrorMatrixIter x = xvec->begin();
910 // x[i] is used as helper storage, needs to have at least size nrow.
911 pivIter piv = pivv->begin();
912 // piv[i] is used to store details of exchanges
913
914 G4double temp1, temp2;
915 G4ErrorMatrixIter ip, mjj, iq;
916 G4double lambda, sigma;
917 const G4double alpha = .6404; // = (1+sqrt(17))/8
918 const G4double epsilon = 32 * DBL_EPSILON;
919 // whenever a sum of two doubles is below or equal to epsilon
920 // it is set to zero.
921 // this constant could be set to zero but then the algorithm
922 // doesn't neccessarily detect that a matrix is singular
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 * L^T
932 // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
933 // L and D^-1 are stored in A = *this, P is stored in piv[]
934
935 for(j = 1; j < nrow; j += ss) // main loop over columns
936 {
937 mjj = m.begin() + j * (j - 1) / 2 + j - 1;
938 lambda = 0; // compute lambda = max of A(j+1:n,j)
939 pivrow = j + 1;
940 ip = m.begin() + (j + 1) * j / 2 + j - 1;
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(pivrow, j:pivrow-1)
969 ip = m.begin() + pivrow * (pivrow - 1) / 2 + j - 1;
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 * lambda * lambda)
977 {
978 ss = 1;
979 pivrow = j;
980 }
981 else if(std::fabs(*(m.begin() + pivrow * (pivrow - 1) / 2 + pivrow -
982 1)) >= alpha * sigma)
983 {
984 ss = 1;
985 }
986 else
987 {
988 ss = 2;
989 }
990 }
991 if(pivrow == j) // no permutation neccessary
992 {
993 piv[j - 1] = pivrow;
994 if(*mjj == 0)
995 {
996 ifail = 1;
997 return;
998 }
999 temp2 = *mjj = 1. / *mjj; // invert D(j,j)
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) / 2 + j - 1) * temp2;
1005 ip = m.begin() + i * (i - 1) / 2 + j;
1006 for(k = j + 1; k <= i; k++)
1007 {
1008 *ip -= temp1 * *(m.begin() + k * (k - 1) / 2 + j - 1);
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 - 1;
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 pivrow in
1028 // submatrix (j:n,j:n)
1029 ip = m.begin() + pivrow * (pivrow - 1) / 2 + j;
1030 for(i = j + 1; i < pivrow; i++, ip++)
1031 {
1032 temp1 = *(m.begin() + i * (i - 1) / 2 + j - 1);
1033 *(m.begin() + i * (i - 1) / 2 + j - 1) = *ip;
1034 *ip = temp1;
1035 }
1036 temp1 = *mjj;
1037 *mjj = *(m.begin() + pivrow * (pivrow - 1) / 2 + pivrow - 1);
1038 *(m.begin() + pivrow * (pivrow - 1) / 2 + pivrow - 1) = temp1;
1039 ip = m.begin() + (pivrow + 1) * pivrow / 2 + j - 1;
1040 iq = ip + pivrow - j;
1041 for(i = pivrow + 1; i <= nrow; ip += i, iq += i++)
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 D(j,j)
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) / 2 + j - 1) * temp2;
1059 ip = m.begin() + i * (i - 1) / 2 + j;
1060 for(k = j + 1; k <= i; k++)
1061 {
1062 *ip -= temp1 * *(m.begin() + k * (k - 1) / 2 + j - 1);
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 - 1;
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 is the second row of a 2x2 pivot
1081
1082 if(j + 1 != pivrow)
1083 {
1084 // interchange rows and columns j+1 and pivrow in
1085 // submatrix (j:n,j:n)
1086 ip = m.begin() + pivrow * (pivrow - 1) / 2 + j + 1;
1087 for(i = j + 2; i < pivrow; i++, ip++)
1088 {
1089 temp1 = *(m.begin() + i * (i - 1) / 2 + j);
1090 *(m.begin() + i * (i - 1) / 2 + j) = *ip;
1091 *ip = temp1;
1092 }
1093 temp1 = *(mjj + j + 1);
1094 *(mjj + j + 1) =
1095 *(m.begin() + pivrow * (pivrow - 1) / 2 + pivrow - 1);
1096 *(m.begin() + pivrow * (pivrow - 1) / 2 + pivrow - 1) = temp1;
1097 temp1 = *(mjj + j);
1098 *(mjj + j) = *(m.begin() + pivrow * (pivrow - 1) / 2 + j - 1);
1099 *(m.begin() + pivrow * (pivrow - 1) / 2 + j - 1) = temp1;
1100 ip = m.begin() + (pivrow + 1) * pivrow / 2 + j;
1101 iq = ip + pivrow - (j + 1);
1102 for(i = pivrow + 1; i <= nrow; ip += i, iq += i++)
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 + j) * *(mjj + j);
1111 if(temp2 == 0)
1112 {
1113 G4Exception("G4ErrorSymMatrix::bunch_invert()",
1114 "GEANT4e-Notification", JustWarning,
1115 "Error in pivot choice!");
1116 }
1117 temp2 = 1. / temp2;
1118
1119 // this quotient is guaranteed to exist by the choice
1120 // of the pivot
1121
1122 temp1 = *mjj;
1123 *mjj = *(mjj + j + 1) * temp2;
1124 *(mjj + j + 1) = temp1 * temp2;
1125 *(mjj + j) = -*(mjj + j) * temp2;
1126
1127 if(j < nrow - 1) // otherwise do nothing
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) / 2 + j - 1;
1133 temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
1134 if(std::fabs(temp1) <= epsilon)
1135 {
1136 temp1 = 0;
1137 }
1138 temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
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) / 2 + k - 1;
1146 iq = m.begin() + k * (k - 1) / 2 + j - 1;
1147 *ip -= temp1 * *iq + temp2 * *(iq + 1);
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) / 2 + j - 1;
1158 temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
1159 if(std::fabs(temp1) <= epsilon)
1160 {
1161 temp1 = 0;
1162 }
1163 *(ip + 1) = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
1164 if(std::fabs(*(ip + 1)) <= epsilon)
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 factorization
1190
1191 for(j = nrow; j >= 1; j -= ss) // loop over columns
1192 {
1193 mjj = m.begin() + j * (j - 1) / 2 + j - 1;
1194 if(piv[j - 1] > 0) // 1x1 pivot, compute column j of inverse
1195 {
1196 ss = 1;
1197 if(j < nrow)
1198 {
1199 ip = m.begin() + (j + 1) * j / 2 + j - 1;
1200 for(i = 0; i < nrow - j; ip += 1 + j + i++)
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 + j;
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 += 1 + j + k++)
1213 {
1214 temp2 += *ip * x[k];
1215 }
1216 *(m.begin() + i * (i - 1) / 2 + j - 1) = -temp2;
1217 }
1218 temp2 = 0;
1219 ip = m.begin() + (j + 1) * j / 2 + j - 1;
1220 for(k = 0; k < nrow - j; ip += 1 + j + k++)
1221 {
1222 temp2 += x[k] * *ip;
1223 }
1224 *mjj -= temp2;
1225 }
1226 }
1227 else // 2x2 pivot, compute columns j and j-1 of the inverse
1228 {
1229 if(piv[j - 1] != 0)
1230 {
1231 std::ostringstream message;
1232 message << "Error in pivot: " << piv[j - 1];
1233 G4Exception("G4ErrorSymMatrix::invertBunchKaufman()",
1234 "GEANT4e-Notification", JustWarning, message);
1235 }
1236 ss = 2;
1237 if(j < nrow)
1238 {
1239 ip = m.begin() + (j + 1) * j / 2 + j - 1;
1240 for(i = 0; i < nrow - j; ip += 1 + j + i++)
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 + j;
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 += 1 + j + k++)
1253 {
1254 temp2 += *ip * x[k];
1255 }
1256 *(m.begin() + i * (i - 1) / 2 + j - 1) = -temp2;
1257 }
1258 temp2 = 0;
1259 ip = m.begin() + (j + 1) * j / 2 + j - 1;
1260 for(k = 0; k < nrow - j; ip += 1 + j + k++)
1261 {
1262 temp2 += x[k] * *ip;
1263 }
1264 *mjj -= temp2;
1265 temp2 = 0;
1266 ip = m.begin() + (j + 1) * j / 2 + 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 - 2;
1273 for(i = 0; i < nrow - j; ip += 1 + j + i++)
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 + j;
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 += 1 + j + k++)
1286 {
1287 temp2 += *ip * x[k];
1288 }
1289 *(m.begin() + i * (i - 1) / 2 + j - 2) = -temp2;
1290 }
1291 temp2 = 0;
1292 ip = m.begin() + (j + 1) * j / 2 + j - 2;
1293 for(k = 0; k < nrow - j; ip += 1 + j + k++)
1294 {
1295 temp2 += x[k] * *ip;
1296 }
1297 *(mjj - j) -= temp2;
1298 }
1299 }
1300
1301 // interchange rows and columns j and piv[j-1]
1302 // or rows and columns j and -piv[j-2]
1303
1304 pivrow = (piv[j - 1] == 0) ? -piv[j - 2] : piv[j - 1];
1305 ip = m.begin() + pivrow * (pivrow - 1) / 2 + j;
1306 for(i = j + 1; i < pivrow; i++, ip++)
1307 {
1308 temp1 = *(m.begin() + i * (i - 1) / 2 + j - 1);
1309 *(m.begin() + i * (i - 1) / 2 + j - 1) = *ip;
1310 *ip = temp1;
1311 }
1312 temp1 = *mjj;
1313 *mjj = *(m.begin() + pivrow * (pivrow - 1) / 2 + pivrow - 1);
1314 *(m.begin() + pivrow * (pivrow - 1) / 2 + pivrow - 1) = temp1;
1315 if(ss == 2)
1316 {
1317 temp1 = *(mjj - 1);
1318 *(mjj - 1) = *(m.begin() + pivrow * (pivrow - 1) / 2 + j - 2);
1319 *(m.begin() + pivrow * (pivrow - 1) / 2 + j - 2) = temp1;
1320 }
1321
1322 ip = m.begin() + (pivrow + 1) * pivrow / 2 + j - 1; // &A(i,j)
1323 iq = ip + pivrow - j;
1324 for(i = pivrow + 1; i <= nrow; ip += i, iq += i++)
1325 {
1326 temp1 = *iq;
1327 *iq = *ip;
1328 *ip = temp1;
1329 }
1330 } // end of loop over columns (in computing inverse from factorization)
1331
1332 return; // inversion successful
1333}
epsilon
G4double epsilon(G4double density, G4double temperature)
Definition: G4DNAElectronHoleRecombination.cc:95
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:65
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
alpha
static const G4double alpha
Definition: G4LivermoreBremsstrahlungModel.cc:80
DBL_EPSILON
#define DBL_EPSILON
Definition: templates.hh:66
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:77

References alpha, DBL_EPSILON, epsilon(), G4Exception(), G4ThreadLocal, JustWarning, G4InuclParticleNames::lambda, m, and nrow.

Referenced by invert().

◆ invertCholesky5()

void G4ErrorSymMatrix::invertCholesky5 ( G4int &  ifail)

Definition at line 2023 of file G4ErrorSymMatrix.cc.

2024{
2025 // Invert by
2026 //
2027 // a) decomposing M = G*G^T with G lower triangular
2028 // (if M is not positive definite this will fail, leaving this unchanged)
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 inverted and 1 is returned.
2033 // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2034
2035 G4double h10; // below-diagonal elements of H
2036 G4double h20, h21;
2037 G4double h30, h31, h32;
2038 G4double h40, h41, h42, h43;
2039
2040 G4double h00, h11, h22, h33, h44; // 1/diagonal elements of G =
2041 // diagonal elements of H
2042
2043 G4double g10; // below-diagonal elements of G
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 succeed...
2049
2050 // Form G -- compute diagonal members of H directly rather than of G
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 products from rest of column 1 and
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 products from rest of column 2 and
2093 // scale to get column 2 of G
2094
2095 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2096 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2097
2098 // Form G33 (actually, just h33)
2099
2100 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2101 if(h33 <= 0)
2102 {
2103 return;
2104 }
2105 h33 = 1.0 / std::sqrt(h33);
2106
2107 // Subtract inter-column column dot product from A43 and scale to get G43
2108
2109 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2110
2111 // Finally form h44 - if this is possible inversion succeeds
2112
2113 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
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 are already correct
2121 //-------------
2122
2123 // The order here is dictated by speed considerations
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 * h44);
2131 h10 = -h00 * g10 * h11;
2132 h20 = -h00 * (g10 * h21 + g20 * h22);
2133 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2134 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2135
2136 // Change this to its inverse = H^T*H
2137 //------------------------------------
2138
2139 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40;
2140 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41;
2141 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41;
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}
A41
#define A41
Definition: G4ErrorSymMatrix.cc:1377
A40
#define A40
Definition: G4ErrorSymMatrix.cc:1376
A24
#define A24
Definition: G4ErrorSymMatrix.cc:1366
A04
#define A04
Definition: G4ErrorSymMatrix.cc:1352
A14
#define A14
Definition: G4ErrorSymMatrix.cc:1359
A44
#define A44
Definition: G4ErrorSymMatrix.cc:1380
A42
#define A42
Definition: G4ErrorSymMatrix.cc:1378
A34
#define A34
Definition: G4ErrorSymMatrix.cc:1373
A43
#define A43
Definition: G4ErrorSymMatrix.cc:1379

References A00, A01, A02, A03, A04, A10, A11, A12, A13, A14, A20, A21, A22, A23, A24, A30, A31, A32, A33, A34, A40, A41, A42, A43, A44, and m.

Referenced by invert5().

◆ invertCholesky6()

void G4ErrorSymMatrix::invertCholesky6 ( G4int &  ifail)

Definition at line 2159 of file G4ErrorSymMatrix.cc.

2160{
2161 // Invert by
2162 //
2163 // a) decomposing M = G*G^T with G lower triangular
2164 // (if M is not positive definite this will fail, leaving this unchanged)
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 inverted and 1 is returned.
2169 // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2170
2171 G4double h10; // below-diagonal elements of H
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; // 1/diagonal elements of G =
2178 // diagonal elements of H
2179
2180 G4double g10; // below-diagonal elements of G
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 succeed...
2187
2188 // Form G -- compute diagonal members of H directly rather than of G
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 products from rest of column 1 and
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 products from rest of column 2 and
2233 // scale to get column 2 of G
2234
2235 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2236 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2237 g52 = (m[A52] - (g20 * g50) - (g21 * g51)) * h22;
2238
2239 // Form G33 (actually, just h33)
2240
2241 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2242 if(h33 <= 0)
2243 {
2244 return;
2245 }
2246 h33 = 1.0 / std::sqrt(h33);
2247
2248 // Subtract inter-column column dot products from rest of column 3 and
2249 // scale to get column 3 of G
2250
2251 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2252 g53 = (m[A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33;
2253
2254 // Form G44 (actually, just h44)
2255
2256 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
2257 if(h44 <= 0)
2258 {
2259 return;
2260 }
2261 h44 = 1.0 / std::sqrt(h44);
2262
2263 // Subtract inter-column column dot product from M54 and scale to get G54
2264
2265 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;
2266
2267 // Finally form h55 - if this is possible inversion succeeds
2268
2269 h55 = m[A55] - (g50 * g50) - (g51 * g51) - (g52 * g52) - (g53 * g53) -
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 are already correct
2278 //-------------
2279
2280 // The order here is dictated by speed considerations
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 * h55);
2288 h21 = -h11 * g21 * h22;
2289 h31 = -h11 * (g21 * h32 + g31 * h33);
2290 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2291 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55);
2292 h10 = -h00 * g10 * h11;
2293 h20 = -h00 * (g10 * h21 + g20 * h22);
2294 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2295 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2296 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55);
2297
2298 // Change this to its inverse = H^T*H
2299 //------------------------------------
2300
2301 m[A00] =
2302 h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40 + h50 * h50;
2303 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51;
2304 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51;
2305 m[A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52;
2306 m[A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52;
2307 m[A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52;
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}
A53
#define A53
Definition: G4ErrorSymMatrix.cc:1386
A45
#define A45
Definition: G4ErrorSymMatrix.cc:1381
A54
#define A54
Definition: G4ErrorSymMatrix.cc:1387
A25
#define A25
Definition: G4ErrorSymMatrix.cc:1367
A55
#define A55
Definition: G4ErrorSymMatrix.cc:1388
A35
#define A35
Definition: G4ErrorSymMatrix.cc:1374
A50
#define A50
Definition: G4ErrorSymMatrix.cc:1383
A51
#define A51
Definition: G4ErrorSymMatrix.cc:1384
A05
#define A05
Definition: G4ErrorSymMatrix.cc:1353
A52
#define A52
Definition: G4ErrorSymMatrix.cc:1385
A15
#define A15
Definition: G4ErrorSymMatrix.cc:1360

References A00, A01, A02, A03, A04, A05, A10, A11, A12, A13, A14, A15, A20, A21, A22, A23, A24, A25, A30, A31, A32, A33, A34, A35, A40, A41, A42, A43, A44, A45, A50, A51, A52, A53, A54, A55, and m.

Referenced by invert6().

◆ invertHaywood4()

void G4ErrorSymMatrix::invertHaywood4 ( G4int &  ifail)

Definition at line 2403 of file G4ErrorSymMatrix.cc.

2404{
2405 invert4(ifail); // For the 4x4 case, the method we use for invert is already
2406 // the Haywood method.
2407}

References invert4().

◆ invertHaywood5()

void G4ErrorSymMatrix::invertHaywood5 ( G4int &  ifail)

Definition at line 1454 of file G4ErrorSymMatrix.cc.

1455{
1456 ifail = 0;
1457
1458 // Find all NECESSARY 2x2 dets: (25 of them)
1459
1460 G4double Det2_23_01 = m[A20] * m[A31] - m[A21] * m[A30];
1461 G4double Det2_23_02 = m[A20] * m[A32] - m[A22] * m[A30];
1462 G4double Det2_23_03 = m[A20] * m[A33] - m[A23] * m[A30];
1463 G4double Det2_23_12 = m[A21] * m[A32] - m[A22] * m[A31];
1464 G4double Det2_23_13 = m[A21] * m[A33] - m[A23] * m[A31];
1465 G4double Det2_23_23 = m[A22] * m[A33] - m[A23] * m[A32];
1466 G4double Det2_24_01 = m[A20] * m[A41] - m[A21] * m[A40];
1467 G4double Det2_24_02 = m[A20] * m[A42] - m[A22] * m[A40];
1468 G4double Det2_24_03 = m[A20] * m[A43] - m[A23] * m[A40];
1469 G4double Det2_24_04 = m[A20] * m[A44] - m[A24] * m[A40];
1470 G4double Det2_24_12 = m[A21] * m[A42] - m[A22] * m[A41];
1471 G4double Det2_24_13 = m[A21] * m[A43] - m[A23] * m[A41];
1472 G4double Det2_24_14 = m[A21] * m[A44] - m[A24] * m[A41];
1473 G4double Det2_24_23 = m[A22] * m[A43] - m[A23] * m[A42];
1474 G4double Det2_24_24 = m[A22] * m[A44] - m[A24] * m[A42];
1475 G4double Det2_34_01 = m[A30] * m[A41] - m[A31] * m[A40];
1476 G4double Det2_34_02 = m[A30] * m[A42] - m[A32] * m[A40];
1477 G4double Det2_34_03 = m[A30] * m[A43] - m[A33] * m[A40];
1478 G4double Det2_34_04 = m[A30] * m[A44] - m[A34] * m[A40];
1479 G4double Det2_34_12 = m[A31] * m[A42] - m[A32] * m[A41];
1480 G4double Det2_34_13 = m[A31] * m[A43] - m[A33] * m[A41];
1481 G4double Det2_34_14 = m[A31] * m[A44] - m[A34] * m[A41];
1482 G4double Det2_34_23 = m[A32] * m[A43] - m[A33] * m[A42];
1483 G4double Det2_34_24 = m[A32] * m[A44] - m[A34] * m[A42];
1484 G4double Det2_34_34 = m[A33] * m[A44] - m[A34] * m[A43];
1485
1486 // Find all NECESSARY 3x3 dets: (30 of them)
1487
1488 G4double Det3_123_012 =
1489 m[A10] * Det2_23_12 - m[A11] * Det2_23_02 + m[A12] * Det2_23_01;
1490 G4double Det3_123_013 =
1491 m[A10] * Det2_23_13 - m[A11] * Det2_23_03 + m[A13] * Det2_23_01;
1492 G4double Det3_123_023 =
1493 m[A10] * Det2_23_23 - m[A12] * Det2_23_03 + m[A13] * Det2_23_02;
1494 G4double Det3_123_123 =
1495 m[A11] * Det2_23_23 - m[A12] * Det2_23_13 + m[A13] * Det2_23_12;
1496 G4double Det3_124_012 =
1497 m[A10] * Det2_24_12 - m[A11] * Det2_24_02 + m[A12] * Det2_24_01;
1498 G4double Det3_124_013 =
1499 m[A10] * Det2_24_13 - m[A11] * Det2_24_03 + m[A13] * Det2_24_01;
1500 G4double Det3_124_014 =
1501 m[A10] * Det2_24_14 - m[A11] * Det2_24_04 + m[A14] * Det2_24_01;
1502 G4double Det3_124_023 =
1503 m[A10] * Det2_24_23 - m[A12] * Det2_24_03 + m[A13] * Det2_24_02;
1504 G4double Det3_124_024 =
1505 m[A10] * Det2_24_24 - m[A12] * Det2_24_04 + m[A14] * Det2_24_02;
1506 G4double Det3_124_123 =
1507 m[A11] * Det2_24_23 - m[A12] * Det2_24_13 + m[A13] * Det2_24_12;
1508 G4double Det3_124_124 =
1509 m[A11] * Det2_24_24 - m[A12] * Det2_24_14 + m[A14] * Det2_24_12;
1510 G4double Det3_134_012 =
1511 m[A10] * Det2_34_12 - m[A11] * Det2_34_02 + m[A12] * Det2_34_01;
1512 G4double Det3_134_013 =
1513 m[A10] * Det2_34_13 - m[A11] * Det2_34_03 + m[A13] * Det2_34_01;
1514 G4double Det3_134_014 =
1515 m[A10] * Det2_34_14 - m[A11] * Det2_34_04 + m[A14] * Det2_34_01;
1516 G4double Det3_134_023 =
1517 m[A10] * Det2_34_23 - m[A12] * Det2_34_03 + m[A13] * Det2_34_02;
1518 G4double Det3_134_024 =
1519 m[A10] * Det2_34_24 - m[A12] * Det2_34_04 + m[A14] * Det2_34_02;
1520 G4double Det3_134_034 =
1521 m[A10] * Det2_34_34 - m[A13] * Det2_34_04 + m[A14] * Det2_34_03;
1522 G4double Det3_134_123 =
1523 m[A11] * Det2_34_23 - m[A12] * Det2_34_13 + m[A13] * Det2_34_12;
1524 G4double Det3_134_124 =
1525 m[A11] * Det2_34_24 - m[A12] * Det2_34_14 + m[A14] * Det2_34_12;
1526 G4double Det3_134_134 =
1527 m[A11] * Det2_34_34 - m[A13] * Det2_34_14 + m[A14] * Det2_34_13;
1528 G4double Det3_234_012 =
1529 m[A20] * Det2_34_12 - m[A21] * Det2_34_02 + m[A22] * Det2_34_01;
1530 G4double Det3_234_013 =
1531 m[A20] * Det2_34_13 - m[A21] * Det2_34_03 + m[A23] * Det2_34_01;
1532 G4double Det3_234_014 =
1533 m[A20] * Det2_34_14 - m[A21] * Det2_34_04 + m[A24] * Det2_34_01;
1534 G4double Det3_234_023 =
1535 m[A20] * Det2_34_23 - m[A22] * Det2_34_03 + m[A23] * Det2_34_02;
1536 G4double Det3_234_024 =
1537 m[A20] * Det2_34_24 - m[A22] * Det2_34_04 + m[A24] * Det2_34_02;
1538 G4double Det3_234_034 =
1539 m[A20] * Det2_34_34 - m[A23] * Det2_34_04 + m[A24] * Det2_34_03;
1540 G4double Det3_234_123 =
1541 m[A21] * Det2_34_23 - m[A22] * Det2_34_13 + m[A23] * Det2_34_12;
1542 G4double Det3_234_124 =
1543 m[A21] * Det2_34_24 - m[A22] * Det2_34_14 + m[A24] * Det2_34_12;
1544 G4double Det3_234_134 =
1545 m[A21] * Det2_34_34 - m[A23] * Det2_34_14 + m[A24] * Det2_34_13;
1546 G4double Det3_234_234 =
1547 m[A22] * Det2_34_34 - m[A23] * Det2_34_24 + m[A24] * Det2_34_23;
1548
1549 // Find all NECESSARY 4x4 dets: (15 of them)
1550
1551 G4double Det4_0123_0123 = m[A00] * Det3_123_123 - m[A01] * Det3_123_023 +
1552 m[A02] * Det3_123_013 - m[A03] * Det3_123_012;
1553 G4double Det4_0124_0123 = m[A00] * Det3_124_123 - m[A01] * Det3_124_023 +
1554 m[A02] * Det3_124_013 - m[A03] * Det3_124_012;
1555 G4double Det4_0124_0124 = m[A00] * Det3_124_124 - m[A01] * Det3_124_024 +
1556 m[A02] * Det3_124_014 - m[A04] * Det3_124_012;
1557 G4double Det4_0134_0123 = m[A00] * Det3_134_123 - m[A01] * Det3_134_023 +
1558 m[A02] * Det3_134_013 - m[A03] * Det3_134_012;
1559 G4double Det4_0134_0124 = m[A00] * Det3_134_124 - m[A01] * Det3_134_024 +
1560 m[A02] * Det3_134_014 - m[A04] * Det3_134_012;
1561 G4double Det4_0134_0134 = m[A00] * Det3_134_134 - m[A01] * Det3_134_034 +
1562 m[A03] * Det3_134_014 - m[A04] * Det3_134_013;
1563 G4double Det4_0234_0123 = m[A00] * Det3_234_123 - m[A01] * Det3_234_023 +
1564 m[A02] * Det3_234_013 - m[A03] * Det3_234_012;
1565 G4double Det4_0234_0124 = m[A00] * Det3_234_124 - m[A01] * Det3_234_024 +
1566 m[A02] * Det3_234_014 - m[A04] * Det3_234_012;
1567 G4double Det4_0234_0134 = m[A00] * Det3_234_134 - m[A01] * Det3_234_034 +
1568 m[A03] * Det3_234_014 - m[A04] * Det3_234_013;
1569 G4double Det4_0234_0234 = m[A00] * Det3_234_234 - m[A02] * Det3_234_034 +
1570 m[A03] * Det3_234_024 - m[A04] * Det3_234_023;
1571 G4double Det4_1234_0123 = m[A10] * Det3_234_123 - m[A11] * Det3_234_023 +
1572 m[A12] * Det3_234_013 - m[A13] * Det3_234_012;
1573 G4double Det4_1234_0124 = m[A10] * Det3_234_124 - m[A11] * Det3_234_024 +
1574 m[A12] * Det3_234_014 - m[A14] * Det3_234_012;
1575 G4double Det4_1234_0134 = m[A10] * Det3_234_134 - m[A11] * Det3_234_034 +
1576 m[A13] * Det3_234_014 - m[A14] * Det3_234_013;
1577 G4double Det4_1234_0234 = m[A10] * Det3_234_234 - m[A12] * Det3_234_034 +
1578 m[A13] * Det3_234_024 - m[A14] * Det3_234_023;
1579 G4double Det4_1234_1234 = m[A11] * Det3_234_234 - m[A12] * Det3_234_134 +
1580 m[A13] * Det3_234_124 - m[A14] * Det3_234_123;
1581
1582 // Find the 5x5 det:
1583
1584 G4double det = m[A00] * Det4_1234_1234 - m[A01] * Det4_1234_0234 +
1585 m[A02] * Det4_1234_0134 - m[A03] * Det4_1234_0124 +
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}

References A00, A01, A02, A03, A04, A10, A11, A12, A13, A14, A20, A21, A22, A23, A24, A30, A31, A32, A33, A34, A40, A41, A42, A43, A44, and m.

Referenced by invert5().

◆ invertHaywood6()

void G4ErrorSymMatrix::invertHaywood6 ( G4int &  ifail)

Definition at line 1620 of file G4ErrorSymMatrix.cc.

1621{
1622 ifail = 0;
1623
1624 // Find all NECESSARY 2x2 dets: (39 of them)
1625
1626 G4double Det2_34_01 = m[A30] * m[A41] - m[A31] * m[A40];
1627 G4double Det2_34_02 = m[A30] * m[A42] - m[A32] * m[A40];
1628 G4double Det2_34_03 = m[A30] * m[A43] - m[A33] * m[A40];
1629 G4double Det2_34_04 = m[A30] * m[A44] - m[A34] * m[A40];
1630 G4double Det2_34_12 = m[A31] * m[A42] - m[A32] * m[A41];
1631 G4double Det2_34_13 = m[A31] * m[A43] - m[A33] * m[A41];
1632 G4double Det2_34_14 = m[A31] * m[A44] - m[A34] * m[A41];
1633 G4double Det2_34_23 = m[A32] * m[A43] - m[A33] * m[A42];
1634 G4double Det2_34_24 = m[A32] * m[A44] - m[A34] * m[A42];
1635 G4double Det2_34_34 = m[A33] * m[A44] - m[A34] * m[A43];
1636 G4double Det2_35_01 = m[A30] * m[A51] - m[A31] * m[A50];
1637 G4double Det2_35_02 = m[A30] * m[A52] - m[A32] * m[A50];
1638 G4double Det2_35_03 = m[A30] * m[A53] - m[A33] * m[A50];
1639 G4double Det2_35_04 = m[A30] * m[A54] - m[A34] * m[A50];
1640 G4double Det2_35_05 = m[A30] * m[A55] - m[A35] * m[A50];
1641 G4double Det2_35_12 = m[A31] * m[A52] - m[A32] * m[A51];
1642 G4double Det2_35_13 = m[A31] * m[A53] - m[A33] * m[A51];
1643 G4double Det2_35_14 = m[A31] * m[A54] - m[A34] * m[A51];
1644 G4double Det2_35_15 = m[A31] * m[A55] - m[A35] * m[A51];
1645 G4double Det2_35_23 = m[A32] * m[A53] - m[A33] * m[A52];
1646 G4double Det2_35_24 = m[A32] * m[A54] - m[A34] * m[A52];
1647 G4double Det2_35_25 = m[A32] * m[A55] - m[A35] * m[A52];
1648 G4double Det2_35_34 = m[A33] * m[A54] - m[A34] * m[A53];
1649 G4double Det2_35_35 = m[A33] * m[A55] - m[A35] * m[A53];
1650 G4double Det2_45_01 = m[A40] * m[A51] - m[A41] * m[A50];
1651 G4double Det2_45_02 = m[A40] * m[A52] - m[A42] * m[A50];
1652 G4double Det2_45_03 = m[A40] * m[A53] - m[A43] * m[A50];
1653 G4double Det2_45_04 = m[A40] * m[A54] - m[A44] * m[A50];
1654 G4double Det2_45_05 = m[A40] * m[A55] - m[A45] * m[A50];
1655 G4double Det2_45_12 = m[A41] * m[A52] - m[A42] * m[A51];
1656 G4double Det2_45_13 = m[A41] * m[A53] - m[A43] * m[A51];
1657 G4double Det2_45_14 = m[A41] * m[A54] - m[A44] * m[A51];
1658 G4double Det2_45_15 = m[A41] * m[A55] - m[A45] * m[A51];
1659 G4double Det2_45_23 = m[A42] * m[A53] - m[A43] * m[A52];
1660 G4double Det2_45_24 = m[A42] * m[A54] - m[A44] * m[A52];
1661 G4double Det2_45_25 = m[A42] * m[A55] - m[A45] * m[A52];
1662 G4double Det2_45_34 = m[A43] * m[A54] - m[A44] * m[A53];
1663 G4double Det2_45_35 = m[A43] * m[A55] - m[A45] * m[A53];
1664 G4double Det2_45_45 = m[A44] * m[A55] - m[A45] * m[A54];
1665
1666 // Find all NECESSARY 3x3 dets: (65 of them)
1667
1668 G4double Det3_234_012 =
1669 m[A20] * Det2_34_12 - m[A21] * Det2_34_02 + m[A22] * Det2_34_01;
1670 G4double Det3_234_013 =
1671 m[A20] * Det2_34_13 - m[A21] * Det2_34_03 + m[A23] * Det2_34_01;
1672 G4double Det3_234_014 =
1673 m[A20] * Det2_34_14 - m[A21] * Det2_34_04 + m[A24] * Det2_34_01;
1674 G4double Det3_234_023 =
1675 m[A20] * Det2_34_23 - m[A22] * Det2_34_03 + m[A23] * Det2_34_02;
1676 G4double Det3_234_024 =
1677 m[A20] * Det2_34_24 - m[A22] * Det2_34_04 + m[A24] * Det2_34_02;
1678 G4double Det3_234_034 =
1679 m[A20] * Det2_34_34 - m[A23] * Det2_34_04 + m[A24] * Det2_34_03;
1680 G4double Det3_234_123 =
1681 m[A21] * Det2_34_23 - m[A22] * Det2_34_13 + m[A23] * Det2_34_12;
1682 G4double Det3_234_124 =
1683 m[A21] * Det2_34_24 - m[A22] * Det2_34_14 + m[A24] * Det2_34_12;
1684 G4double Det3_234_134 =
1685 m[A21] * Det2_34_34 - m[A23] * Det2_34_14 + m[A24] * Det2_34_13;
1686 G4double Det3_234_234 =
1687 m[A22] * Det2_34_34 - m[A23] * Det2_34_24 + m[A24] * Det2_34_23;
1688 G4double Det3_235_012 =
1689 m[A20] * Det2_35_12 - m[A21] * Det2_35_02 + m[A22] * Det2_35_01;
1690 G4double Det3_235_013 =
1691 m[A20] * Det2_35_13 - m[A21] * Det2_35_03 + m[A23] * Det2_35_01;
1692 G4double Det3_235_014 =
1693 m[A20] * Det2_35_14 - m[A21] * Det2_35_04 + m[A24] * Det2_35_01;
1694 G4double Det3_235_015 =
1695 m[A20] * Det2_35_15 - m[A21] * Det2_35_05 + m[A25] * Det2_35_01;
1696 G4double Det3_235_023 =
1697 m[A20] * Det2_35_23 - m[A22] * Det2_35_03 + m[A23] * Det2_35_02;
1698 G4double Det3_235_024 =
1699 m[A20] * Det2_35_24 - m[A22] * Det2_35_04 + m[A24] * Det2_35_02;
1700 G4double Det3_235_025 =
1701 m[A20] * Det2_35_25 - m[A22] * Det2_35_05 + m[A25] * Det2_35_02;
1702 G4double Det3_235_034 =
1703 m[A20] * Det2_35_34 - m[A23] * Det2_35_04 + m[A24] * Det2_35_03;
1704 G4double Det3_235_035 =
1705 m[A20] * Det2_35_35 - m[A23] * Det2_35_05 + m[A25] * Det2_35_03;
1706 G4double Det3_235_123 =
1707 m[A21] * Det2_35_23 - m[A22] * Det2_35_13 + m[A23] * Det2_35_12;
1708 G4double Det3_235_124 =
1709 m[A21] * Det2_35_24 - m[A22] * Det2_35_14 + m[A24] * Det2_35_12;
1710 G4double Det3_235_125 =
1711 m[A21] * Det2_35_25 - m[A22] * Det2_35_15 + m[A25] * Det2_35_12;
1712 G4double Det3_235_134 =
1713 m[A21] * Det2_35_34 - m[A23] * Det2_35_14 + m[A24] * Det2_35_13;
1714 G4double Det3_235_135 =
1715 m[A21] * Det2_35_35 - m[A23] * Det2_35_15 + m[A25] * Det2_35_13;
1716 G4double Det3_235_234 =
1717 m[A22] * Det2_35_34 - m[A23] * Det2_35_24 + m[A24] * Det2_35_23;
1718 G4double Det3_235_235 =
1719 m[A22] * Det2_35_35 - m[A23] * Det2_35_25 + m[A25] * Det2_35_23;
1720 G4double Det3_245_012 =
1721 m[A20] * Det2_45_12 - m[A21] * Det2_45_02 + m[A22] * Det2_45_01;
1722 G4double Det3_245_013 =
1723 m[A20] * Det2_45_13 - m[A21] * Det2_45_03 + m[A23] * Det2_45_01;
1724 G4double Det3_245_014 =
1725 m[A20] * Det2_45_14 - m[A21] * Det2_45_04 + m[A24] * Det2_45_01;
1726 G4double Det3_245_015 =
1727 m[A20] * Det2_45_15 - m[A21] * Det2_45_05 + m[A25] * Det2_45_01;
1728 G4double Det3_245_023 =
1729 m[A20] * Det2_45_23 - m[A22] * Det2_45_03 + m[A23] * Det2_45_02;
1730 G4double Det3_245_024 =
1731 m[A20] * Det2_45_24 - m[A22] * Det2_45_04 + m[A24] * Det2_45_02;
1732 G4double Det3_245_025 =
1733 m[A20] * Det2_45_25 - m[A22] * Det2_45_05 + m[A25] * Det2_45_02;
1734 G4double Det3_245_034 =
1735 m[A20] * Det2_45_34 - m[A23] * Det2_45_04 + m[A24] * Det2_45_03;
1736 G4double Det3_245_035 =
1737 m[A20] * Det2_45_35 - m[A23] * Det2_45_05 + m[A25] * Det2_45_03;
1738 G4double Det3_245_045 =
1739 m[A20] * Det2_45_45 - m[A24] * Det2_45_05 + m[A25] * Det2_45_04;
1740 G4double Det3_245_123 =
1741 m[A21] * Det2_45_23 - m[A22] * Det2_45_13 + m[A23] * Det2_45_12;
1742 G4double Det3_245_124 =
1743 m[A21] * Det2_45_24 - m[A22] * Det2_45_14 + m[A24] * Det2_45_12;
1744 G4double Det3_245_125 =
1745 m[A21] * Det2_45_25 - m[A22] * Det2_45_15 + m[A25] * Det2_45_12;
1746 G4double Det3_245_134 =
1747 m[A21] * Det2_45_34 - m[A23] * Det2_45_14 + m[A24] * Det2_45_13;
1748 G4double Det3_245_135 =
1749 m[A21] * Det2_45_35 - m[A23] * Det2_45_15 + m[A25] * Det2_45_13;
1750 G4double Det3_245_145 =
1751 m[A21] * Det2_45_45 - m[A24] * Det2_45_15 + m[A25] * Det2_45_14;
1752 G4double Det3_245_234 =
1753 m[A22] * Det2_45_34 - m[A23] * Det2_45_24 + m[A24] * Det2_45_23;
1754 G4double Det3_245_235 =
1755 m[A22] * Det2_45_35 - m[A23] * Det2_45_25 + m[A25] * Det2_45_23;
1756 G4double Det3_245_245 =
1757 m[A22] * Det2_45_45 - m[A24] * Det2_45_25 + m[A25] * Det2_45_24;
1758 G4double Det3_345_012 =
1759 m[A30] * Det2_45_12 - m[A31] * Det2_45_02 + m[A32] * Det2_45_01;
1760 G4double Det3_345_013 =
1761 m[A30] * Det2_45_13 - m[A31] * Det2_45_03 + m[A33] * Det2_45_01;
1762 G4double Det3_345_014 =
1763 m[A30] * Det2_45_14 - m[A31] * Det2_45_04 + m[A34] * Det2_45_01;
1764 G4double Det3_345_015 =
1765 m[A30] * Det2_45_15 - m[A31] * Det2_45_05 + m[A35] * Det2_45_01;
1766 G4double Det3_345_023 =
1767 m[A30] * Det2_45_23 - m[A32] * Det2_45_03 + m[A33] * Det2_45_02;
1768 G4double Det3_345_024 =
1769 m[A30] * Det2_45_24 - m[A32] * Det2_45_04 + m[A34] * Det2_45_02;
1770 G4double Det3_345_025 =
1771 m[A30] * Det2_45_25 - m[A32] * Det2_45_05 + m[A35] * Det2_45_02;
1772 G4double Det3_345_034 =
1773 m[A30] * Det2_45_34 - m[A33] * Det2_45_04 + m[A34] * Det2_45_03;
1774 G4double Det3_345_035 =
1775 m[A30] * Det2_45_35 - m[A33] * Det2_45_05 + m[A35] * Det2_45_03;
1776 G4double Det3_345_045 =
1777 m[A30] * Det2_45_45 - m[A34] * Det2_45_05 + m[A35] * Det2_45_04;
1778 G4double Det3_345_123 =
1779 m[A31] * Det2_45_23 - m[A32] * Det2_45_13 + m[A33] * Det2_45_12;
1780 G4double Det3_345_124 =
1781 m[A31] * Det2_45_24 - m[A32] * Det2_45_14 + m[A34] * Det2_45_12;
1782 G4double Det3_345_125 =
1783 m[A31] * Det2_45_25 - m[A32] * Det2_45_15 + m[A35] * Det2_45_12;
1784 G4double Det3_345_134 =
1785 m[A31] * Det2_45_34 - m[A33] * Det2_45_14 + m[A34] * Det2_45_13;
1786 G4double Det3_345_135 =
1787 m[A31] * Det2_45_35 - m[A33] * Det2_45_15 + m[A35] * Det2_45_13;
1788 G4double Det3_345_145 =
1789 m[A31] * Det2_45_45 - m[A34] * Det2_45_15 + m[A35] * Det2_45_14;
1790 G4double Det3_345_234 =
1791 m[A32] * Det2_45_34 - m[A33] * Det2_45_24 + m[A34] * Det2_45_23;
1792 G4double Det3_345_235 =
1793 m[A32] * Det2_45_35 - m[A33] * Det2_45_25 + m[A35] * Det2_45_23;
1794 G4double Det3_345_245 =
1795 m[A32] * Det2_45_45 - m[A34] * Det2_45_25 + m[A35] * Det2_45_24;
1796 G4double Det3_345_345 =
1797 m[A33] * Det2_45_45 - m[A34] * Det2_45_35 + m[A35] * Det2_45_34;
1798
1799 // Find all NECESSARY 4x4 dets: (55 of them)
1800
1801 G4double Det4_1234_0123 = m[A10] * Det3_234_123 - m[A11] * Det3_234_023 +
1802 m[A12] * Det3_234_013 - m[A13] * Det3_234_012;
1803 G4double Det4_1234_0124 = m[A10] * Det3_234_124 - m[A11] * Det3_234_024 +
1804 m[A12] * Det3_234_014 - m[A14] * Det3_234_012;
1805 G4double Det4_1234_0134 = m[A10] * Det3_234_134 - m[A11] * Det3_234_034 +
1806 m[A13] * Det3_234_014 - m[A14] * Det3_234_013;
1807 G4double Det4_1234_0234 = m[A10] * Det3_234_234 - m[A12] * Det3_234_034 +
1808 m[A13] * Det3_234_024 - m[A14] * Det3_234_023;
1809 G4double Det4_1234_1234 = m[A11] * Det3_234_234 - m[A12] * Det3_234_134 +
1810 m[A13] * Det3_234_124 - m[A14] * Det3_234_123;
1811 G4double Det4_1235_0123 = m[A10] * Det3_235_123 - m[A11] * Det3_235_023 +
1812 m[A12] * Det3_235_013 - m[A13] * Det3_235_012;
1813 G4double Det4_1235_0124 = m[A10] * Det3_235_124 - m[A11] * Det3_235_024 +
1814 m[A12] * Det3_235_014 - m[A14] * Det3_235_012;
1815 G4double Det4_1235_0125 = m[A10] * Det3_235_125 - m[A11] * Det3_235_025 +
1816 m[A12] * Det3_235_015 - m[A15] * Det3_235_012;
1817 G4double Det4_1235_0134 = m[A10] * Det3_235_134 - m[A11] * Det3_235_034 +
1818 m[A13] * Det3_235_014 - m[A14] * Det3_235_013;
1819 G4double Det4_1235_0135 = m[A10] * Det3_235_135 - m[A11] * Det3_235_035 +
1820 m[A13] * Det3_235_015 - m[A15] * Det3_235_013;
1821 G4double Det4_1235_0234 = m[A10] * Det3_235_234 - m[A12] * Det3_235_034 +
1822 m[A13] * Det3_235_024 - m[A14] * Det3_235_023;
1823 G4double Det4_1235_0235 = m[A10] * Det3_235_235 - m[A12] * Det3_235_035 +
1824 m[A13] * Det3_235_025 - m[A15] * Det3_235_023;
1825 G4double Det4_1235_1234 = m[A11] * Det3_235_234 - m[A12] * Det3_235_134 +
1826 m[A13] * Det3_235_124 - m[A14] * Det3_235_123;
1827 G4double Det4_1235_1235 = m[A11] * Det3_235_235 - m[A12] * Det3_235_135 +
1828 m[A13] * Det3_235_125 - m[A15] * Det3_235_123;
1829 G4double Det4_1245_0123 = m[A10] * Det3_245_123 - m[A11] * Det3_245_023 +
1830 m[A12] * Det3_245_013 - m[A13] * Det3_245_012;
1831 G4double Det4_1245_0124 = m[A10] * Det3_245_124 - m[A11] * Det3_245_024 +
1832 m[A12] * Det3_245_014 - m[A14] * Det3_245_012;
1833 G4double Det4_1245_0125 = m[A10] * Det3_245_125 - m[A11] * Det3_245_025 +
1834 m[A12] * Det3_245_015 - m[A15] * Det3_245_012;
1835 G4double Det4_1245_0134 = m[A10] * Det3_245_134 - m[A11] * Det3_245_034 +
1836 m[A13] * Det3_245_014 - m[A14] * Det3_245_013;
1837 G4double Det4_1245_0135 = m[A10] * Det3_245_135 - m[A11] * Det3_245_035 +
1838 m[A13] * Det3_245_015 - m[A15] * Det3_245_013;
1839 G4double Det4_1245_0145 = m[A10] * Det3_245_145 - m[A11] * Det3_245_045 +
1840 m[A14] * Det3_245_015 - m[A15] * Det3_245_014;
1841 G4double Det4_1245_0234 = m[A10] * Det3_245_234 - m[A12] * Det3_245_034 +
1842 m[A13] * Det3_245_024 - m[A14] * Det3_245_023;
1843 G4double Det4_1245_0235 = m[A10] * Det3_245_235 - m[A12] * Det3_245_035 +
1844 m[A13] * Det3_245_025 - m[A15] * Det3_245_023;
1845 G4double Det4_1245_0245 = m[A10] * Det3_245_245 - m[A12] * Det3_245_045 +
1846 m[A14] * Det3_245_025 - m[A15] * Det3_245_024;
1847 G4double Det4_1245_1234 = m[A11] * Det3_245_234 - m[A12] * Det3_245_134 +
1848 m[A13] * Det3_245_124 - m[A14] * Det3_245_123;
1849 G4double Det4_1245_1235 = m[A11] * Det3_245_235 - m[A12] * Det3_245_135 +
1850 m[A13] * Det3_245_125 - m[A15] * Det3_245_123;
1851 G4double Det4_1245_1245 = m[A11] * Det3_245_245 - m[A12] * Det3_245_145 +
1852 m[A14] * Det3_245_125 - m[A15] * Det3_245_124;
1853 G4double Det4_1345_0123 = m[A10] * Det3_345_123 - m[A11] * Det3_345_023 +
1854 m[A12] * Det3_345_013 - m[A13] * Det3_345_012;
1855 G4double Det4_1345_0124 = m[A10] * Det3_345_124 - m[A11] * Det3_345_024 +
1856 m[A12] * Det3_345_014 - m[A14] * Det3_345_012;
1857 G4double Det4_1345_0125 = m[A10] * Det3_345_125 - m[A11] * Det3_345_025 +
1858 m[A12] * Det3_345_015 - m[A15] * Det3_345_012;
1859 G4double Det4_1345_0134 = m[A10] * Det3_345_134 - m[A11] * Det3_345_034 +
1860 m[A13] * Det3_345_014 - m[A14] * Det3_345_013;
1861 G4double Det4_1345_0135 = m[A10] * Det3_345_135 - m[A11] * Det3_345_035 +
1862 m[A13] * Det3_345_015 - m[A15] * Det3_345_013;
1863 G4double Det4_1345_0145 = m[A10] * Det3_345_145 - m[A11] * Det3_345_045 +
1864 m[A14] * Det3_345_015 - m[A15] * Det3_345_014;
1865 G4double Det4_1345_0234 = m[A10] * Det3_345_234 - m[A12] * Det3_345_034 +
1866 m[A13] * Det3_345_024 - m[A14] * Det3_345_023;
1867 G4double Det4_1345_0235 = m[A10] * Det3_345_235 - m[A12] * Det3_345_035 +
1868 m[A13] * Det3_345_025 - m[A15] * Det3_345_023;
1869 G4double Det4_1345_0245 = m[A10] * Det3_345_245 - m[A12] * Det3_345_045 +
1870 m[A14] * Det3_345_025 - m[A15] * Det3_345_024;
1871 G4double Det4_1345_0345 = m[A10] * Det3_345_345 - m[A13] * Det3_345_045 +
1872 m[A14] * Det3_345_035 - m[A15] * Det3_345_034;
1873 G4double Det4_1345_1234 = m[A11] * Det3_345_234 - m[A12] * Det3_345_134 +
1874 m[A13] * Det3_345_124 - m[A14] * Det3_345_123;
1875 G4double Det4_1345_1235 = m[A11] * Det3_345_235 - m[A12] * Det3_345_135 +
1876 m[A13] * Det3_345_125 - m[A15] * Det3_345_123;
1877 G4double Det4_1345_1245 = m[A11] * Det3_345_245 - m[A12] * Det3_345_145 +
1878 m[A14] * Det3_345_125 - m[A15] * Det3_345_124;
1879 G4double Det4_1345_1345 = m[A11] * Det3_345_345 - m[A13] * Det3_345_145 +
1880 m[A14] * Det3_345_135 - m[A15] * Det3_345_134;
1881 G4double Det4_2345_0123 = m[A20] * Det3_345_123 - m[A21] * Det3_345_023 +
1882 m[A22] * Det3_345_013 - m[A23] * Det3_345_012;
1883 G4double Det4_2345_0124 = m[A20] * Det3_345_124 - m[A21] * Det3_345_024 +
1884 m[A22] * Det3_345_014 - m[A24] * Det3_345_012;
1885 G4double Det4_2345_0125 = m[A20] * Det3_345_125 - m[A21] * Det3_345_025 +
1886 m[A22] * Det3_345_015 - m[A25] * Det3_345_012;
1887 G4double Det4_2345_0134 = m[A20] * Det3_345_134 - m[A21] * Det3_345_034 +
1888 m[A23] * Det3_345_014 - m[A24] * Det3_345_013;
1889 G4double Det4_2345_0135 = m[A20] * Det3_345_135 - m[A21] * Det3_345_035 +
1890 m[A23] * Det3_345_015 - m[A25] * Det3_345_013;
1891 G4double Det4_2345_0145 = m[A20] * Det3_345_145 - m[A21] * Det3_345_045 +
1892 m[A24] * Det3_345_015 - m[A25] * Det3_345_014;
1893 G4double Det4_2345_0234 = m[A20] * Det3_345_234 - m[A22] * Det3_345_034 +
1894 m[A23] * Det3_345_024 - m[A24] * Det3_345_023;
1895 G4double Det4_2345_0235 = m[A20] * Det3_345_235 - m[A22] * Det3_345_035 +
1896 m[A23] * Det3_345_025 - m[A25] * Det3_345_023;
1897 G4double Det4_2345_0245 = m[A20] * Det3_345_245 - m[A22] * Det3_345_045 +
1898 m[A24] * Det3_345_025 - m[A25] * Det3_345_024;
1899 G4double Det4_2345_0345 = m[A20] * Det3_345_345 - m[A23] * Det3_345_045 +
1900 m[A24] * Det3_345_035 - m[A25] * Det3_345_034;
1901 G4double Det4_2345_1234 = m[A21] * Det3_345_234 - m[A22] * Det3_345_134 +
1902 m[A23] * Det3_345_124 - m[A24] * Det3_345_123;
1903 G4double Det4_2345_1235 = m[A21] * Det3_345_235 - m[A22] * Det3_345_135 +
1904 m[A23] * Det3_345_125 - m[A25] * Det3_345_123;
1905 G4double Det4_2345_1245 = m[A21] * Det3_345_245 - m[A22] * Det3_345_145 +
1906 m[A24] * Det3_345_125 - m[A25] * Det3_345_124;
1907 G4double Det4_2345_1345 = m[A21] * Det3_345_345 - m[A23] * Det3_345_145 +
1908 m[A24] * Det3_345_135 - m[A25] * Det3_345_134;
1909 G4double Det4_2345_2345 = m[A22] * Det3_345_345 - m[A23] * Det3_345_245 +
1910 m[A24] * Det3_345_235 - m[A25] * Det3_345_234;
1911
1912 // Find all NECESSARY 5x5 dets: (19 of them)
1913
1914 G4double Det5_01234_01234 =
1915 m[A00] * Det4_1234_1234 - m[A01] * Det4_1234_0234 +
1916 m[A02] * Det4_1234_0134 - m[A03] * Det4_1234_0124 + m[A04] * Det4_1234_0123;
1917 G4double Det5_01235_01234 =
1918 m[A00] * Det4_1235_1234 - m[A01] * Det4_1235_0234 +
1919 m[A02] * Det4_1235_0134 - m[A03] * Det4_1235_0124 + m[A04] * Det4_1235_0123;
1920 G4double Det5_01235_01235 =
1921 m[A00] * Det4_1235_1235 - m[A01] * Det4_1235_0235 +
1922 m[A02] * Det4_1235_0135 - m[A03] * Det4_1235_0125 + m[A05] * Det4_1235_0123;
1923 G4double Det5_01245_01234 =
1924 m[A00] * Det4_1245_1234 - m[A01] * Det4_1245_0234 +
1925 m[A02] * Det4_1245_0134 - m[A03] * Det4_1245_0124 + m[A04] * Det4_1245_0123;
1926 G4double Det5_01245_01235 =
1927 m[A00] * Det4_1245_1235 - m[A01] * Det4_1245_0235 +
1928 m[A02] * Det4_1245_0135 - m[A03] * Det4_1245_0125 + m[A05] * Det4_1245_0123;
1929 G4double Det5_01245_01245 =
1930 m[A00] * Det4_1245_1245 - m[A01] * Det4_1245_0245 +
1931 m[A02] * Det4_1245_0145 - m[A04] * Det4_1245_0125 + m[A05] * Det4_1245_0124;
1932 G4double Det5_01345_01234 =
1933 m[A00] * Det4_1345_1234 - m[A01] * Det4_1345_0234 +
1934 m[A02] * Det4_1345_0134 - m[A03] * Det4_1345_0124 + m[A04] * Det4_1345_0123;
1935 G4double Det5_01345_01235 =
1936 m[A00] * Det4_1345_1235 - m[A01] * Det4_1345_0235 +
1937 m[A02] * Det4_1345_0135 - m[A03] * Det4_1345_0125 + m[A05] * Det4_1345_0123;
1938 G4double Det5_01345_01245 =
1939 m[A00] * Det4_1345_1245 - m[A01] * Det4_1345_0245 +
1940 m[A02] * Det4_1345_0145 - m[A04] * Det4_1345_0125 + m[A05] * Det4_1345_0124;
1941 G4double Det5_01345_01345 =
1942 m[A00] * Det4_1345_1345 - m[A01] * Det4_1345_0345 +
1943 m[A03] * Det4_1345_0145 - m[A04] * Det4_1345_0135 + m[A05] * Det4_1345_0134;
1944 G4double Det5_02345_01234 =
1945 m[A00] * Det4_2345_1234 - m[A01] * Det4_2345_0234 +
1946 m[A02] * Det4_2345_0134 - m[A03] * Det4_2345_0124 + m[A04] * Det4_2345_0123;
1947 G4double Det5_02345_01235 =
1948 m[A00] * Det4_2345_1235 - m[A01] * Det4_2345_0235 +
1949 m[A02] * Det4_2345_0135 - m[A03] * Det4_2345_0125 + m[A05] * Det4_2345_0123;
1950 G4double Det5_02345_01245 =
1951 m[A00] * Det4_2345_1245 - m[A01] * Det4_2345_0245 +
1952 m[A02] * Det4_2345_0145 - m[A04] * Det4_2345_0125 + m[A05] * Det4_2345_0124;
1953 G4double Det5_02345_01345 =
1954 m[A00] * Det4_2345_1345 - m[A01] * Det4_2345_0345 +
1955 m[A03] * Det4_2345_0145 - m[A04] * Det4_2345_0135 + m[A05] * Det4_2345_0134;
1956 G4double Det5_02345_02345 =
1957 m[A00] * Det4_2345_2345 - m[A02] * Det4_2345_0345 +
1958 m[A03] * Det4_2345_0245 - m[A04] * Det4_2345_0235 + m[A05] * Det4_2345_0234;
1959 G4double Det5_12345_01234 =
1960 m[A10] * Det4_2345_1234 - m[A11] * Det4_2345_0234 +
1961 m[A12] * Det4_2345_0134 - m[A13] * Det4_2345_0124 + m[A14] * Det4_2345_0123;
1962 G4double Det5_12345_01235 =
1963 m[A10] * Det4_2345_1235 - m[A11] * Det4_2345_0235 +
1964 m[A12] * Det4_2345_0135 - m[A13] * Det4_2345_0125 + m[A15] * Det4_2345_0123;
1965 G4double Det5_12345_01245 =
1966 m[A10] * Det4_2345_1245 - m[A11] * Det4_2345_0245 +
1967 m[A12] * Det4_2345_0145 - m[A14] * Det4_2345_0125 + m[A15] * Det4_2345_0124;
1968 G4double Det5_12345_01345 =
1969 m[A10] * Det4_2345_1345 - m[A11] * Det4_2345_0345 +
1970 m[A13] * Det4_2345_0145 - m[A14] * Det4_2345_0135 + m[A15] * Det4_2345_0134;
1971 G4double Det5_12345_02345 =
1972 m[A10] * Det4_2345_2345 - m[A12] * Det4_2345_0345 +
1973 m[A13] * Det4_2345_0245 - m[A14] * Det4_2345_0235 + m[A15] * Det4_2345_0234;
1974 G4double Det5_12345_12345 =
1975 m[A11] * Det4_2345_2345 - m[A12] * Det4_2345_1345 +
1976 m[A13] * Det4_2345_1245 - m[A14] * Det4_2345_1235 + m[A15] * Det4_2345_1234;
1977
1978 // Find the determinant
1979
1980 G4double det = m[A00] * Det5_12345_12345 - m[A01] * Det5_12345_02345 +
1981 m[A02] * Det5_12345_01345 - m[A03] * Det5_12345_01245 +
1982 m[A04] * Det5_12345_01235 - m[A05] * Det5_12345_01234;
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}

References A00, A01, A02, A03, A04, A05, A10, A11, A12, A13, A14, A15, A20, A21, A22, A23, A24, A25, A30, A31, A32, A33, A34, A35, A40, A41, A42, A43, A44, A45, A50, A51, A52, A53, A54, A55, and m.

Referenced by invert6().

◆ num_col()

G4int G4ErrorSymMatrix::num_col ( ) const
inline

Referenced by operator+(), G4ErrorMatrix::operator+=(), operator+=(), operator-(), G4ErrorMatrix::operator-=(), operator-=(), operator<<(), and similarity().

◆ num_row()

G4int G4ErrorSymMatrix::num_row ( ) const
inline

Referenced by apply(), dsum(), operator+(), G4ErrorMatrix::operator+=(), operator+=(), operator-(), G4ErrorMatrix::operator-=(), operator-=(), operator<<(), similarity(), similarityT(), and sub().

◆ num_size()

G4int G4ErrorSymMatrix::num_size ( ) const
inlineprotected

Referenced by operator-().

◆ operator()() [1/2]

G4double & G4ErrorSymMatrix::operator() ( G4int  row,
G4int  col 
)

◆ operator()() [2/2]

const G4double & G4ErrorSymMatrix::operator() ( G4int  row,
G4int  col 
) const

◆ operator*=()

G4ErrorSymMatrix & G4ErrorSymMatrix::operator*= ( G4double  t)

Definition at line 508 of file G4ErrorSymMatrix.cc.

509{
510 SIMPLE_UOP(*=)
511 return (*this);
512}
SIMPLE_UOP
#define SIMPLE_UOP(OPER)
Definition: G4ErrorSymMatrix.cc:40

References SIMPLE_UOP.

◆ operator+=()

G4ErrorSymMatrix & G4ErrorSymMatrix::operator+= ( const G4ErrorSymMatrix &  m2)

Definition at line 462 of file G4ErrorSymMatrix.cc.

463{
464 CHK_DIM_2(num_row(), mat2.num_row(), num_col(), mat2.num_col(), +=);
465 SIMPLE_BOP(+=)
466 return (*this);
467}
CHK_DIM_2
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition: G4ErrorSymMatrix.cc:61
SIMPLE_BOP
#define SIMPLE_BOP(OPER)
Definition: G4ErrorSymMatrix.cc:46
G4ErrorSymMatrix::num_col
G4int num_col() const

References CHK_DIM_2, num_col(), num_row(), and SIMPLE_BOP.

◆ operator-()

G4ErrorSymMatrix G4ErrorSymMatrix::operator- ( ) const

Definition at line 209 of file G4ErrorSymMatrix.cc.

210{
211 G4ErrorSymMatrix mat2(nrow);
212 G4ErrorMatrixConstIter a = m.begin();
213 G4ErrorMatrixIter b = mat2.m.begin();
214 G4ErrorMatrixConstIter e = m.begin() + num_size();
215 for(; a < e; a++, b++)
216 {
217 (*b) = -(*a);
218 }
219 return mat2;
220}
G4ErrorSymMatrix::num_size
G4int num_size() const

References m, nrow, and num_size().

◆ operator-=()

G4ErrorSymMatrix & G4ErrorSymMatrix::operator-= ( const G4ErrorSymMatrix &  m2)

Definition at line 495 of file G4ErrorSymMatrix.cc.

496{
497 CHK_DIM_2(num_row(), mat2.num_row(), num_col(), mat2.num_col(), -=);
498 SIMPLE_BOP(-=)
499 return (*this);
500}

References CHK_DIM_2, num_col(), num_row(), and SIMPLE_BOP.

◆ operator/=()

G4ErrorSymMatrix & G4ErrorSymMatrix::operator/= ( G4double  t)

Definition at line 502 of file G4ErrorSymMatrix.cc.

503{
504 SIMPLE_UOP(/=)
505 return (*this);
506}

References SIMPLE_UOP.

◆ operator=() [1/2]

G4ErrorSymMatrix & G4ErrorSymMatrix::operator= ( const G4ErrorSymMatrix &  m2)

Definition at line 546 of file G4ErrorSymMatrix.cc.

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}

References m, nrow, and size.

◆ operator=() [2/2]

G4ErrorSymMatrix & G4ErrorSymMatrix::operator= ( G4ErrorSymMatrix &&  )
default

◆ operator[]() [1/2]

G4ErrorSymMatrix_row G4ErrorSymMatrix::operator[] ( G4int  )
inline

◆ operator[]() [2/2]

G4ErrorSymMatrix_row_const G4ErrorSymMatrix::operator[] ( G4int  ) const
inline

◆ similarity() [1/2]

G4ErrorSymMatrix G4ErrorSymMatrix::similarity ( const G4ErrorMatrix &  m1) const

Definition at line 629 of file G4ErrorSymMatrix.cc.

630{
631 G4ErrorSymMatrix mret(mat1.num_row());
632 G4ErrorMatrix temp = mat1 * (*this);
633
634 // If mat1*(*this) has correct dimensions, then so will the mat1.T
635 // multiplication. So there is no need to check dimensions again.
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}
G4ErrorMatrix::m
std::vector< G4double > m
Definition: G4ErrorMatrix.hh:222

References G4ErrorMatrix::m, m, CLHEP::detail::n, G4ErrorMatrix::num_col(), G4ErrorMatrix::num_row(), and num_row().

Referenced by G4ErrorSurfaceTrajState::BuildErrorMatrix(), G4ErrorFreeTrajState::G4ErrorFreeTrajState(), and G4ErrorFreeTrajState::PropagateError().

◆ similarity() [2/2]

G4ErrorSymMatrix G4ErrorSymMatrix::similarity ( const G4ErrorSymMatrix &  m1) const

Definition at line 660 of file G4ErrorSymMatrix.cc.

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}

References G4ErrorMatrix::m, m, CLHEP::detail::n, num_col(), and num_row().

◆ similarityT()

G4ErrorSymMatrix G4ErrorSymMatrix::similarityT ( const G4ErrorMatrix &  m1) const

Definition at line 695 of file G4ErrorSymMatrix.cc.

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}

References G4ErrorMatrix::m, m, CLHEP::detail::n, G4ErrorMatrix::num_col(), G4ErrorMatrix::num_row(), and num_row().

◆ sub() [1/3]

G4ErrorSymMatrix G4ErrorSymMatrix::sub ( G4int  min_row,
G4int  max_row 
)

Definition at line 151 of file G4ErrorSymMatrix.cc.

152{
153 G4ErrorSymMatrix mret(max_row - min_row + 1);
154 if(max_row > num_row())
155 {
156 G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range");
157 }
158 G4ErrorMatrixIter a = mret.m.begin();
159 G4ErrorMatrixIter b1 = m.begin() + (min_row + 2) * (min_row - 1) / 2;
160 for(G4int irow = 1; irow <= mret.num_row(); irow++)
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}

References G4ErrorMatrix::error(), m, and num_row().

◆ sub() [2/3]

G4ErrorSymMatrix G4ErrorSymMatrix::sub ( G4int  min_row,
G4int  max_row 
) const

Definition at line 130 of file G4ErrorSymMatrix.cc.

131{
132 G4ErrorSymMatrix mret(max_row - min_row + 1);
133 if(max_row > num_row())
134 {
135 G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range");
136 }
137 G4ErrorMatrixIter a = mret.m.begin();
138 G4ErrorMatrixConstIter b1 = m.begin() + (min_row + 2) * (min_row - 1) / 2;
139 for(G4int irow = 1; irow <= mret.num_row(); irow++)
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}

References G4ErrorMatrix::error(), m, and num_row().

Referenced by dsum().

◆ sub() [3/3]

void G4ErrorSymMatrix::sub ( G4int  row,
const G4ErrorSymMatrix &  m1 
)

Definition at line 172 of file G4ErrorSymMatrix.cc.

173{
174 if(row < 1 || row + mat1.num_row() - 1 > num_row())
175 {
176 G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range");
177 }
178 G4ErrorMatrixConstIter a = mat1.m.begin();
179 G4ErrorMatrixIter b1 = m.begin() + (row + 2) * (row - 1) / 2;
180 for(G4int irow = 1; irow <= mat1.num_row(); irow++)
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}

References G4ErrorMatrix::error(), m, and num_row().

◆ T()

G4ErrorSymMatrix G4ErrorSymMatrix::T ( ) const

Referenced by G4ErrorFreeTrajState::PropagateError().

◆ trace()

G4double G4ErrorSymMatrix::trace ( ) const

Definition at line 864 of file G4ErrorSymMatrix.cc.

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}

References m, and nrow.

Friends And Related Function Documentation

◆ condition

G4double condition ( const G4ErrorSymMatrix &  m)
friend

◆ diag_step [1/2]

void diag_step ( G4ErrorSymMatrix *  t,
G4ErrorMatrix *  u,
G4int  begin,
G4int  end 
)
friend

◆ diag_step [2/2]

void diag_step ( G4ErrorSymMatrix *  t,
G4int  begin,
G4int  end 
)
friend

◆ diagonalize

G4ErrorMatrix diagonalize ( G4ErrorSymMatrix *  s)
friend

◆ G4ErrorMatrix

friend class G4ErrorMatrix
friend

Definition at line 194 of file G4ErrorSymMatrix.hh.

◆ G4ErrorSymMatrix_row

friend class G4ErrorSymMatrix_row
friend

Definition at line 192 of file G4ErrorSymMatrix.hh.

◆ G4ErrorSymMatrix_row_const

friend class G4ErrorSymMatrix_row_const
friend

Definition at line 193 of file G4ErrorSymMatrix.hh.

◆ house_with_update2

void house_with_update2 ( G4ErrorSymMatrix *  a,
G4ErrorMatrix *  v,
G4int  row,
G4int  col 
)
friend

◆ operator* [1/3]

G4ErrorMatrix operator* ( const G4ErrorMatrix &  m1,
const G4ErrorSymMatrix &  m2 
)
friend

Definition at line 302 of file G4ErrorSymMatrix.cc.

303{
304 G4ErrorMatrix mret(mat1.num_row(), mat2.num_col());
305 CHK_DIM_1(mat1.num_col(), mat2.num_row(), *);
306 G4ErrorMatrixConstIter mit1, mit2, sp, snp; // mit2=0
307 G4double temp;
308 G4ErrorMatrixIter mir = mret.m.begin();
309 for(mit1 = mat1.m.begin();
310 mit1 < mat1.m.begin() + mat1.num_row() * mat1.num_col(); mit1 = mit2)
311 {
312 snp = mat2.m.begin();
313 for(int step = 1; step <= mat2.num_row(); ++step)
314 {
315 mit2 = mit1;
316 sp = snp;
317 snp += step;
318 temp = 0;
319 while(sp < snp) // Loop checking, 06.08.2015, G.Cosmo
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 <= mat2.num_row(); stept++)
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}
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: G4ErrorSymMatrix.cc:67

◆ operator* [2/3]

G4ErrorMatrix operator* ( const G4ErrorSymMatrix &  m1,
const G4ErrorMatrix &  m2 
)
friend

Definition at line 339 of file G4ErrorSymMatrix.cc.

340{
341 G4ErrorMatrix mret(mat1.num_row(), mat2.num_col());
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 <= mat1.num_row(); snp += step++)
348 {
349 for(mit1 = mat2.m.begin(); mit1 < mat2.m.begin() + mat2.num_col(); mit1++)
350 {
351 mit2 = mit1;
352 sp = snp;
353 temp = 0;
354 while(sp < snp + step) // Loop checking, 06.08.2015, G.Cosmo
355 {
356 temp += *mit2 * (*(sp++));
357 mit2 += mat2.num_col();
358 }
359 sp += step - 1;
360 for(stept = step + 1; stept <= mat1.num_row(); stept++)
361 {
362 temp += *mit2 * (*sp);
363 mit2 += mat2.num_col();
364 sp += stept;
365 }
366 *(mir++) = temp;
367 }
368 }
369 return mret;
370}

◆ operator* [3/3]

G4ErrorMatrix operator* ( const G4ErrorSymMatrix &  m1,
const G4ErrorSymMatrix &  m2 
)
friend

Definition at line 372 of file G4ErrorSymMatrix.cc.

374{
375 G4ErrorMatrix mret(mat1.num_row(), mat1.num_row());
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 <= mat1.num_row();
382 snp1 += step1++)
383 {
384 for(step2 = 1, snp2 = mat2.m.begin(); step2 <= mat2.num_row();)
385 {
386 sp1 = snp1;
387 sp2 = snp2;
388 snp2 += step2;
389 temp = 0;
390 if(step1 < step2)
391 {
392 while(sp1 < snp1 + step1) // Loop checking, 06.08.2015, G.Cosmo
393 {
394 temp += (*(sp1++)) * (*(sp2++));
395 }
396 sp1 += step1 - 1;
397 for(stept1 = step1 + 1; stept1 != step2 + 1; sp1 += stept1++)
398 {
399 temp += (*sp1) * (*(sp2++));
400 }
401 sp2 += step2 - 1;
402 for(stept2 = ++step2; stept2 <= mat2.num_row();
403 sp1 += stept1++, sp2 += stept2++)
404 {
405 temp += (*sp1) * (*sp2);
406 }
407 }
408 else
409 {
410 while(sp2 < snp2) // Loop checking, 06.08.2015, G.Cosmo
411 {
412 temp += (*(sp1++)) * (*(sp2++));
413 }
414 sp2 += step2 - 1;
415 for(stept2 = ++step2; stept2 != step1 + 1; sp2 += stept2++)
416 {
417 temp += (*(sp1++)) * (*sp2);
418 }
419 sp1 += step1 - 1;
420 for(stept1 = step1 + 1; stept1 <= mat1.num_row();
421 sp1 += stept1++, sp2 += stept2++)
422 {
423 temp += (*sp1) * (*sp2);
424 }
425 }
426 *(mr++) = temp;
427 }
428 }
429 return mret;
430}

◆ operator+

G4ErrorSymMatrix operator+ ( const G4ErrorSymMatrix &  m1,
const G4ErrorSymMatrix &  m2 
)
friend

Definition at line 238 of file G4ErrorSymMatrix.cc.

240{
241 G4ErrorSymMatrix mret(mat1.nrow);
242 CHK_DIM_1(mat1.nrow, mat2.nrow, +);
243 SIMPLE_TOP(+)
244 return mret;
245}
SIMPLE_TOP
#define SIMPLE_TOP(OPER)
Definition: G4ErrorSymMatrix.cc:53

◆ operator-

G4ErrorSymMatrix operator- ( const G4ErrorSymMatrix &  m1,
const G4ErrorSymMatrix &  m2 
)
friend

Definition at line 267 of file G4ErrorSymMatrix.cc.

269{
270 G4ErrorSymMatrix mret(mat1.num_row());
271 CHK_DIM_1(mat1.num_row(), mat2.num_row(), -);
272 SIMPLE_TOP(-)
273 return mret;
274}

◆ tridiagonal

void tridiagonal ( G4ErrorSymMatrix *  a,
G4ErrorMatrix *  hsm 
)
friend

Field Documentation

◆ adjustment5x5

G4ThreadLocal G4double G4ErrorSymMatrix::adjustment5x5 = 0.0
staticprivate

Definition at line 225 of file G4ErrorSymMatrix.hh.

Referenced by invert5().

◆ adjustment6x6

G4ThreadLocal G4double G4ErrorSymMatrix::adjustment6x6 = 0.0
staticprivate

Definition at line 230 of file G4ErrorSymMatrix.hh.

Referenced by invert6().

◆ CHOLESKY_CREEP_5x5

const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_5x5 = .005
staticprivate

Definition at line 227 of file G4ErrorSymMatrix.hh.

Referenced by invert5().

◆ CHOLESKY_CREEP_6x6

const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_6x6 = .002
staticprivate

Definition at line 232 of file G4ErrorSymMatrix.hh.

Referenced by invert6().

◆ CHOLESKY_THRESHOLD_5x5

const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_5x5 = .5
staticprivate

Definition at line 226 of file G4ErrorSymMatrix.hh.

Referenced by invert5().

◆ CHOLESKY_THRESHOLD_6x6

const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_6x6 = .2
staticprivate

Definition at line 231 of file G4ErrorSymMatrix.hh.

Referenced by invert6().

◆ m

std::vector<G4double> G4ErrorSymMatrix::m
private

Definition at line 219 of file G4ErrorSymMatrix.hh.

Referenced by apply(), assign(), G4ErrorMatrix::G4ErrorMatrix(), G4ErrorSymMatrix(), invert(), invert4(), invertBunchKaufman(), invertCholesky5(), invertCholesky6(), invertHaywood5(), invertHaywood6(), G4ErrorMatrix::operator+=(), operator-(), G4ErrorMatrix::operator-=(), G4ErrorMatrix::operator=(), operator=(), similarity(), similarityT(), sub(), and trace().

◆ nrow

G4int G4ErrorSymMatrix::nrow
private

Definition at line 221 of file G4ErrorSymMatrix.hh.

Referenced by assign(), determinant(), G4ErrorSymMatrix(), invert(), invertBunchKaufman(), operator-(), G4ErrorMatrix::operator=(), operator=(), and trace().

◆ posDefFraction5x5

G4ThreadLocal G4double G4ErrorSymMatrix::posDefFraction5x5 = 1.0
staticprivate

Definition at line 224 of file G4ErrorSymMatrix.hh.

Referenced by invert5().

◆ posDefFraction6x6

G4ThreadLocal G4double G4ErrorSymMatrix::posDefFraction6x6 = 1.0
staticprivate

Definition at line 229 of file G4ErrorSymMatrix.hh.

Referenced by invert6().

◆ size

G4int G4ErrorSymMatrix::size
private

Definition at line 222 of file G4ErrorSymMatrix.hh.

Referenced by assign(), G4ErrorSymMatrix(), and operator=().

The documentation for this class was generated from the following files:
Generated by doxygen 1.9.3