Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Classes | Public Member Functions | Protected Member Functions | Friends | List of all members
G4ErrorSymMatrix Class Reference

#include <G4ErrorSymMatrix.hh>

Classes

class  G4ErrorSymMatrix_row
 
class  G4ErrorSymMatrix_row_const
 

Public Member Functions

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

Protected Member Functions

G4int num_size () const
 

Friends

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

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 70 of file G4ErrorSymMatrix.cc.

71 : m(p*(p+1)/2), nrow(p)
72{
73 size = nrow * (nrow+1) / 2;
74 m.assign(size,0);
75}

◆ G4ErrorSymMatrix() [3/5]

G4ErrorSymMatrix::G4ErrorSymMatrix ( G4int  p,
G4int  init 
)

Definition at line 77 of file G4ErrorSymMatrix.cc.

78 : m(p*(p+1)/2), nrow(p)
79{
80 size = nrow * (nrow+1) / 2;
81
82 m.assign(size,0);
83 switch(init)
84 {
85 case 0:
86 break;
87
88 case 1:
89 {
90 G4ErrorMatrixIter a = m.begin();
91 for(G4int i=1;i<=nrow;i++)
92 {
93 *a = 1.0;
94 a += (i+1);
95 }
96 break;
97 }
98 default:
99 G4ErrorMatrix::error("G4ErrorSymMatrix: initialization must be 0 or 1.");
100 }
101}
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:43
G4int
int G4int
Definition: G4Types.hh:85
G4ErrorMatrix::error
static void error(const char *s)
Definition: G4ErrorMatrix.cc:804

◆ G4ErrorSymMatrix() [4/5]

G4ErrorSymMatrix::G4ErrorSymMatrix ( const G4ErrorSymMatrix &  m1)

Definition at line 111 of file G4ErrorSymMatrix.cc.

112 : m(mat1.size), nrow(mat1.nrow), size(mat1.size)
113{
114 m = mat1.m;
115}

◆ G4ErrorSymMatrix() [5/5]

G4ErrorSymMatrix::G4ErrorSymMatrix ( G4ErrorSymMatrix &&  )
default

◆ ~G4ErrorSymMatrix()

G4ErrorSymMatrix::~G4ErrorSymMatrix ( )
virtual

Definition at line 107 of file G4ErrorSymMatrix.cc.

108{
109}

Member Function Documentation

◆ apply()

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

Definition at line 551 of file G4ErrorSymMatrix.cc.

553{
554 G4ErrorSymMatrix mret(num_row());
555 G4ErrorMatrixConstIter a = m.begin();
556 G4ErrorMatrixIter b = mret.m.begin();
557 for(G4int ir=1;ir<=num_row();ir++)
558 {
559 for(G4int ic=1;ic<=ir;ic++)
560 {
561 *(b++) = (*f)(*(a++), ir, ic);
562 }
563 }
564 return mret;
565}
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
G4ErrorSymMatrix::num_row
G4int num_row() const

◆ assign() [1/2]

void G4ErrorSymMatrix::assign ( const G4ErrorMatrix &  m2)

Definition at line 567 of file G4ErrorSymMatrix.cc.

568{
569 if(mat1.nrow != nrow)
570 {
571 nrow = mat1.nrow;
572 size = nrow * (nrow+1) / 2;
573 m.resize(size);
574 }
575 G4ErrorMatrixConstIter a = mat1.m.begin();
576 G4ErrorMatrixIter b = m.begin();
577 for(G4int r=1;r<=nrow;r++)
578 {
579 G4ErrorMatrixConstIter d = a;
580 for(G4int c=1;c<=r;c++)
581 {
582 *(b++) = *(d++);
583 }
584 a += nrow;
585 }
586}

◆ assign() [2/2]

void G4ErrorSymMatrix::assign ( const G4ErrorSymMatrix &  m2)

◆ determinant()

G4double G4ErrorSymMatrix::determinant ( ) const

Definition at line 798 of file G4ErrorSymMatrix.cc.

799{
800 static const G4int max_array = 20;
801
802 // ir must point to an array which is ***1 longer than*** nrow
803
804 static std::vector<G4int> ir_vec (max_array+1);
805 if (ir_vec.size() <= static_cast<unsigned int>(nrow))
806 {
807 ir_vec.resize(nrow+1);
808 }
809 G4int * ir = &ir_vec[0];
810
811 G4double det;
812 G4ErrorMatrix mt(*this);
813 G4int i = mt.dfact_matrix(det, ir);
814 if(i==0) { return det; }
815 return 0.0;
816}
G4double
double G4double
Definition: G4Types.hh:83

◆ 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 682 of file G4ErrorSymMatrix.cc.

683{
684 ifail = 0;
685
686 switch(nrow)
687 {
688 case 3:
689 {
690 G4double det, temp;
691 G4double t1, t2, t3;
692 G4double c11,c12,c13,c22,c23,c33;
693 c11 = (*(m.begin()+2)) * (*(m.begin()+5))
694 - (*(m.begin()+4)) * (*(m.begin()+4));
695 c12 = (*(m.begin()+4)) * (*(m.begin()+3))
696 - (*(m.begin()+1)) * (*(m.begin()+5));
697 c13 = (*(m.begin()+1)) * (*(m.begin()+4))
698 - (*(m.begin()+2)) * (*(m.begin()+3));
699 c22 = (*(m.begin()+5)) * (*m.begin())
700 - (*(m.begin()+3)) * (*(m.begin()+3));
701 c23 = (*(m.begin()+3)) * (*(m.begin()+1))
702 - (*(m.begin()+4)) * (*m.begin());
703 c33 = (*m.begin()) * (*(m.begin()+2))
704 - (*(m.begin()+1)) * (*(m.begin()+1));
705 t1 = std::fabs(*m.begin());
706 t2 = std::fabs(*(m.begin()+1));
707 t3 = std::fabs(*(m.begin()+3));
708 if (t1 >= t2)
709 {
710 if (t3 >= t1)
711 {
712 temp = *(m.begin()+3);
713 det = c23*c12-c22*c13;
714 }
715 else
716 {
717 temp = *m.begin();
718 det = c22*c33-c23*c23;
719 }
720 }
721 else if (t3 >= t2)
722 {
723 temp = *(m.begin()+3);
724 det = c23*c12-c22*c13;
725 }
726 else
727 {
728 temp = *(m.begin()+1);
729 det = c13*c23-c12*c33;
730 }
731 if (det==0)
732 {
733 ifail = 1;
734 return;
735 }
736 {
737 G4double ss = temp/det;
738 G4ErrorMatrixIter mq = m.begin();
739 *(mq++) = ss*c11;
740 *(mq++) = ss*c12;
741 *(mq++) = ss*c22;
742 *(mq++) = ss*c13;
743 *(mq++) = ss*c23;
744 *(mq) = ss*c33;
745 }
746 }
747 break;
748 case 2:
749 {
750 G4double det, temp, ss;
751 det = (*m.begin())*(*(m.begin()+2)) - (*(m.begin()+1))*(*(m.begin()+1));
752 if (det==0)
753 {
754 ifail = 1;
755 return;
756 }
757 ss = 1.0/det;
758 *(m.begin()+1) *= -ss;
759 temp = ss*(*(m.begin()+2));
760 *(m.begin()+2) = ss*(*m.begin());
761 *m.begin() = temp;
762 break;
763 }
764 case 1:
765 {
766 if ((*m.begin())==0)
767 {
768 ifail = 1;
769 return;
770 }
771 *m.begin() = 1.0/(*m.begin());
772 break;
773 }
774 case 5:
775 {
776 invert5(ifail);
777 return;
778 }
779 case 6:
780 {
781 invert6(ifail);
782 return;
783 }
784 case 4:
785 {
786 invert4(ifail);
787 return;
788 }
789 default:
790 {
791 invertBunchKaufman(ifail);
792 return;
793 }
794 }
795 return; // inversion successful
796}
G4ErrorSymMatrix::invertBunchKaufman
void invertBunchKaufman(G4int &ifail)
Definition: G4ErrorSymMatrix.cc:826

◆ invertBunchKaufman()

void G4ErrorSymMatrix::invertBunchKaufman ( G4int &  ifail)

Definition at line 826 of file G4ErrorSymMatrix.cc.

827{
828 // Bunch-Kaufman diagonal pivoting method
829 // It is decribed in J.R. Bunch, L. Kaufman (1977).
830 // "Some Stable Methods for Calculating Inertia and Solving Symmetric
831 // Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub,
832 // Charles F. van Loan, "Matrix Computations" (the second edition
833 // has a bug.) and implemented in "lapack"
834 // Mario Stanke, 09/97
835
836 G4int i, j, k, ss;
837 G4int pivrow;
838
839 // Establish the two working-space arrays needed: x and piv are
840 // used as pointers to arrays of doubles and ints respectively, each
841 // of length nrow. We do not want to reallocate each time through
842 // unless the size needs to grow. We do not want to leak memory, even
843 // by having a new without a delete that is only done once.
844
845 static const G4int max_array = 25;
846 static G4ThreadLocal std::vector<G4double> *xvec = 0;
847 if (!xvec) xvec = new std::vector<G4double> (max_array) ;
848 static G4ThreadLocal std::vector<G4int> *pivv = 0;
849 if (!pivv) pivv = new std::vector<G4int> (max_array) ;
850 typedef std::vector<G4int>::iterator pivIter;
851 if (xvec->size() < static_cast<unsigned int>(nrow)) xvec->resize(nrow);
852 if (pivv->size() < static_cast<unsigned int>(nrow)) pivv->resize(nrow);
853 // Note - resize should do nothing if the size is already larger than nrow,
854 // but on VC++ there are indications that it does so we check.
855 // Note - the data elements in a vector are guaranteed to be contiguous,
856 // so x[i] and piv[i] are optimally fast.
857 G4ErrorMatrixIter x = xvec->begin();
858 // x[i] is used as helper storage, needs to have at least size nrow.
859 pivIter piv = pivv->begin();
860 // piv[i] is used to store details of exchanges
861
862 G4double temp1, temp2;
863 G4ErrorMatrixIter ip, mjj, iq;
864 G4double lambda, sigma;
865 const G4double alpha = .6404; // = (1+sqrt(17))/8
866 const G4double epsilon = 32*DBL_EPSILON;
867 // whenever a sum of two doubles is below or equal to epsilon
868 // it is set to zero.
869 // this constant could be set to zero but then the algorithm
870 // doesn't neccessarily detect that a matrix is singular
871
872 for (i = 0; i < nrow; i++)
873 {
874 piv[i] = i+1;
875 }
876
877 ifail = 0;
878
879 // compute the factorization P*A*P^T = L * D * L^T
880 // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
881 // L and D^-1 are stored in A = *this, P is stored in piv[]
882
883 for (j=1; j < nrow; j+=ss) // main loop over columns
884 {
885 mjj = m.begin() + j*(j-1)/2 + j-1;
886 lambda = 0; // compute lambda = max of A(j+1:n,j)
887 pivrow = j+1;
888 ip = m.begin() + (j+1)*j/2 + j-1;
889 for (i=j+1; i <= nrow ; ip += i++)
890 {
891 if (std::fabs(*ip) > lambda)
892 {
893 lambda = std::fabs(*ip);
894 pivrow = i;
895 }
896 }
897 if (lambda == 0 )
898 {
899 if (*mjj == 0)
900 {
901 ifail = 1;
902 return;
903 }
904 ss=1;
905 *mjj = 1./ *mjj;
906 }
907 else
908 {
909 if (std::fabs(*mjj) >= lambda*alpha)
910 {
911 ss=1;
912 pivrow=j;
913 }
914 else
915 {
916 sigma = 0; // compute sigma = max A(pivrow, j:pivrow-1)
917 ip = m.begin() + pivrow*(pivrow-1)/2+j-1;
918 for (k=j; k < pivrow; k++)
919 {
920 if (std::fabs(*ip) > sigma)
921 sigma = std::fabs(*ip);
922 ip++;
923 }
924 if (sigma * std::fabs(*mjj) >= alpha * lambda * lambda)
925 {
926 ss=1;
927 pivrow = j;
928 }
929 else if (std::fabs(*(m.begin()+pivrow*(pivrow-1)/2+pivrow-1))
930 >= alpha * sigma)
931 { ss=1; }
932 else
933 { ss=2; }
934 }
935 if (pivrow == j) // no permutation neccessary
936 {
937 piv[j-1] = pivrow;
938 if (*mjj == 0)
939 {
940 ifail=1;
941 return;
942 }
943 temp2 = *mjj = 1./ *mjj; // invert D(j,j)
944
945 // update A(j+1:n, j+1,n)
946 for (i=j+1; i <= nrow; i++)
947 {
948 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
949 ip = m.begin()+i*(i-1)/2+j;
950 for (k=j+1; k<=i; k++)
951 {
952 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
953 if (std::fabs(*ip) <= epsilon)
954 { *ip=0; }
955 ip++;
956 }
957 }
958 // update L
959 ip = m.begin() + (j+1)*j/2 + j-1;
960 for (i=j+1; i <= nrow; ip += i++)
961 {
962 *ip *= temp2;
963 }
964 }
965 else if (ss==1) // 1x1 pivot
966 {
967 piv[j-1] = pivrow;
968
969 // interchange rows and columns j and pivrow in
970 // submatrix (j:n,j:n)
971 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
972 for (i=j+1; i < pivrow; i++, ip++)
973 {
974 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
975 *(m.begin() + i*(i-1)/2 + j-1)= *ip;
976 *ip = temp1;
977 }
978 temp1 = *mjj;
979 *mjj = *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1);
980 *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1) = temp1;
981 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;
982 iq = ip + pivrow-j;
983 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
984 {
985 temp1 = *iq;
986 *iq = *ip;
987 *ip = temp1;
988 }
989
990 if (*mjj == 0)
991 {
992 ifail = 1;
993 return;
994 }
995 temp2 = *mjj = 1./ *mjj; // invert D(j,j)
996
997 // update A(j+1:n, j+1:n)
998 for (i = j+1; i <= nrow; i++)
999 {
1000 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
1001 ip = m.begin()+i*(i-1)/2+j;
1002 for (k=j+1; k<=i; k++)
1003 {
1004 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
1005 if (std::fabs(*ip) <= epsilon)
1006 { *ip=0; }
1007 ip++;
1008 }
1009 }
1010 // update L
1011 ip = m.begin() + (j+1)*j/2 + j-1;
1012 for (i=j+1; i<=nrow; ip += i++)
1013 {
1014 *ip *= temp2;
1015 }
1016 }
1017 else // ss=2, ie use a 2x2 pivot
1018 {
1019 piv[j-1] = -pivrow;
1020 piv[j] = 0; // that means this is the second row of a 2x2 pivot
1021
1022 if (j+1 != pivrow)
1023 {
1024 // interchange rows and columns j+1 and pivrow in
1025 // submatrix (j:n,j:n)
1026 ip = m.begin() + pivrow*(pivrow-1)/2 + j+1;
1027 for (i=j+2; i < pivrow; i++, ip++)
1028 {
1029 temp1 = *(m.begin() + i*(i-1)/2 + j);
1030 *(m.begin() + i*(i-1)/2 + j) = *ip;
1031 *ip = temp1;
1032 }
1033 temp1 = *(mjj + j + 1);
1034 *(mjj + j + 1) =
1035 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1036 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1037 temp1 = *(mjj + j);
1038 *(mjj + j) = *(m.begin() + pivrow*(pivrow-1)/2 + j-1);
1039 *(m.begin() + pivrow*(pivrow-1)/2 + j-1) = temp1;
1040 ip = m.begin() + (pivrow+1)*pivrow/2 + j;
1041 iq = ip + pivrow-(j+1);
1042 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
1043 {
1044 temp1 = *iq;
1045 *iq = *ip;
1046 *ip = temp1;
1047 }
1048 }
1049 // invert D(j:j+1,j:j+1)
1050 temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j);
1051 if (temp2 == 0)
1052 {
1053 G4Exception("G4ErrorSymMatrix::bunch_invert()",
1054 "GEANT4e-Notification", JustWarning,
1055 "Error in pivot choice!");
1056 }
1057 temp2 = 1. / temp2;
1058
1059 // this quotient is guaranteed to exist by the choice
1060 // of the pivot
1061
1062 temp1 = *mjj;
1063 *mjj = *(mjj + j + 1) * temp2;
1064 *(mjj + j + 1) = temp1 * temp2;
1065 *(mjj + j) = - *(mjj + j) * temp2;
1066
1067 if (j < nrow-1) // otherwise do nothing
1068 {
1069 // update A(j+2:n, j+2:n)
1070 for (i=j+2; i <= nrow ; i++)
1071 {
1072 ip = m.begin() + i*(i-1)/2 + j-1;
1073 temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
1074 if (std::fabs(temp1 ) <= epsilon)
1075 { temp1 = 0; }
1076 temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
1077 if (std::fabs(temp2 ) <= epsilon)
1078 { temp2 = 0; }
1079 for (k = j+2; k <= i ; k++)
1080 {
1081 ip = m.begin() + i*(i-1)/2 + k-1;
1082 iq = m.begin() + k*(k-1)/2 + j-1;
1083 *ip -= temp1 * *iq + temp2 * *(iq+1);
1084 if (std::fabs(*ip) <= epsilon)
1085 { *ip = 0; }
1086 }
1087 }
1088 // update L
1089 for (i=j+2; i <= nrow ; i++)
1090 {
1091 ip = m.begin() + i*(i-1)/2 + j-1;
1092 temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
1093 if (std::fabs(temp1) <= epsilon)
1094 { temp1 = 0; }
1095 *(ip+1) = *ip * *(mjj + j) + *(ip+1) * *(mjj + j + 1);
1096 if (std::fabs(*(ip+1)) <= epsilon)
1097 { *(ip+1) = 0; }
1098 *ip = temp1;
1099 }
1100 }
1101 }
1102 }
1103 } // end of main loop over columns
1104
1105 if (j == nrow) // the the last pivot is 1x1
1106 {
1107 mjj = m.begin() + j*(j-1)/2 + j-1;
1108 if (*mjj == 0)
1109 {
1110 ifail = 1;
1111 return;
1112 }
1113 else
1114 {
1115 *mjj = 1. / *mjj;
1116 }
1117 } // end of last pivot code
1118
1119 // computing the inverse from the factorization
1120
1121 for (j = nrow ; j >= 1 ; j -= ss) // loop over columns
1122 {
1123 mjj = m.begin() + j*(j-1)/2 + j-1;
1124 if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse
1125 {
1126 ss = 1;
1127 if (j < nrow)
1128 {
1129 ip = m.begin() + (j+1)*j/2 + j-1;
1130 for (i=0; i < nrow-j; ip += 1+j+i++)
1131 {
1132 x[i] = *ip;
1133 }
1134 for (i=j+1; i<=nrow ; i++)
1135 {
1136 temp2=0;
1137 ip = m.begin() + i*(i-1)/2 + j;
1138 for (k=0; k <= i-j-1; k++)
1139 { temp2 += *ip++ * x[k]; }
1140 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1141 { temp2 += *ip * x[k]; }
1142 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1143 }
1144 temp2 = 0;
1145 ip = m.begin() + (j+1)*j/2 + j-1;
1146 for (k=0; k < nrow-j; ip += 1+j+k++)
1147 { temp2 += x[k] * *ip; }
1148 *mjj -= temp2;
1149 }
1150 }
1151 else //2x2 pivot, compute columns j and j-1 of the inverse
1152 {
1153 if (piv[j-1] != 0)
1154 {
1155 std::ostringstream message;
1156 message << "Error in pivot: " << piv[j-1];
1157 G4Exception("G4ErrorSymMatrix::invertBunchKaufman()",
1158 "GEANT4e-Notification", JustWarning, message);
1159 }
1160 ss=2;
1161 if (j < nrow)
1162 {
1163 ip = m.begin() + (j+1)*j/2 + j-1;
1164 for (i=0; i < nrow-j; ip += 1+j+i++)
1165 {
1166 x[i] = *ip;
1167 }
1168 for (i=j+1; i<=nrow ; i++)
1169 {
1170 temp2 = 0;
1171 ip = m.begin() + i*(i-1)/2 + j;
1172 for (k=0; k <= i-j-1; k++)
1173 { temp2 += *ip++ * x[k]; }
1174 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1175 { temp2 += *ip * x[k]; }
1176 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1177 }
1178 temp2 = 0;
1179 ip = m.begin() + (j+1)*j/2 + j-1;
1180 for (k=0; k < nrow-j; ip += 1+j+k++)
1181 { temp2 += x[k] * *ip; }
1182 *mjj -= temp2;
1183 temp2 = 0;
1184 ip = m.begin() + (j+1)*j/2 + j-2;
1185 for (i=j+1; i <= nrow; ip += i++)
1186 { temp2 += *ip * *(ip+1); }
1187 *(mjj-1) -= temp2;
1188 ip = m.begin() + (j+1)*j/2 + j-2;
1189 for (i=0; i < nrow-j; ip += 1+j+i++)
1190 {
1191 x[i] = *ip;
1192 }
1193 for (i=j+1; i <= nrow ; i++)
1194 {
1195 temp2 = 0;
1196 ip = m.begin() + i*(i-1)/2 + j;
1197 for (k=0; k <= i-j-1; k++)
1198 { temp2 += *ip++ * x[k]; }
1199 for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1200 { temp2 += *ip * x[k]; }
1201 *(m.begin()+ i*(i-1)/2 + j-2)= -temp2;
1202 }
1203 temp2 = 0;
1204 ip = m.begin() + (j+1)*j/2 + j-2;
1205 for (k=0; k < nrow-j; ip += 1+j+k++)
1206 { temp2 += x[k] * *ip; }
1207 *(mjj-j) -= temp2;
1208 }
1209 }
1210
1211 // interchange rows and columns j and piv[j-1]
1212 // or rows and columns j and -piv[j-2]
1213
1214 pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
1215 ip = m.begin() + pivrow*(pivrow-1)/2 + j;
1216 for (i=j+1;i < pivrow; i++, ip++)
1217 {
1218 temp1 = *(m.begin() + i*(i-1)/2 + j-1);
1219 *(m.begin() + i*(i-1)/2 + j-1) = *ip;
1220 *ip = temp1;
1221 }
1222 temp1 = *mjj;
1223 *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1224 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1225 if (ss==2)
1226 {
1227 temp1 = *(mjj-1);
1228 *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2);
1229 *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1;
1230 }
1231
1232 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1; // &A(i,j)
1233 iq = ip + pivrow-j;
1234 for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
1235 {
1236 temp1 = *iq;
1237 *iq = *ip;
1238 *ip = temp1;
1239 }
1240 } // end of loop over columns (in computing inverse from factorization)
1241
1242 return; // inversion successful
1243}
epsilon
double epsilon(double density, double temperature)
Definition: G4DNAElectronHoleRecombination.cc:96
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:65
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
DBL_EPSILON
#define DBL_EPSILON
Definition: templates.hh:66
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:77

Referenced by invert().

◆ invertCholesky5()

void G4ErrorSymMatrix::invertCholesky5 ( G4int &  ifail)

Definition at line 1917 of file G4ErrorSymMatrix.cc.

1918{
1919
1920 // Invert by
1921 //
1922 // a) decomposing M = G*G^T with G lower triangular
1923 // (if M is not positive definite this will fail, leaving this unchanged)
1924 // b) inverting G to form H
1925 // c) multiplying H^T * H to get M^-1.
1926 //
1927 // If the matrix is pos. def. it is inverted and 1 is returned.
1928 // If the matrix is not pos. def. it remains unaltered and 0 is returned.
1929
1930 G4double h10; // below-diagonal elements of H
1931 G4double h20, h21;
1932 G4double h30, h31, h32;
1933 G4double h40, h41, h42, h43;
1934
1935 G4double h00, h11, h22, h33, h44; // 1/diagonal elements of G =
1936 // diagonal elements of H
1937
1938 G4double g10; // below-diagonal elements of G
1939 G4double g20, g21;
1940 G4double g30, g31, g32;
1941 G4double g40, g41, g42, g43;
1942
1943 ifail = 1; // We start by assuing we won't succeed...
1944
1945 // Form G -- compute diagonal members of H directly rather than of G
1946 //-------
1947
1948 // Scale first column by 1/sqrt(A00)
1949
1950 h00 = m[A00];
1951 if (h00 <= 0) { return; }
1952 h00 = 1.0 / std::sqrt(h00);
1953
1954 g10 = m[A10] * h00;
1955 g20 = m[A20] * h00;
1956 g30 = m[A30] * h00;
1957 g40 = m[A40] * h00;
1958
1959 // Form G11 (actually, just h11)
1960
1961 h11 = m[A11] - (g10 * g10);
1962 if (h11 <= 0) { return; }
1963 h11 = 1.0 / std::sqrt(h11);
1964
1965 // Subtract inter-column column dot products from rest of column 1 and
1966 // scale to get column 1 of G
1967
1968 g21 = (m[A21] - (g10 * g20)) * h11;
1969 g31 = (m[A31] - (g10 * g30)) * h11;
1970 g41 = (m[A41] - (g10 * g40)) * h11;
1971
1972 // Form G22 (actually, just h22)
1973
1974 h22 = m[A22] - (g20 * g20) - (g21 * g21);
1975 if (h22 <= 0) { return; }
1976 h22 = 1.0 / std::sqrt(h22);
1977
1978 // Subtract inter-column column dot products from rest of column 2 and
1979 // scale to get column 2 of G
1980
1981 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
1982 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
1983
1984 // Form G33 (actually, just h33)
1985
1986 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
1987 if (h33 <= 0) { return; }
1988 h33 = 1.0 / std::sqrt(h33);
1989
1990 // Subtract inter-column column dot product from A43 and scale to get G43
1991
1992 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
1993
1994 // Finally form h44 - if this is possible inversion succeeds
1995
1996 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
1997 if (h44 <= 0) { return; }
1998 h44 = 1.0 / std::sqrt(h44);
1999
2000 // Form H = 1/G -- diagonal members of H are already correct
2001 //-------------
2002
2003 // The order here is dictated by speed considerations
2004
2005 h43 = -h33 * g43 * h44;
2006 h32 = -h22 * g32 * h33;
2007 h42 = -h22 * (g32 * h43 + g42 * h44);
2008 h21 = -h11 * g21 * h22;
2009 h31 = -h11 * (g21 * h32 + g31 * h33);
2010 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2011 h10 = -h00 * g10 * h11;
2012 h20 = -h00 * (g10 * h21 + g20 * h22);
2013 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2014 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2015
2016 // Change this to its inverse = H^T*H
2017 //------------------------------------
2018
2019 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40;
2020 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41;
2021 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41;
2022 m[A02] = h20 * h22 + h30 * h32 + h40 * h42;
2023 m[A12] = h21 * h22 + h31 * h32 + h41 * h42;
2024 m[A22] = h22 * h22 + h32 * h32 + h42 * h42;
2025 m[A03] = h30 * h33 + h40 * h43;
2026 m[A13] = h31 * h33 + h41 * h43;
2027 m[A23] = h32 * h33 + h42 * h43;
2028 m[A33] = h33 * h33 + h43 * h43;
2029 m[A04] = h40 * h44;
2030 m[A14] = h41 * h44;
2031 m[A24] = h42 * h44;
2032 m[A34] = h43 * h44;
2033 m[A44] = h44 * h44;
2034
2035 ifail = 0;
2036 return;
2037}
A41
#define A41
Definition: G4ErrorMatrix.cc:846
A20
#define A20
Definition: G4ErrorMatrix.cc:831
A01
#define A01
Definition: G4ErrorMatrix.cc:818
A23
#define A23
Definition: G4ErrorMatrix.cc:834
A13
#define A13
Definition: G4ErrorMatrix.cc:827
A00
#define A00
Definition: G4ErrorMatrix.cc:817
A40
#define A40
Definition: G4ErrorMatrix.cc:845
A02
#define A02
Definition: G4ErrorMatrix.cc:819
A24
#define A24
Definition: G4ErrorMatrix.cc:835
A22
#define A22
Definition: G4ErrorMatrix.cc:833
A04
#define A04
Definition: G4ErrorMatrix.cc:821
A30
#define A30
Definition: G4ErrorMatrix.cc:838
A12
#define A12
Definition: G4ErrorMatrix.cc:826
A03
#define A03
Definition: G4ErrorMatrix.cc:820
A14
#define A14
Definition: G4ErrorMatrix.cc:828
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
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
A34
#define A34
Definition: G4ErrorMatrix.cc:842
A43
#define A43
Definition: G4ErrorMatrix.cc:848

◆ invertCholesky6()

void G4ErrorSymMatrix::invertCholesky6 ( G4int &  ifail)

Definition at line 2039 of file G4ErrorSymMatrix.cc.

2040{
2041 // Invert by
2042 //
2043 // a) decomposing M = G*G^T with G lower triangular
2044 // (if M is not positive definite this will fail, leaving this unchanged)
2045 // b) inverting G to form H
2046 // c) multiplying H^T * H to get M^-1.
2047 //
2048 // If the matrix is pos. def. it is inverted and 1 is returned.
2049 // If the matrix is not pos. def. it remains unaltered and 0 is returned.
2050
2051 G4double h10; // below-diagonal elements of H
2052 G4double h20, h21;
2053 G4double h30, h31, h32;
2054 G4double h40, h41, h42, h43;
2055 G4double h50, h51, h52, h53, h54;
2056
2057 G4double h00, h11, h22, h33, h44, h55; // 1/diagonal elements of G =
2058 // diagonal elements of H
2059
2060 G4double g10; // below-diagonal elements of G
2061 G4double g20, g21;
2062 G4double g30, g31, g32;
2063 G4double g40, g41, g42, g43;
2064 G4double g50, g51, g52, g53, g54;
2065
2066 ifail = 1; // We start by assuing we won't succeed...
2067
2068 // Form G -- compute diagonal members of H directly rather than of G
2069 //-------
2070
2071 // Scale first column by 1/sqrt(A00)
2072
2073 h00 = m[A00];
2074 if (h00 <= 0) { return; }
2075 h00 = 1.0 / std::sqrt(h00);
2076
2077 g10 = m[A10] * h00;
2078 g20 = m[A20] * h00;
2079 g30 = m[A30] * h00;
2080 g40 = m[A40] * h00;
2081 g50 = m[A50] * h00;
2082
2083 // Form G11 (actually, just h11)
2084
2085 h11 = m[A11] - (g10 * g10);
2086 if (h11 <= 0) { return; }
2087 h11 = 1.0 / std::sqrt(h11);
2088
2089 // Subtract inter-column column dot products from rest of column 1 and
2090 // scale to get column 1 of G
2091
2092 g21 = (m[A21] - (g10 * g20)) * h11;
2093 g31 = (m[A31] - (g10 * g30)) * h11;
2094 g41 = (m[A41] - (g10 * g40)) * h11;
2095 g51 = (m[A51] - (g10 * g50)) * h11;
2096
2097 // Form G22 (actually, just h22)
2098
2099 h22 = m[A22] - (g20 * g20) - (g21 * g21);
2100 if (h22 <= 0) { return; }
2101 h22 = 1.0 / std::sqrt(h22);
2102
2103 // Subtract inter-column column dot products from rest of column 2 and
2104 // scale to get column 2 of G
2105
2106 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
2107 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
2108 g52 = (m[A52] - (g20 * g50) - (g21 * g51)) * h22;
2109
2110 // Form G33 (actually, just h33)
2111
2112 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
2113 if (h33 <= 0) { return; }
2114 h33 = 1.0 / std::sqrt(h33);
2115
2116 // Subtract inter-column column dot products from rest of column 3 and
2117 // scale to get column 3 of G
2118
2119 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
2120 g53 = (m[A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33;
2121
2122 // Form G44 (actually, just h44)
2123
2124 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
2125 if (h44 <= 0) { return; }
2126 h44 = 1.0 / std::sqrt(h44);
2127
2128 // Subtract inter-column column dot product from M54 and scale to get G54
2129
2130 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;
2131
2132 // Finally form h55 - if this is possible inversion succeeds
2133
2134 h55 = m[A55] - (g50*g50) - (g51*g51) - (g52*g52) - (g53*g53) - (g54*g54);
2135 if (h55 <= 0) { return; }
2136 h55 = 1.0 / std::sqrt(h55);
2137
2138 // Form H = 1/G -- diagonal members of H are already correct
2139 //-------------
2140
2141 // The order here is dictated by speed considerations
2142
2143 h54 = -h44 * g54 * h55;
2144 h43 = -h33 * g43 * h44;
2145 h53 = -h33 * (g43 * h54 + g53 * h55);
2146 h32 = -h22 * g32 * h33;
2147 h42 = -h22 * (g32 * h43 + g42 * h44);
2148 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * h55);
2149 h21 = -h11 * g21 * h22;
2150 h31 = -h11 * (g21 * h32 + g31 * h33);
2151 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
2152 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55);
2153 h10 = -h00 * g10 * h11;
2154 h20 = -h00 * (g10 * h21 + g20 * h22);
2155 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
2156 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
2157 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55);
2158
2159 // Change this to its inverse = H^T*H
2160 //------------------------------------
2161
2162 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40 + h50*h50;
2163 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51;
2164 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51;
2165 m[A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52;
2166 m[A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52;
2167 m[A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52;
2168 m[A03] = h30 * h33 + h40 * h43 + h50 * h53;
2169 m[A13] = h31 * h33 + h41 * h43 + h51 * h53;
2170 m[A23] = h32 * h33 + h42 * h43 + h52 * h53;
2171 m[A33] = h33 * h33 + h43 * h43 + h53 * h53;
2172 m[A04] = h40 * h44 + h50 * h54;
2173 m[A14] = h41 * h44 + h51 * h54;
2174 m[A24] = h42 * h44 + h52 * h54;
2175 m[A34] = h43 * h44 + h53 * h54;
2176 m[A44] = h44 * h44 + h54 * h54;
2177 m[A05] = h50 * h55;
2178 m[A15] = h51 * h55;
2179 m[A25] = h52 * h55;
2180 m[A35] = h53 * h55;
2181 m[A45] = h54 * h55;
2182 m[A55] = h55 * h55;
2183
2184 ifail = 0;
2185 return;
2186}
A53
#define A53
Definition: G4ErrorMatrix.cc:855
A45
#define A45
Definition: G4ErrorMatrix.cc:850
A54
#define A54
Definition: G4ErrorMatrix.cc:856
A25
#define A25
Definition: G4ErrorMatrix.cc:836
A55
#define A55
Definition: G4ErrorMatrix.cc:857
A35
#define A35
Definition: G4ErrorMatrix.cc:843
A50
#define A50
Definition: G4ErrorMatrix.cc:852
A51
#define A51
Definition: G4ErrorMatrix.cc:853
A05
#define A05
Definition: G4ErrorMatrix.cc:822
A52
#define A52
Definition: G4ErrorMatrix.cc:854
A15
#define A15
Definition: G4ErrorMatrix.cc:829

◆ invertHaywood4()

void G4ErrorSymMatrix::invertHaywood4 ( G4int &  ifail)

Definition at line 2266 of file G4ErrorSymMatrix.cc.

2267{
2268 invert4(ifail); // For the 4x4 case, the method we use for invert is already
2269 // the Haywood method.
2270}

◆ invertHaywood5()

void G4ErrorSymMatrix::invertHaywood5 ( G4int &  ifail)

Definition at line 1364 of file G4ErrorSymMatrix.cc.

1365{
1366 ifail = 0;
1367
1368 // Find all NECESSARY 2x2 dets: (25 of them)
1369
1370 G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
1371 G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
1372 G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
1373 G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
1374 G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
1375 G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
1376 G4double Det2_24_01 = m[A20]*m[A41] - m[A21]*m[A40];
1377 G4double Det2_24_02 = m[A20]*m[A42] - m[A22]*m[A40];
1378 G4double Det2_24_03 = m[A20]*m[A43] - m[A23]*m[A40];
1379 G4double Det2_24_04 = m[A20]*m[A44] - m[A24]*m[A40];
1380 G4double Det2_24_12 = m[A21]*m[A42] - m[A22]*m[A41];
1381 G4double Det2_24_13 = m[A21]*m[A43] - m[A23]*m[A41];
1382 G4double Det2_24_14 = m[A21]*m[A44] - m[A24]*m[A41];
1383 G4double Det2_24_23 = m[A22]*m[A43] - m[A23]*m[A42];
1384 G4double Det2_24_24 = m[A22]*m[A44] - m[A24]*m[A42];
1385 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1386 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1387 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1388 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1389 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1390 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1391 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1392 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1393 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1394 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1395
1396 // Find all NECESSARY 3x3 dets: (30 of them)
1397
1398 G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02
1399 + m[A12]*Det2_23_01;
1400 G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03
1401 + m[A13]*Det2_23_01;
1402 G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03
1403 + m[A13]*Det2_23_02;
1404 G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13
1405 + m[A13]*Det2_23_12;
1406 G4double Det3_124_012 = m[A10]*Det2_24_12 - m[A11]*Det2_24_02
1407 + m[A12]*Det2_24_01;
1408 G4double Det3_124_013 = m[A10]*Det2_24_13 - m[A11]*Det2_24_03
1409 + m[A13]*Det2_24_01;
1410 G4double Det3_124_014 = m[A10]*Det2_24_14 - m[A11]*Det2_24_04
1411 + m[A14]*Det2_24_01;
1412 G4double Det3_124_023 = m[A10]*Det2_24_23 - m[A12]*Det2_24_03
1413 + m[A13]*Det2_24_02;
1414 G4double Det3_124_024 = m[A10]*Det2_24_24 - m[A12]*Det2_24_04
1415 + m[A14]*Det2_24_02;
1416 G4double Det3_124_123 = m[A11]*Det2_24_23 - m[A12]*Det2_24_13
1417 + m[A13]*Det2_24_12;
1418 G4double Det3_124_124 = m[A11]*Det2_24_24 - m[A12]*Det2_24_14
1419 + m[A14]*Det2_24_12;
1420 G4double Det3_134_012 = m[A10]*Det2_34_12 - m[A11]*Det2_34_02
1421 + m[A12]*Det2_34_01;
1422 G4double Det3_134_013 = m[A10]*Det2_34_13 - m[A11]*Det2_34_03
1423 + m[A13]*Det2_34_01;
1424 G4double Det3_134_014 = m[A10]*Det2_34_14 - m[A11]*Det2_34_04
1425 + m[A14]*Det2_34_01;
1426 G4double Det3_134_023 = m[A10]*Det2_34_23 - m[A12]*Det2_34_03
1427 + m[A13]*Det2_34_02;
1428 G4double Det3_134_024 = m[A10]*Det2_34_24 - m[A12]*Det2_34_04
1429 + m[A14]*Det2_34_02;
1430 G4double Det3_134_034 = m[A10]*Det2_34_34 - m[A13]*Det2_34_04
1431 + m[A14]*Det2_34_03;
1432 G4double Det3_134_123 = m[A11]*Det2_34_23 - m[A12]*Det2_34_13
1433 + m[A13]*Det2_34_12;
1434 G4double Det3_134_124 = m[A11]*Det2_34_24 - m[A12]*Det2_34_14
1435 + m[A14]*Det2_34_12;
1436 G4double Det3_134_134 = m[A11]*Det2_34_34 - m[A13]*Det2_34_14
1437 + m[A14]*Det2_34_13;
1438 G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
1439 + m[A22]*Det2_34_01;
1440 G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
1441 + m[A23]*Det2_34_01;
1442 G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
1443 + m[A24]*Det2_34_01;
1444 G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
1445 + m[A23]*Det2_34_02;
1446 G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
1447 + m[A24]*Det2_34_02;
1448 G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
1449 + m[A24]*Det2_34_03;
1450 G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
1451 + m[A23]*Det2_34_12;
1452 G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
1453 + m[A24]*Det2_34_12;
1454 G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
1455 + m[A24]*Det2_34_13;
1456 G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
1457 + m[A24]*Det2_34_23;
1458
1459 // Find all NECESSARY 4x4 dets: (15 of them)
1460
1461 G4double Det4_0123_0123 = m[A00]*Det3_123_123 - m[A01]*Det3_123_023
1462 + m[A02]*Det3_123_013 - m[A03]*Det3_123_012;
1463 G4double Det4_0124_0123 = m[A00]*Det3_124_123 - m[A01]*Det3_124_023
1464 + m[A02]*Det3_124_013 - m[A03]*Det3_124_012;
1465 G4double Det4_0124_0124 = m[A00]*Det3_124_124 - m[A01]*Det3_124_024
1466 + m[A02]*Det3_124_014 - m[A04]*Det3_124_012;
1467 G4double Det4_0134_0123 = m[A00]*Det3_134_123 - m[A01]*Det3_134_023
1468 + m[A02]*Det3_134_013 - m[A03]*Det3_134_012;
1469 G4double Det4_0134_0124 = m[A00]*Det3_134_124 - m[A01]*Det3_134_024
1470 + m[A02]*Det3_134_014 - m[A04]*Det3_134_012;
1471 G4double Det4_0134_0134 = m[A00]*Det3_134_134 - m[A01]*Det3_134_034
1472 + m[A03]*Det3_134_014 - m[A04]*Det3_134_013;
1473 G4double Det4_0234_0123 = m[A00]*Det3_234_123 - m[A01]*Det3_234_023
1474 + m[A02]*Det3_234_013 - m[A03]*Det3_234_012;
1475 G4double Det4_0234_0124 = m[A00]*Det3_234_124 - m[A01]*Det3_234_024
1476 + m[A02]*Det3_234_014 - m[A04]*Det3_234_012;
1477 G4double Det4_0234_0134 = m[A00]*Det3_234_134 - m[A01]*Det3_234_034
1478 + m[A03]*Det3_234_014 - m[A04]*Det3_234_013;
1479 G4double Det4_0234_0234 = m[A00]*Det3_234_234 - m[A02]*Det3_234_034
1480 + m[A03]*Det3_234_024 - m[A04]*Det3_234_023;
1481 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
1482 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1483 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
1484 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1485 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
1486 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1487 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
1488 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1489 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
1490 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1491
1492 // Find the 5x5 det:
1493
1494 G4double det = m[A00]*Det4_1234_1234
1495 - m[A01]*Det4_1234_0234
1496 + m[A02]*Det4_1234_0134
1497 - m[A03]*Det4_1234_0124
1498 + m[A04]*Det4_1234_0123;
1499
1500 if ( det == 0 )
1501 {
1502 ifail = 1;
1503 return;
1504 }
1505
1506 G4double oneOverDet = 1.0/det;
1507 G4double mn1OverDet = - oneOverDet;
1508
1509 m[A00] = Det4_1234_1234 * oneOverDet;
1510 m[A01] = Det4_1234_0234 * mn1OverDet;
1511 m[A02] = Det4_1234_0134 * oneOverDet;
1512 m[A03] = Det4_1234_0124 * mn1OverDet;
1513 m[A04] = Det4_1234_0123 * oneOverDet;
1514
1515 m[A11] = Det4_0234_0234 * oneOverDet;
1516 m[A12] = Det4_0234_0134 * mn1OverDet;
1517 m[A13] = Det4_0234_0124 * oneOverDet;
1518 m[A14] = Det4_0234_0123 * mn1OverDet;
1519
1520 m[A22] = Det4_0134_0134 * oneOverDet;
1521 m[A23] = Det4_0134_0124 * mn1OverDet;
1522 m[A24] = Det4_0134_0123 * oneOverDet;
1523
1524 m[A33] = Det4_0124_0124 * oneOverDet;
1525 m[A34] = Det4_0124_0123 * mn1OverDet;
1526
1527 m[A44] = Det4_0123_0123 * oneOverDet;
1528
1529 return;
1530}

◆ invertHaywood6()

void G4ErrorSymMatrix::invertHaywood6 ( G4int &  ifail)

Definition at line 1532 of file G4ErrorSymMatrix.cc.

1533{
1534 ifail = 0;
1535
1536 // Find all NECESSARY 2x2 dets: (39 of them)
1537
1538 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
1539 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
1540 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
1541 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
1542 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
1543 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
1544 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
1545 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
1546 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
1547 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
1548 G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
1549 G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
1550 G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
1551 G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
1552 G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
1553 G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
1554 G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
1555 G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
1556 G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
1557 G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
1558 G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
1559 G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
1560 G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
1561 G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
1562 G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
1563 G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
1564 G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
1565 G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
1566 G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
1567 G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
1568 G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
1569 G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
1570 G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
1571 G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
1572 G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
1573 G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
1574 G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
1575 G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
1576 G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
1577
1578 // Find all NECESSARY 3x3 dets: (65 of them)
1579
1580 G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02
1581 + m[A22]*Det2_34_01;
1582 G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03
1583 + m[A23]*Det2_34_01;
1584 G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04
1585 + m[A24]*Det2_34_01;
1586 G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03
1587 + m[A23]*Det2_34_02;
1588 G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04
1589 + m[A24]*Det2_34_02;
1590 G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04
1591 + m[A24]*Det2_34_03;
1592 G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13
1593 + m[A23]*Det2_34_12;
1594 G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14
1595 + m[A24]*Det2_34_12;
1596 G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14
1597 + m[A24]*Det2_34_13;
1598 G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24
1599 + m[A24]*Det2_34_23;
1600 G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02
1601 + m[A22]*Det2_35_01;
1602 G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03
1603 + m[A23]*Det2_35_01;
1604 G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04
1605 + m[A24]*Det2_35_01;
1606 G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05
1607 + m[A25]*Det2_35_01;
1608 G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03
1609 + m[A23]*Det2_35_02;
1610 G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04
1611 + m[A24]*Det2_35_02;
1612 G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05
1613 + m[A25]*Det2_35_02;
1614 G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04
1615 + m[A24]*Det2_35_03;
1616 G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05
1617 + m[A25]*Det2_35_03;
1618 G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13
1619 + m[A23]*Det2_35_12;
1620 G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14
1621 + m[A24]*Det2_35_12;
1622 G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15
1623 + m[A25]*Det2_35_12;
1624 G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14
1625 + m[A24]*Det2_35_13;
1626 G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15
1627 + m[A25]*Det2_35_13;
1628 G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24
1629 + m[A24]*Det2_35_23;
1630 G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25
1631 + m[A25]*Det2_35_23;
1632 G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02
1633 + m[A22]*Det2_45_01;
1634 G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03
1635 + m[A23]*Det2_45_01;
1636 G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04
1637 + m[A24]*Det2_45_01;
1638 G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05
1639 + m[A25]*Det2_45_01;
1640 G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03
1641 + m[A23]*Det2_45_02;
1642 G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04
1643 + m[A24]*Det2_45_02;
1644 G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05
1645 + m[A25]*Det2_45_02;
1646 G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04
1647 + m[A24]*Det2_45_03;
1648 G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05
1649 + m[A25]*Det2_45_03;
1650 G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05
1651 + m[A25]*Det2_45_04;
1652 G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13
1653 + m[A23]*Det2_45_12;
1654 G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14
1655 + m[A24]*Det2_45_12;
1656 G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15
1657 + m[A25]*Det2_45_12;
1658 G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14
1659 + m[A24]*Det2_45_13;
1660 G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15
1661 + m[A25]*Det2_45_13;
1662 G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15
1663 + m[A25]*Det2_45_14;
1664 G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24
1665 + m[A24]*Det2_45_23;
1666 G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25
1667 + m[A25]*Det2_45_23;
1668 G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25
1669 + m[A25]*Det2_45_24;
1670 G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02
1671 + m[A32]*Det2_45_01;
1672 G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03
1673 + m[A33]*Det2_45_01;
1674 G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04
1675 + m[A34]*Det2_45_01;
1676 G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05
1677 + m[A35]*Det2_45_01;
1678 G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03
1679 + m[A33]*Det2_45_02;
1680 G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04
1681 + m[A34]*Det2_45_02;
1682 G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05
1683 + m[A35]*Det2_45_02;
1684 G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04
1685 + m[A34]*Det2_45_03;
1686 G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05
1687 + m[A35]*Det2_45_03;
1688 G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05
1689 + m[A35]*Det2_45_04;
1690 G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13
1691 + m[A33]*Det2_45_12;
1692 G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14
1693 + m[A34]*Det2_45_12;
1694 G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15
1695 + m[A35]*Det2_45_12;
1696 G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14
1697 + m[A34]*Det2_45_13;
1698 G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15
1699 + m[A35]*Det2_45_13;
1700 G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15
1701 + m[A35]*Det2_45_14;
1702 G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24
1703 + m[A34]*Det2_45_23;
1704 G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25
1705 + m[A35]*Det2_45_23;
1706 G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25
1707 + m[A35]*Det2_45_24;
1708 G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35
1709 + m[A35]*Det2_45_34;
1710
1711 // Find all NECESSARY 4x4 dets: (55 of them)
1712
1713 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023
1714 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
1715 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024
1716 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
1717 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034
1718 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
1719 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034
1720 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
1721 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134
1722 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
1723 G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023
1724 + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
1725 G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024
1726 + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
1727 G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025
1728 + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
1729 G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034
1730 + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
1731 G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035
1732 + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
1733 G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034
1734 + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
1735 G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035
1736 + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
1737 G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134
1738 + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
1739 G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135
1740 + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
1741 G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023
1742 + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
1743 G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024
1744 + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
1745 G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025
1746 + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
1747 G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034
1748 + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
1749 G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035
1750 + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
1751 G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045
1752 + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
1753 G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034
1754 + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
1755 G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035
1756 + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
1757 G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045
1758 + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
1759 G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134
1760 + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
1761 G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135
1762 + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
1763 G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145
1764 + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
1765 G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023
1766 + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
1767 G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024
1768 + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
1769 G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025
1770 + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
1771 G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034
1772 + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
1773 G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035
1774 + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
1775 G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045
1776 + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
1777 G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034
1778 + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
1779 G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035
1780 + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
1781 G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045
1782 + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
1783 G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045
1784 + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
1785 G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134
1786 + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
1787 G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135
1788 + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
1789 G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145
1790 + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
1791 G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145
1792 + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
1793 G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023
1794 + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
1795 G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024
1796 + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
1797 G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025
1798 + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
1799 G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034
1800 + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
1801 G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035
1802 + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
1803 G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045
1804 + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
1805 G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034
1806 + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
1807 G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035
1808 + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
1809 G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045
1810 + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
1811 G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045
1812 + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
1813 G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134
1814 + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
1815 G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135
1816 + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
1817 G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145
1818 + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
1819 G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145
1820 + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
1821 G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245
1822 + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
1823
1824 // Find all NECESSARY 5x5 dets: (19 of them)
1825
1826 G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234
1827 + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
1828 G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234
1829 + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
1830 G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235
1831 + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
1832 G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234
1833 + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
1834 G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235
1835 + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
1836 G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245
1837 + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
1838 G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234
1839 + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
1840 G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235
1841 + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
1842 G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245
1843 + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
1844 G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345
1845 + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
1846 G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234
1847 + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
1848 G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235
1849 + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
1850 G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245
1851 + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
1852 G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345
1853 + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
1854 G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345
1855 + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
1856 G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234
1857 + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
1858 G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235
1859 + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
1860 G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245
1861 + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
1862 G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345
1863 + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
1864 G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345
1865 + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
1866 G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345
1867 + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
1868
1869 // Find the determinant
1870
1871 G4double det = m[A00]*Det5_12345_12345
1872 - m[A01]*Det5_12345_02345
1873 + m[A02]*Det5_12345_01345
1874 - m[A03]*Det5_12345_01245
1875 + m[A04]*Det5_12345_01235
1876 - m[A05]*Det5_12345_01234;
1877
1878 if ( det == 0 )
1879 {
1880 ifail = 1;
1881 return;
1882 }
1883
1884 G4double oneOverDet = 1.0/det;
1885 G4double mn1OverDet = - oneOverDet;
1886
1887 m[A00] = Det5_12345_12345*oneOverDet;
1888 m[A01] = Det5_12345_02345*mn1OverDet;
1889 m[A02] = Det5_12345_01345*oneOverDet;
1890 m[A03] = Det5_12345_01245*mn1OverDet;
1891 m[A04] = Det5_12345_01235*oneOverDet;
1892 m[A05] = Det5_12345_01234*mn1OverDet;
1893
1894 m[A11] = Det5_02345_02345*oneOverDet;
1895 m[A12] = Det5_02345_01345*mn1OverDet;
1896 m[A13] = Det5_02345_01245*oneOverDet;
1897 m[A14] = Det5_02345_01235*mn1OverDet;
1898 m[A15] = Det5_02345_01234*oneOverDet;
1899
1900 m[A22] = Det5_01345_01345*oneOverDet;
1901 m[A23] = Det5_01345_01245*mn1OverDet;
1902 m[A24] = Det5_01345_01235*oneOverDet;
1903 m[A25] = Det5_01345_01234*mn1OverDet;
1904
1905 m[A33] = Det5_01245_01245*oneOverDet;
1906 m[A34] = Det5_01245_01235*mn1OverDet;
1907 m[A35] = Det5_01245_01234*oneOverDet;
1908
1909 m[A44] = Det5_01235_01235*oneOverDet;
1910 m[A45] = Det5_01235_01234*mn1OverDet;
1911
1912 m[A55] = Det5_01234_01234*oneOverDet;
1913
1914 return;
1915}

◆ 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 471 of file G4ErrorSymMatrix.cc.

472{
473 SIMPLE_UOP(*=)
474 return (*this);
475}
SIMPLE_UOP
#define SIMPLE_UOP(OPER)
Definition: G4ErrorMatrix.cc:41

◆ operator+=()

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

Definition at line 426 of file G4ErrorSymMatrix.cc.

427{
428 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=);
429 SIMPLE_BOP(+=)
430 return (*this);
431}
CHK_DIM_2
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition: G4ErrorMatrix.cc:61
SIMPLE_BOP
#define SIMPLE_BOP(OPER)
Definition: G4ErrorMatrix.cc:46
G4ErrorSymMatrix::num_col
G4int num_col() const

◆ operator-()

G4ErrorSymMatrix G4ErrorSymMatrix::operator- ( ) const

Definition at line 196 of file G4ErrorSymMatrix.cc.

197{
198 G4ErrorSymMatrix mat2(nrow);
199 G4ErrorMatrixConstIter a=m.begin();
200 G4ErrorMatrixIter b=mat2.m.begin();
201 G4ErrorMatrixConstIter e=m.begin()+num_size();
202 for(;a<e; a++, b++) { (*b) = -(*a); }
203 return mat2;
204}
G4ErrorSymMatrix::num_size
G4int num_size() const

◆ operator-=()

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

Definition at line 458 of file G4ErrorSymMatrix.cc.

459{
460 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
461 SIMPLE_BOP(-=)
462 return (*this);
463}

◆ operator/=()

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

Definition at line 465 of file G4ErrorSymMatrix.cc.

466{
467 SIMPLE_UOP(/=)
468 return (*this);
469}

◆ operator=() [1/2]

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

Definition at line 508 of file G4ErrorSymMatrix.cc.

509{
510 if (&mat1 == this) { return *this; }
511 if(mat1.nrow != nrow)
512 {
513 nrow = mat1.nrow;
514 size = mat1.size;
515 m.resize(size);
516 }
517 m = mat1.m;
518 return (*this);
519}

◆ 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 588 of file G4ErrorSymMatrix.cc.

589{
590 G4ErrorSymMatrix mret(mat1.num_row());
591 G4ErrorMatrix temp = mat1*(*this);
592
593// If mat1*(*this) has correct dimensions, then so will the mat1.T multiplication.
594// So there is no need to check dimensions again.
595
596 G4int n = mat1.num_col();
597 G4ErrorMatrixIter mr = mret.m.begin();
598 G4ErrorMatrixIter tempr1 = temp.m.begin();
599 for(G4int r=1;r<=mret.num_row();r++)
600 {
601 G4ErrorMatrixConstIter m1c1 = mat1.m.begin();
602 for(G4int c=1;c<=r;c++)
603 {
604 G4double tmp = 0.0;
605 G4ErrorMatrixIter tempri = tempr1;
606 G4ErrorMatrixConstIter m1ci = m1c1;
607 for(G4int i=1;i<=mat1.num_col();i++)
608 {
609 tmp+=(*(tempri++))*(*(m1ci++));
610 }
611 *(mr++) = tmp;
612 m1c1 += n;
613 }
614 tempr1 += n;
615 }
616 return mret;
617}

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

◆ similarity() [2/2]

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

Definition at line 619 of file G4ErrorSymMatrix.cc.

620{
621 G4ErrorSymMatrix mret(mat1.num_row());
622 G4ErrorMatrix temp = mat1*(*this);
623 G4int n = mat1.num_col();
624 G4ErrorMatrixIter mr = mret.m.begin();
625 G4ErrorMatrixIter tempr1 = temp.m.begin();
626 for(G4int r=1;r<=mret.num_row();r++)
627 {
628 G4ErrorMatrixConstIter m1c1 = mat1.m.begin();
629 G4int c;
630 for(c=1;c<=r;c++)
631 {
632 G4double tmp = 0.0;
633 G4ErrorMatrixIter tempri = tempr1;
634 G4ErrorMatrixConstIter m1ci = m1c1;
635 G4int i;
636 for(i=1;i<c;i++)
637 {
638 tmp+=(*(tempri++))*(*(m1ci++));
639 }
640 for(i=c;i<=mat1.num_col();i++)
641 {
642 tmp+=(*(tempri++))*(*(m1ci));
643 m1ci += i;
644 }
645 *(mr++) = tmp;
646 m1c1 += c;
647 }
648 tempr1 += n;
649 }
650 return mret;
651}

◆ similarityT()

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

Definition at line 653 of file G4ErrorSymMatrix.cc.

654{
655 G4ErrorSymMatrix mret(mat1.num_col());
656 G4ErrorMatrix temp = (*this)*mat1;
657 G4int n = mat1.num_col();
658 G4ErrorMatrixIter mrc = mret.m.begin();
659 G4ErrorMatrixIter temp1r = temp.m.begin();
660 for(G4int r=1;r<=mret.num_row();r++)
661 {
662 G4ErrorMatrixConstIter m11c = mat1.m.begin();
663 for(G4int c=1;c<=r;c++)
664 {
665 G4double tmp = 0.0;
666 G4ErrorMatrixIter tempir = temp1r;
667 G4ErrorMatrixConstIter m1ic = m11c;
668 for(G4int i=1;i<=mat1.num_row();i++)
669 {
670 tmp+=(*(tempir))*(*(m1ic));
671 tempir += n;
672 m1ic += n;
673 }
674 *(mrc++) = tmp;
675 m11c++;
676 }
677 temp1r++;
678 }
679 return mret;
680}

◆ sub() [1/3]

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

Definition at line 142 of file G4ErrorSymMatrix.cc.

143{
144 G4ErrorSymMatrix mret(max_row-min_row+1);
145 if(max_row > num_row())
146 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
147 G4ErrorMatrixIter a = mret.m.begin();
148 G4ErrorMatrixIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
149 for(G4int irow=1; irow<=mret.num_row(); irow++)
150 {
151 G4ErrorMatrixIter b = b1;
152 for(G4int icol=1; icol<=irow; icol++)
153 {
154 *(a++) = *(b++);
155 }
156 b1 += irow+min_row-1;
157 }
158 return mret;
159}

◆ sub() [2/3]

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

Definition at line 123 of file G4ErrorSymMatrix.cc.

124{
125 G4ErrorSymMatrix mret(max_row-min_row+1);
126 if(max_row > num_row())
127 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
128 G4ErrorMatrixIter a = mret.m.begin();
129 G4ErrorMatrixConstIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
130 for(G4int irow=1; irow<=mret.num_row(); irow++)
131 {
132 G4ErrorMatrixConstIter b = b1;
133 for(G4int icol=1; icol<=irow; icol++)
134 {
135 *(a++) = *(b++);
136 }
137 b1 += irow+min_row-1;
138 }
139 return mret;
140}

Referenced by dsum().

◆ sub() [3/3]

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

Definition at line 161 of file G4ErrorSymMatrix.cc.

162{
163 if(row <1 || row+mat1.num_row()-1 > num_row() )
164 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
165 G4ErrorMatrixConstIter a = mat1.m.begin();
166 G4ErrorMatrixIter b1 = m.begin() + (row+2)*(row-1)/2;
167 for(G4int irow=1; irow<=mat1.num_row(); irow++)
168 {
169 G4ErrorMatrixIter b = b1;
170 for(G4int icol=1; icol<=irow; icol++)
171 {
172 *(b++) = *(a++);
173 }
174 b1 += irow+row-1;
175 }
176}

◆ T()

G4ErrorSymMatrix G4ErrorSymMatrix::T ( ) const

Referenced by G4ErrorFreeTrajState::PropagateError().

◆ trace()

G4double G4ErrorSymMatrix::trace ( ) const

Definition at line 818 of file G4ErrorSymMatrix.cc.

819{
820 G4double t = 0.0;
821 for (G4int i=0; i<nrow; i++)
822 { t += *(m.begin() + (i+3)*i/2); }
823 return t;
824}

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 195 of file G4ErrorSymMatrix.hh.

◆ G4ErrorSymMatrix_row

friend class G4ErrorSymMatrix_row
friend

Definition at line 193 of file G4ErrorSymMatrix.hh.

◆ G4ErrorSymMatrix_row_const

friend class G4ErrorSymMatrix_row_const
friend

Definition at line 194 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 286 of file G4ErrorSymMatrix.cc.

287{
288 G4ErrorMatrix mret(mat1.num_row(),mat2.num_col());
289 CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
290 G4ErrorMatrixConstIter mit1, mit2, sp,snp; //mit2=0
291 G4double temp;
292 G4ErrorMatrixIter mir=mret.m.begin();
293 for(mit1=mat1.m.begin();
294 mit1<mat1.m.begin()+mat1.num_row()*mat1.num_col();
295 mit1 = mit2)
296 {
297 snp=mat2.m.begin();
298 for(int step=1;step<=mat2.num_row();++step)
299 {
300 mit2=mit1;
301 sp=snp;
302 snp+=step;
303 temp=0;
304 while(sp<snp) // Loop checking, 06.08.2015, G.Cosmo
305 { temp+=*(sp++)*(*(mit2++)); }
306 if( step<mat2.num_row() ) { // only if we aren't on the last row
307 sp+=step-1;
308 for(int stept=step+1;stept<=mat2.num_row();stept++)
309 {
310 temp+=*sp*(*(mit2++));
311 if(stept<mat2.num_row()) sp+=stept;
312 }
313 } // if(step
314 *(mir++)=temp;
315 } // for(step
316 } // for(mit1
317 return mret;
318}
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: G4ErrorMatrix.cc:66

◆ operator* [2/3]

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

Definition at line 320 of file G4ErrorSymMatrix.cc.

321{
322 G4ErrorMatrix mret(mat1.num_row(),mat2.num_col());
323 CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
324 G4int step,stept;
325 G4ErrorMatrixConstIter mit1,mit2,sp,snp;
326 G4double temp;
327 G4ErrorMatrixIter mir=mret.m.begin();
328 for(step=1,snp=mat1.m.begin();step<=mat1.num_row();snp+=step++)
329 {
330 for(mit1=mat2.m.begin();mit1<mat2.m.begin()+mat2.num_col();mit1++)
331 {
332 mit2=mit1;
333 sp=snp;
334 temp=0;
335 while(sp<snp+step) // Loop checking, 06.08.2015, G.Cosmo
336 {
337 temp+=*mit2*(*(sp++));
338 mit2+=mat2.num_col();
339 }
340 sp+=step-1;
341 for(stept=step+1;stept<=mat1.num_row();stept++)
342 {
343 temp+=*mit2*(*sp);
344 mit2+=mat2.num_col();
345 sp+=stept;
346 }
347 *(mir++)=temp;
348 }
349 }
350 return mret;
351}

◆ operator* [3/3]

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

Definition at line 353 of file G4ErrorSymMatrix.cc.

354{
355 G4ErrorMatrix mret(mat1.num_row(),mat1.num_row());
356 CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
357 G4int step1,stept1,step2,stept2;
358 G4ErrorMatrixConstIter snp1,sp1,snp2,sp2;
359 G4double temp;
360 G4ErrorMatrixIter mr = mret.m.begin();
361 for(step1=1,snp1=mat1.m.begin();step1<=mat1.num_row();snp1+=step1++)
362 {
363 for(step2=1,snp2=mat2.m.begin();step2<=mat2.num_row();)
364 {
365 sp1=snp1;
366 sp2=snp2;
367 snp2+=step2;
368 temp=0;
369 if(step1<step2)
370 {
371 while(sp1<snp1+step1) // Loop checking, 06.08.2015, G.Cosmo
372 { temp+=(*(sp1++))*(*(sp2++)); }
373 sp1+=step1-1;
374 for(stept1=step1+1;stept1!=step2+1;sp1+=stept1++)
375 { temp+=(*sp1)*(*(sp2++)); }
376 sp2+=step2-1;
377 for(stept2=++step2;stept2<=mat2.num_row();sp1+=stept1++,sp2+=stept2++)
378 { temp+=(*sp1)*(*sp2); }
379 }
380 else
381 {
382 while(sp2<snp2) // Loop checking, 06.08.2015, G.Cosmo
383 { temp+=(*(sp1++))*(*(sp2++)); }
384 sp2+=step2-1;
385 for(stept2=++step2;stept2!=step1+1;sp2+=stept2++)
386 { temp+=(*(sp1++))*(*sp2); }
387 sp1+=step1-1;
388 for(stept1=step1+1;stept1<=mat1.num_row();sp1+=stept1++,sp2+=stept2++)
389 { temp+=(*sp1)*(*sp2); }
390 }
391 *(mr++)=temp;
392 }
393 }
394 return mret;
395}

◆ operator+

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

Definition at line 222 of file G4ErrorSymMatrix.cc.

224{
225 G4ErrorSymMatrix mret(mat1.nrow);
226 CHK_DIM_1(mat1.nrow, mat2.nrow,+);
227 SIMPLE_TOP(+)
228 return mret;
229}
SIMPLE_TOP
#define SIMPLE_TOP(OPER)
Definition: G4ErrorMatrix.cc:52

◆ operator-

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

Definition at line 251 of file G4ErrorSymMatrix.cc.

253{
254 G4ErrorSymMatrix mret(mat1.num_row());
255 CHK_DIM_1(mat1.num_row(),mat2.num_row(),-);
256 SIMPLE_TOP(-)
257 return mret;
258}

◆ tridiagonal

void tridiagonal ( G4ErrorSymMatrix *  a,
G4ErrorMatrix *  hsm 
)
friend
The documentation for this class was generated from the following files: