#include <G4ErrorSymMatrix.hh>
Definition at line 43 of file G4ErrorSymMatrix.hh.
G4ErrorSymMatrix::G4ErrorSymMatrix | ( | ) | [inline] |
G4ErrorSymMatrix::G4ErrorSymMatrix | ( | G4int | p | ) | [explicit] |
Definition at line 71 of file G4ErrorSymMatrix.cc.
00072 : m(p*(p+1)/2), nrow(p) 00073 { 00074 size = nrow * (nrow+1) / 2; 00075 m.assign(size,0); 00076 }
Definition at line 78 of file G4ErrorSymMatrix.cc.
References G4ErrorMatrix::error().
00079 : m(p*(p+1)/2), nrow(p) 00080 { 00081 size = nrow * (nrow+1) / 2; 00082 00083 m.assign(size,0); 00084 switch(init) 00085 { 00086 case 0: 00087 break; 00088 00089 case 1: 00090 { 00091 G4ErrorMatrixIter a = m.begin(); 00092 for(G4int i=1;i<=nrow;i++) 00093 { 00094 *a = 1.0; 00095 a += (i+1); 00096 } 00097 break; 00098 } 00099 default: 00100 G4ErrorMatrix::error("G4ErrorSymMatrix: initialization must be 0 or 1."); 00101 } 00102 }
G4ErrorSymMatrix::G4ErrorSymMatrix | ( | const G4ErrorSymMatrix & | m1 | ) |
G4ErrorSymMatrix::~G4ErrorSymMatrix | ( | ) | [virtual] |
G4ErrorSymMatrix G4ErrorSymMatrix::apply | ( | G4double(*)(G4double, G4int, G4int) | f | ) | const |
Definition at line 553 of file G4ErrorSymMatrix.cc.
00554 { 00555 G4ErrorSymMatrix mret(num_row()); 00556 G4ErrorMatrixConstIter a = m.begin(); 00557 G4ErrorMatrixIter b = mret.m.begin(); 00558 for(G4int ir=1;ir<=num_row();ir++) 00559 { 00560 for(G4int ic=1;ic<=ir;ic++) 00561 { 00562 *(b++) = (*f)(*(a++), ir, ic); 00563 } 00564 } 00565 return mret; 00566 }
void G4ErrorSymMatrix::assign | ( | const G4ErrorSymMatrix & | m2 | ) | [inline] |
void G4ErrorSymMatrix::assign | ( | const G4ErrorMatrix & | m2 | ) |
Definition at line 568 of file G4ErrorSymMatrix.cc.
References G4ErrorMatrix::m, and G4ErrorMatrix::nrow.
00569 { 00570 if(mat1.nrow != nrow) 00571 { 00572 nrow = mat1.nrow; 00573 size = nrow * (nrow+1) / 2; 00574 m.resize(size); 00575 } 00576 G4ErrorMatrixConstIter a = mat1.m.begin(); 00577 G4ErrorMatrixIter b = m.begin(); 00578 for(G4int r=1;r<=nrow;r++) 00579 { 00580 G4ErrorMatrixConstIter d = a; 00581 for(G4int c=1;c<=r;c++) 00582 { 00583 *(b++) = *(d++); 00584 } 00585 a += nrow; 00586 } 00587 }
G4double G4ErrorSymMatrix::determinant | ( | ) | const |
Definition at line 799 of file G4ErrorSymMatrix.cc.
References G4ErrorMatrix::dfact_matrix().
00800 { 00801 static const G4int max_array = 20; 00802 00803 // ir must point to an array which is ***1 longer than*** nrow 00804 00805 static std::vector<G4int> ir_vec (max_array+1); 00806 if (ir_vec.size() <= static_cast<unsigned int>(nrow)) 00807 { 00808 ir_vec.resize(nrow+1); 00809 } 00810 G4int * ir = &ir_vec[0]; 00811 00812 G4double det; 00813 G4ErrorMatrix mt(*this); 00814 G4int i = mt.dfact_matrix(det, ir); 00815 if(i==0) { return det; } 00816 return 0.0; 00817 }
G4ErrorSymMatrix G4ErrorSymMatrix::inverse | ( | G4int & | ifail | ) | const [inline] |
Definition at line 135 of file G4ErrorSymMatrix.icc.
References invert().
00136 { 00137 G4ErrorSymMatrix mTmp(*this); 00138 mTmp.invert(ifail); 00139 return mTmp; 00140 }
void G4ErrorSymMatrix::invert | ( | G4int & | ifail | ) |
Definition at line 683 of file G4ErrorSymMatrix.cc.
References invertBunchKaufman().
Referenced by inverse().
00684 { 00685 ifail = 0; 00686 00687 switch(nrow) 00688 { 00689 case 3: 00690 { 00691 G4double det, temp; 00692 G4double t1, t2, t3; 00693 G4double c11,c12,c13,c22,c23,c33; 00694 c11 = (*(m.begin()+2)) * (*(m.begin()+5)) 00695 - (*(m.begin()+4)) * (*(m.begin()+4)); 00696 c12 = (*(m.begin()+4)) * (*(m.begin()+3)) 00697 - (*(m.begin()+1)) * (*(m.begin()+5)); 00698 c13 = (*(m.begin()+1)) * (*(m.begin()+4)) 00699 - (*(m.begin()+2)) * (*(m.begin()+3)); 00700 c22 = (*(m.begin()+5)) * (*m.begin()) 00701 - (*(m.begin()+3)) * (*(m.begin()+3)); 00702 c23 = (*(m.begin()+3)) * (*(m.begin()+1)) 00703 - (*(m.begin()+4)) * (*m.begin()); 00704 c33 = (*m.begin()) * (*(m.begin()+2)) 00705 - (*(m.begin()+1)) * (*(m.begin()+1)); 00706 t1 = std::fabs(*m.begin()); 00707 t2 = std::fabs(*(m.begin()+1)); 00708 t3 = std::fabs(*(m.begin()+3)); 00709 if (t1 >= t2) 00710 { 00711 if (t3 >= t1) 00712 { 00713 temp = *(m.begin()+3); 00714 det = c23*c12-c22*c13; 00715 } 00716 else 00717 { 00718 temp = *m.begin(); 00719 det = c22*c33-c23*c23; 00720 } 00721 } 00722 else if (t3 >= t2) 00723 { 00724 temp = *(m.begin()+3); 00725 det = c23*c12-c22*c13; 00726 } 00727 else 00728 { 00729 temp = *(m.begin()+1); 00730 det = c13*c23-c12*c33; 00731 } 00732 if (det==0) 00733 { 00734 ifail = 1; 00735 return; 00736 } 00737 { 00738 G4double ss = temp/det; 00739 G4ErrorMatrixIter mq = m.begin(); 00740 *(mq++) = ss*c11; 00741 *(mq++) = ss*c12; 00742 *(mq++) = ss*c22; 00743 *(mq++) = ss*c13; 00744 *(mq++) = ss*c23; 00745 *(mq) = ss*c33; 00746 } 00747 } 00748 break; 00749 case 2: 00750 { 00751 G4double det, temp, ss; 00752 det = (*m.begin())*(*(m.begin()+2)) - (*(m.begin()+1))*(*(m.begin()+1)); 00753 if (det==0) 00754 { 00755 ifail = 1; 00756 return; 00757 } 00758 ss = 1.0/det; 00759 *(m.begin()+1) *= -ss; 00760 temp = ss*(*(m.begin()+2)); 00761 *(m.begin()+2) = ss*(*m.begin()); 00762 *m.begin() = temp; 00763 break; 00764 } 00765 case 1: 00766 { 00767 if ((*m.begin())==0) 00768 { 00769 ifail = 1; 00770 return; 00771 } 00772 *m.begin() = 1.0/(*m.begin()); 00773 break; 00774 } 00775 case 5: 00776 { 00777 invert5(ifail); 00778 return; 00779 } 00780 case 6: 00781 { 00782 invert6(ifail); 00783 return; 00784 } 00785 case 4: 00786 { 00787 invert4(ifail); 00788 return; 00789 } 00790 default: 00791 { 00792 invertBunchKaufman(ifail); 00793 return; 00794 } 00795 } 00796 return; // inversion successful 00797 }
void G4ErrorSymMatrix::invertBunchKaufman | ( | G4int & | ifail | ) |
Definition at line 827 of file G4ErrorSymMatrix.cc.
References G4InuclParticleNames::alpha, DBL_EPSILON, G4cerr, G4endl, and G4InuclParticleNames::lambda.
Referenced by invert().
00828 { 00829 // Bunch-Kaufman diagonal pivoting method 00830 // It is decribed in J.R. Bunch, L. Kaufman (1977). 00831 // "Some Stable Methods for Calculating Inertia and Solving Symmetric 00832 // Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub, 00833 // Charles F. van Loan, "Matrix Computations" (the second edition 00834 // has a bug.) and implemented in "lapack" 00835 // Mario Stanke, 09/97 00836 00837 G4int i, j, k, ss; 00838 G4int pivrow; 00839 00840 // Establish the two working-space arrays needed: x and piv are 00841 // used as pointers to arrays of doubles and ints respectively, each 00842 // of length nrow. We do not want to reallocate each time through 00843 // unless the size needs to grow. We do not want to leak memory, even 00844 // by having a new without a delete that is only done once. 00845 00846 static const G4int max_array = 25; 00847 static std::vector<G4double> xvec (max_array); 00848 static std::vector<G4int> pivv (max_array); 00849 typedef std::vector<G4int>::iterator pivIter; 00850 if (xvec.size() < static_cast<unsigned int>(nrow)) xvec.resize(nrow); 00851 if (pivv.size() < static_cast<unsigned int>(nrow)) pivv.resize(nrow); 00852 // Note - resize should do nothing if the size is already larger than nrow, 00853 // but on VC++ there are indications that it does so we check. 00854 // Note - the data elements in a vector are guaranteed to be contiguous, 00855 // so x[i] and piv[i] are optimally fast. 00856 G4ErrorMatrixIter x = xvec.begin(); 00857 // x[i] is used as helper storage, needs to have at least size nrow. 00858 pivIter piv = pivv.begin(); 00859 // piv[i] is used to store details of exchanges 00860 00861 G4double temp1, temp2; 00862 G4ErrorMatrixIter ip, mjj, iq; 00863 G4double lambda, sigma; 00864 const G4double alpha = .6404; // = (1+sqrt(17))/8 00865 const G4double epsilon = 32*DBL_EPSILON; 00866 // whenever a sum of two doubles is below or equal to epsilon 00867 // it is set to zero. 00868 // this constant could be set to zero but then the algorithm 00869 // doesn't neccessarily detect that a matrix is singular 00870 00871 for (i = 0; i < nrow; i++) 00872 { 00873 piv[i] = i+1; 00874 } 00875 00876 ifail = 0; 00877 00878 // compute the factorization P*A*P^T = L * D * L^T 00879 // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices 00880 // L and D^-1 are stored in A = *this, P is stored in piv[] 00881 00882 for (j=1; j < nrow; j+=ss) // main loop over columns 00883 { 00884 mjj = m.begin() + j*(j-1)/2 + j-1; 00885 lambda = 0; // compute lambda = max of A(j+1:n,j) 00886 pivrow = j+1; 00887 ip = m.begin() + (j+1)*j/2 + j-1; 00888 for (i=j+1; i <= nrow ; ip += i++) 00889 { 00890 if (std::fabs(*ip) > lambda) 00891 { 00892 lambda = std::fabs(*ip); 00893 pivrow = i; 00894 } 00895 } 00896 if (lambda == 0 ) 00897 { 00898 if (*mjj == 0) 00899 { 00900 ifail = 1; 00901 return; 00902 } 00903 ss=1; 00904 *mjj = 1./ *mjj; 00905 } 00906 else 00907 { 00908 if (std::fabs(*mjj) >= lambda*alpha) 00909 { 00910 ss=1; 00911 pivrow=j; 00912 } 00913 else 00914 { 00915 sigma = 0; // compute sigma = max A(pivrow, j:pivrow-1) 00916 ip = m.begin() + pivrow*(pivrow-1)/2+j-1; 00917 for (k=j; k < pivrow; k++) 00918 { 00919 if (std::fabs(*ip) > sigma) 00920 sigma = std::fabs(*ip); 00921 ip++; 00922 } 00923 if (sigma * std::fabs(*mjj) >= alpha * lambda * lambda) 00924 { 00925 ss=1; 00926 pivrow = j; 00927 } 00928 else if (std::fabs(*(m.begin()+pivrow*(pivrow-1)/2+pivrow-1)) 00929 >= alpha * sigma) 00930 { ss=1; } 00931 else 00932 { ss=2; } 00933 } 00934 if (pivrow == j) // no permutation neccessary 00935 { 00936 piv[j-1] = pivrow; 00937 if (*mjj == 0) 00938 { 00939 ifail=1; 00940 return; 00941 } 00942 temp2 = *mjj = 1./ *mjj; // invert D(j,j) 00943 00944 // update A(j+1:n, j+1,n) 00945 for (i=j+1; i <= nrow; i++) 00946 { 00947 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2; 00948 ip = m.begin()+i*(i-1)/2+j; 00949 for (k=j+1; k<=i; k++) 00950 { 00951 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1); 00952 if (std::fabs(*ip) <= epsilon) 00953 { *ip=0; } 00954 ip++; 00955 } 00956 } 00957 // update L 00958 ip = m.begin() + (j+1)*j/2 + j-1; 00959 for (i=j+1; i <= nrow; ip += i++) 00960 { 00961 *ip *= temp2; 00962 } 00963 } 00964 else if (ss==1) // 1x1 pivot 00965 { 00966 piv[j-1] = pivrow; 00967 00968 // interchange rows and columns j and pivrow in 00969 // submatrix (j:n,j:n) 00970 ip = m.begin() + pivrow*(pivrow-1)/2 + j; 00971 for (i=j+1; i < pivrow; i++, ip++) 00972 { 00973 temp1 = *(m.begin() + i*(i-1)/2 + j-1); 00974 *(m.begin() + i*(i-1)/2 + j-1)= *ip; 00975 *ip = temp1; 00976 } 00977 temp1 = *mjj; 00978 *mjj = *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1); 00979 *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1) = temp1; 00980 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1; 00981 iq = ip + pivrow-j; 00982 for (i = pivrow+1; i <= nrow; ip += i, iq += i++) 00983 { 00984 temp1 = *iq; 00985 *iq = *ip; 00986 *ip = temp1; 00987 } 00988 00989 if (*mjj == 0) 00990 { 00991 ifail = 1; 00992 return; 00993 } 00994 temp2 = *mjj = 1./ *mjj; // invert D(j,j) 00995 00996 // update A(j+1:n, j+1:n) 00997 for (i = j+1; i <= nrow; i++) 00998 { 00999 temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2; 01000 ip = m.begin()+i*(i-1)/2+j; 01001 for (k=j+1; k<=i; k++) 01002 { 01003 *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1); 01004 if (std::fabs(*ip) <= epsilon) 01005 { *ip=0; } 01006 ip++; 01007 } 01008 } 01009 // update L 01010 ip = m.begin() + (j+1)*j/2 + j-1; 01011 for (i=j+1; i<=nrow; ip += i++) 01012 { 01013 *ip *= temp2; 01014 } 01015 } 01016 else // ss=2, ie use a 2x2 pivot 01017 { 01018 piv[j-1] = -pivrow; 01019 piv[j] = 0; // that means this is the second row of a 2x2 pivot 01020 01021 if (j+1 != pivrow) 01022 { 01023 // interchange rows and columns j+1 and pivrow in 01024 // submatrix (j:n,j:n) 01025 ip = m.begin() + pivrow*(pivrow-1)/2 + j+1; 01026 for (i=j+2; i < pivrow; i++, ip++) 01027 { 01028 temp1 = *(m.begin() + i*(i-1)/2 + j); 01029 *(m.begin() + i*(i-1)/2 + j) = *ip; 01030 *ip = temp1; 01031 } 01032 temp1 = *(mjj + j + 1); 01033 *(mjj + j + 1) = 01034 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1); 01035 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1; 01036 temp1 = *(mjj + j); 01037 *(mjj + j) = *(m.begin() + pivrow*(pivrow-1)/2 + j-1); 01038 *(m.begin() + pivrow*(pivrow-1)/2 + j-1) = temp1; 01039 ip = m.begin() + (pivrow+1)*pivrow/2 + j; 01040 iq = ip + pivrow-(j+1); 01041 for (i = pivrow+1; i <= nrow; ip += i, iq += i++) 01042 { 01043 temp1 = *iq; 01044 *iq = *ip; 01045 *ip = temp1; 01046 } 01047 } 01048 // invert D(j:j+1,j:j+1) 01049 temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j); 01050 if (temp2 == 0) 01051 { 01052 G4cerr 01053 << "G4ErrorSymMatrix::bunch_invert: error in pivot choice" 01054 << G4endl; 01055 } 01056 temp2 = 1. / temp2; 01057 01058 // this quotient is guaranteed to exist by the choice 01059 // of the pivot 01060 01061 temp1 = *mjj; 01062 *mjj = *(mjj + j + 1) * temp2; 01063 *(mjj + j + 1) = temp1 * temp2; 01064 *(mjj + j) = - *(mjj + j) * temp2; 01065 01066 if (j < nrow-1) // otherwise do nothing 01067 { 01068 // update A(j+2:n, j+2:n) 01069 for (i=j+2; i <= nrow ; i++) 01070 { 01071 ip = m.begin() + i*(i-1)/2 + j-1; 01072 temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j); 01073 if (std::fabs(temp1 ) <= epsilon) 01074 { temp1 = 0; } 01075 temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1); 01076 if (std::fabs(temp2 ) <= epsilon) 01077 { temp2 = 0; } 01078 for (k = j+2; k <= i ; k++) 01079 { 01080 ip = m.begin() + i*(i-1)/2 + k-1; 01081 iq = m.begin() + k*(k-1)/2 + j-1; 01082 *ip -= temp1 * *iq + temp2 * *(iq+1); 01083 if (std::fabs(*ip) <= epsilon) 01084 { *ip = 0; } 01085 } 01086 } 01087 // update L 01088 for (i=j+2; i <= nrow ; i++) 01089 { 01090 ip = m.begin() + i*(i-1)/2 + j-1; 01091 temp1 = *ip * *mjj + *(ip+1) * *(mjj + j); 01092 if (std::fabs(temp1) <= epsilon) 01093 { temp1 = 0; } 01094 *(ip+1) = *ip * *(mjj + j) + *(ip+1) * *(mjj + j + 1); 01095 if (std::fabs(*(ip+1)) <= epsilon) 01096 { *(ip+1) = 0; } 01097 *ip = temp1; 01098 } 01099 } 01100 } 01101 } 01102 } // end of main loop over columns 01103 01104 if (j == nrow) // the the last pivot is 1x1 01105 { 01106 mjj = m.begin() + j*(j-1)/2 + j-1; 01107 if (*mjj == 0) 01108 { 01109 ifail = 1; 01110 return; 01111 } 01112 else 01113 { 01114 *mjj = 1. / *mjj; 01115 } 01116 } // end of last pivot code 01117 01118 // computing the inverse from the factorization 01119 01120 for (j = nrow ; j >= 1 ; j -= ss) // loop over columns 01121 { 01122 mjj = m.begin() + j*(j-1)/2 + j-1; 01123 if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse 01124 { 01125 ss = 1; 01126 if (j < nrow) 01127 { 01128 ip = m.begin() + (j+1)*j/2 + j-1; 01129 for (i=0; i < nrow-j; ip += 1+j+i++) 01130 { 01131 x[i] = *ip; 01132 } 01133 for (i=j+1; i<=nrow ; i++) 01134 { 01135 temp2=0; 01136 ip = m.begin() + i*(i-1)/2 + j; 01137 for (k=0; k <= i-j-1; k++) 01138 { temp2 += *ip++ * x[k]; } 01139 for (ip += i-1; k < nrow-j; ip += 1+j+k++) 01140 { temp2 += *ip * x[k]; } 01141 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2; 01142 } 01143 temp2 = 0; 01144 ip = m.begin() + (j+1)*j/2 + j-1; 01145 for (k=0; k < nrow-j; ip += 1+j+k++) 01146 { temp2 += x[k] * *ip; } 01147 *mjj -= temp2; 01148 } 01149 } 01150 else //2x2 pivot, compute columns j and j-1 of the inverse 01151 { 01152 if (piv[j-1] != 0) 01153 { G4cerr << "error in piv" << piv[j-1] << G4endl; } 01154 ss=2; 01155 if (j < nrow) 01156 { 01157 ip = m.begin() + (j+1)*j/2 + j-1; 01158 for (i=0; i < nrow-j; ip += 1+j+i++) 01159 { 01160 x[i] = *ip; 01161 } 01162 for (i=j+1; i<=nrow ; i++) 01163 { 01164 temp2 = 0; 01165 ip = m.begin() + i*(i-1)/2 + j; 01166 for (k=0; k <= i-j-1; k++) 01167 { temp2 += *ip++ * x[k]; } 01168 for (ip += i-1; k < nrow-j; ip += 1+j+k++) 01169 { temp2 += *ip * x[k]; } 01170 *(m.begin()+ i*(i-1)/2 + j-1) = -temp2; 01171 } 01172 temp2 = 0; 01173 ip = m.begin() + (j+1)*j/2 + j-1; 01174 for (k=0; k < nrow-j; ip += 1+j+k++) 01175 { temp2 += x[k] * *ip; } 01176 *mjj -= temp2; 01177 temp2 = 0; 01178 ip = m.begin() + (j+1)*j/2 + j-2; 01179 for (i=j+1; i <= nrow; ip += i++) 01180 { temp2 += *ip * *(ip+1); } 01181 *(mjj-1) -= temp2; 01182 ip = m.begin() + (j+1)*j/2 + j-2; 01183 for (i=0; i < nrow-j; ip += 1+j+i++) 01184 { 01185 x[i] = *ip; 01186 } 01187 for (i=j+1; i <= nrow ; i++) 01188 { 01189 temp2 = 0; 01190 ip = m.begin() + i*(i-1)/2 + j; 01191 for (k=0; k <= i-j-1; k++) 01192 { temp2 += *ip++ * x[k]; } 01193 for (ip += i-1; k < nrow-j; ip += 1+j+k++) 01194 { temp2 += *ip * x[k]; } 01195 *(m.begin()+ i*(i-1)/2 + j-2)= -temp2; 01196 } 01197 temp2 = 0; 01198 ip = m.begin() + (j+1)*j/2 + j-2; 01199 for (k=0; k < nrow-j; ip += 1+j+k++) 01200 { temp2 += x[k] * *ip; } 01201 *(mjj-j) -= temp2; 01202 } 01203 } 01204 01205 // interchange rows and columns j and piv[j-1] 01206 // or rows and columns j and -piv[j-2] 01207 01208 pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1]; 01209 ip = m.begin() + pivrow*(pivrow-1)/2 + j; 01210 for (i=j+1;i < pivrow; i++, ip++) 01211 { 01212 temp1 = *(m.begin() + i*(i-1)/2 + j-1); 01213 *(m.begin() + i*(i-1)/2 + j-1) = *ip; 01214 *ip = temp1; 01215 } 01216 temp1 = *mjj; 01217 *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1); 01218 *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1; 01219 if (ss==2) 01220 { 01221 temp1 = *(mjj-1); 01222 *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2); 01223 *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1; 01224 } 01225 01226 ip = m.begin() + (pivrow+1)*pivrow/2 + j-1; // &A(i,j) 01227 iq = ip + pivrow-j; 01228 for (i = pivrow+1; i <= nrow; ip += i, iq += i++) 01229 { 01230 temp1 = *iq; 01231 *iq = *ip; 01232 *ip = temp1; 01233 } 01234 } // end of loop over columns (in computing inverse from factorization) 01235 01236 return; // inversion successful 01237 }
void G4ErrorSymMatrix::invertCholesky5 | ( | G4int & | ifail | ) |
Definition at line 1911 of file G4ErrorSymMatrix.cc.
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, and A44.
01912 { 01913 01914 // Invert by 01915 // 01916 // a) decomposing M = G*G^T with G lower triangular 01917 // (if M is not positive definite this will fail, leaving this unchanged) 01918 // b) inverting G to form H 01919 // c) multiplying H^T * H to get M^-1. 01920 // 01921 // If the matrix is pos. def. it is inverted and 1 is returned. 01922 // If the matrix is not pos. def. it remains unaltered and 0 is returned. 01923 01924 G4double h10; // below-diagonal elements of H 01925 G4double h20, h21; 01926 G4double h30, h31, h32; 01927 G4double h40, h41, h42, h43; 01928 01929 G4double h00, h11, h22, h33, h44; // 1/diagonal elements of G = 01930 // diagonal elements of H 01931 01932 G4double g10; // below-diagonal elements of G 01933 G4double g20, g21; 01934 G4double g30, g31, g32; 01935 G4double g40, g41, g42, g43; 01936 01937 ifail = 1; // We start by assuing we won't succeed... 01938 01939 // Form G -- compute diagonal members of H directly rather than of G 01940 //------- 01941 01942 // Scale first column by 1/sqrt(A00) 01943 01944 h00 = m[A00]; 01945 if (h00 <= 0) { return; } 01946 h00 = 1.0 / std::sqrt(h00); 01947 01948 g10 = m[A10] * h00; 01949 g20 = m[A20] * h00; 01950 g30 = m[A30] * h00; 01951 g40 = m[A40] * h00; 01952 01953 // Form G11 (actually, just h11) 01954 01955 h11 = m[A11] - (g10 * g10); 01956 if (h11 <= 0) { return; } 01957 h11 = 1.0 / std::sqrt(h11); 01958 01959 // Subtract inter-column column dot products from rest of column 1 and 01960 // scale to get column 1 of G 01961 01962 g21 = (m[A21] - (g10 * g20)) * h11; 01963 g31 = (m[A31] - (g10 * g30)) * h11; 01964 g41 = (m[A41] - (g10 * g40)) * h11; 01965 01966 // Form G22 (actually, just h22) 01967 01968 h22 = m[A22] - (g20 * g20) - (g21 * g21); 01969 if (h22 <= 0) { return; } 01970 h22 = 1.0 / std::sqrt(h22); 01971 01972 // Subtract inter-column column dot products from rest of column 2 and 01973 // scale to get column 2 of G 01974 01975 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22; 01976 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22; 01977 01978 // Form G33 (actually, just h33) 01979 01980 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32); 01981 if (h33 <= 0) { return; } 01982 h33 = 1.0 / std::sqrt(h33); 01983 01984 // Subtract inter-column column dot product from A43 and scale to get G43 01985 01986 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33; 01987 01988 // Finally form h44 - if this is possible inversion succeeds 01989 01990 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43); 01991 if (h44 <= 0) { return; } 01992 h44 = 1.0 / std::sqrt(h44); 01993 01994 // Form H = 1/G -- diagonal members of H are already correct 01995 //------------- 01996 01997 // The order here is dictated by speed considerations 01998 01999 h43 = -h33 * g43 * h44; 02000 h32 = -h22 * g32 * h33; 02001 h42 = -h22 * (g32 * h43 + g42 * h44); 02002 h21 = -h11 * g21 * h22; 02003 h31 = -h11 * (g21 * h32 + g31 * h33); 02004 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44); 02005 h10 = -h00 * g10 * h11; 02006 h20 = -h00 * (g10 * h21 + g20 * h22); 02007 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33); 02008 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44); 02009 02010 // Change this to its inverse = H^T*H 02011 //------------------------------------ 02012 02013 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40; 02014 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41; 02015 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41; 02016 m[A02] = h20 * h22 + h30 * h32 + h40 * h42; 02017 m[A12] = h21 * h22 + h31 * h32 + h41 * h42; 02018 m[A22] = h22 * h22 + h32 * h32 + h42 * h42; 02019 m[A03] = h30 * h33 + h40 * h43; 02020 m[A13] = h31 * h33 + h41 * h43; 02021 m[A23] = h32 * h33 + h42 * h43; 02022 m[A33] = h33 * h33 + h43 * h43; 02023 m[A04] = h40 * h44; 02024 m[A14] = h41 * h44; 02025 m[A24] = h42 * h44; 02026 m[A34] = h43 * h44; 02027 m[A44] = h44 * h44; 02028 02029 ifail = 0; 02030 return; 02031 }
void G4ErrorSymMatrix::invertCholesky6 | ( | G4int & | ifail | ) |
Definition at line 2033 of file G4ErrorSymMatrix.cc.
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, and A55.
02034 { 02035 // Invert by 02036 // 02037 // a) decomposing M = G*G^T with G lower triangular 02038 // (if M is not positive definite this will fail, leaving this unchanged) 02039 // b) inverting G to form H 02040 // c) multiplying H^T * H to get M^-1. 02041 // 02042 // If the matrix is pos. def. it is inverted and 1 is returned. 02043 // If the matrix is not pos. def. it remains unaltered and 0 is returned. 02044 02045 G4double h10; // below-diagonal elements of H 02046 G4double h20, h21; 02047 G4double h30, h31, h32; 02048 G4double h40, h41, h42, h43; 02049 G4double h50, h51, h52, h53, h54; 02050 02051 G4double h00, h11, h22, h33, h44, h55; // 1/diagonal elements of G = 02052 // diagonal elements of H 02053 02054 G4double g10; // below-diagonal elements of G 02055 G4double g20, g21; 02056 G4double g30, g31, g32; 02057 G4double g40, g41, g42, g43; 02058 G4double g50, g51, g52, g53, g54; 02059 02060 ifail = 1; // We start by assuing we won't succeed... 02061 02062 // Form G -- compute diagonal members of H directly rather than of G 02063 //------- 02064 02065 // Scale first column by 1/sqrt(A00) 02066 02067 h00 = m[A00]; 02068 if (h00 <= 0) { return; } 02069 h00 = 1.0 / std::sqrt(h00); 02070 02071 g10 = m[A10] * h00; 02072 g20 = m[A20] * h00; 02073 g30 = m[A30] * h00; 02074 g40 = m[A40] * h00; 02075 g50 = m[A50] * h00; 02076 02077 // Form G11 (actually, just h11) 02078 02079 h11 = m[A11] - (g10 * g10); 02080 if (h11 <= 0) { return; } 02081 h11 = 1.0 / std::sqrt(h11); 02082 02083 // Subtract inter-column column dot products from rest of column 1 and 02084 // scale to get column 1 of G 02085 02086 g21 = (m[A21] - (g10 * g20)) * h11; 02087 g31 = (m[A31] - (g10 * g30)) * h11; 02088 g41 = (m[A41] - (g10 * g40)) * h11; 02089 g51 = (m[A51] - (g10 * g50)) * h11; 02090 02091 // Form G22 (actually, just h22) 02092 02093 h22 = m[A22] - (g20 * g20) - (g21 * g21); 02094 if (h22 <= 0) { return; } 02095 h22 = 1.0 / std::sqrt(h22); 02096 02097 // Subtract inter-column column dot products from rest of column 2 and 02098 // scale to get column 2 of G 02099 02100 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22; 02101 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22; 02102 g52 = (m[A52] - (g20 * g50) - (g21 * g51)) * h22; 02103 02104 // Form G33 (actually, just h33) 02105 02106 h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32); 02107 if (h33 <= 0) { return; } 02108 h33 = 1.0 / std::sqrt(h33); 02109 02110 // Subtract inter-column column dot products from rest of column 3 and 02111 // scale to get column 3 of G 02112 02113 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33; 02114 g53 = (m[A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33; 02115 02116 // Form G44 (actually, just h44) 02117 02118 h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43); 02119 if (h44 <= 0) { return; } 02120 h44 = 1.0 / std::sqrt(h44); 02121 02122 // Subtract inter-column column dot product from M54 and scale to get G54 02123 02124 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44; 02125 02126 // Finally form h55 - if this is possible inversion succeeds 02127 02128 h55 = m[A55] - (g50*g50) - (g51*g51) - (g52*g52) - (g53*g53) - (g54*g54); 02129 if (h55 <= 0) { return; } 02130 h55 = 1.0 / std::sqrt(h55); 02131 02132 // Form H = 1/G -- diagonal members of H are already correct 02133 //------------- 02134 02135 // The order here is dictated by speed considerations 02136 02137 h54 = -h44 * g54 * h55; 02138 h43 = -h33 * g43 * h44; 02139 h53 = -h33 * (g43 * h54 + g53 * h55); 02140 h32 = -h22 * g32 * h33; 02141 h42 = -h22 * (g32 * h43 + g42 * h44); 02142 h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * h55); 02143 h21 = -h11 * g21 * h22; 02144 h31 = -h11 * (g21 * h32 + g31 * h33); 02145 h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44); 02146 h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55); 02147 h10 = -h00 * g10 * h11; 02148 h20 = -h00 * (g10 * h21 + g20 * h22); 02149 h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33); 02150 h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44); 02151 h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55); 02152 02153 // Change this to its inverse = H^T*H 02154 //------------------------------------ 02155 02156 m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40 + h50*h50; 02157 m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51; 02158 m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51; 02159 m[A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52; 02160 m[A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52; 02161 m[A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52; 02162 m[A03] = h30 * h33 + h40 * h43 + h50 * h53; 02163 m[A13] = h31 * h33 + h41 * h43 + h51 * h53; 02164 m[A23] = h32 * h33 + h42 * h43 + h52 * h53; 02165 m[A33] = h33 * h33 + h43 * h43 + h53 * h53; 02166 m[A04] = h40 * h44 + h50 * h54; 02167 m[A14] = h41 * h44 + h51 * h54; 02168 m[A24] = h42 * h44 + h52 * h54; 02169 m[A34] = h43 * h44 + h53 * h54; 02170 m[A44] = h44 * h44 + h54 * h54; 02171 m[A05] = h50 * h55; 02172 m[A15] = h51 * h55; 02173 m[A25] = h52 * h55; 02174 m[A35] = h53 * h55; 02175 m[A45] = h54 * h55; 02176 m[A55] = h55 * h55; 02177 02178 ifail = 0; 02179 return; 02180 }
void G4ErrorSymMatrix::invertHaywood4 | ( | G4int & | ifail | ) |
Definition at line 2260 of file G4ErrorSymMatrix.cc.
02261 { 02262 invert4(ifail); // For the 4x4 case, the method we use for invert is already 02263 // the Haywood method. 02264 }
void G4ErrorSymMatrix::invertHaywood5 | ( | G4int & | ifail | ) |
Definition at line 1358 of file G4ErrorSymMatrix.cc.
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, and A44.
01359 { 01360 ifail = 0; 01361 01362 // Find all NECESSARY 2x2 dets: (25 of them) 01363 01364 G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30]; 01365 G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30]; 01366 G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30]; 01367 G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31]; 01368 G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31]; 01369 G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32]; 01370 G4double Det2_24_01 = m[A20]*m[A41] - m[A21]*m[A40]; 01371 G4double Det2_24_02 = m[A20]*m[A42] - m[A22]*m[A40]; 01372 G4double Det2_24_03 = m[A20]*m[A43] - m[A23]*m[A40]; 01373 G4double Det2_24_04 = m[A20]*m[A44] - m[A24]*m[A40]; 01374 G4double Det2_24_12 = m[A21]*m[A42] - m[A22]*m[A41]; 01375 G4double Det2_24_13 = m[A21]*m[A43] - m[A23]*m[A41]; 01376 G4double Det2_24_14 = m[A21]*m[A44] - m[A24]*m[A41]; 01377 G4double Det2_24_23 = m[A22]*m[A43] - m[A23]*m[A42]; 01378 G4double Det2_24_24 = m[A22]*m[A44] - m[A24]*m[A42]; 01379 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40]; 01380 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40]; 01381 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40]; 01382 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40]; 01383 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41]; 01384 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41]; 01385 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41]; 01386 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42]; 01387 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42]; 01388 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43]; 01389 01390 // Find all NECESSARY 3x3 dets: (30 of them) 01391 01392 G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02 01393 + m[A12]*Det2_23_01; 01394 G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03 01395 + m[A13]*Det2_23_01; 01396 G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03 01397 + m[A13]*Det2_23_02; 01398 G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13 01399 + m[A13]*Det2_23_12; 01400 G4double Det3_124_012 = m[A10]*Det2_24_12 - m[A11]*Det2_24_02 01401 + m[A12]*Det2_24_01; 01402 G4double Det3_124_013 = m[A10]*Det2_24_13 - m[A11]*Det2_24_03 01403 + m[A13]*Det2_24_01; 01404 G4double Det3_124_014 = m[A10]*Det2_24_14 - m[A11]*Det2_24_04 01405 + m[A14]*Det2_24_01; 01406 G4double Det3_124_023 = m[A10]*Det2_24_23 - m[A12]*Det2_24_03 01407 + m[A13]*Det2_24_02; 01408 G4double Det3_124_024 = m[A10]*Det2_24_24 - m[A12]*Det2_24_04 01409 + m[A14]*Det2_24_02; 01410 G4double Det3_124_123 = m[A11]*Det2_24_23 - m[A12]*Det2_24_13 01411 + m[A13]*Det2_24_12; 01412 G4double Det3_124_124 = m[A11]*Det2_24_24 - m[A12]*Det2_24_14 01413 + m[A14]*Det2_24_12; 01414 G4double Det3_134_012 = m[A10]*Det2_34_12 - m[A11]*Det2_34_02 01415 + m[A12]*Det2_34_01; 01416 G4double Det3_134_013 = m[A10]*Det2_34_13 - m[A11]*Det2_34_03 01417 + m[A13]*Det2_34_01; 01418 G4double Det3_134_014 = m[A10]*Det2_34_14 - m[A11]*Det2_34_04 01419 + m[A14]*Det2_34_01; 01420 G4double Det3_134_023 = m[A10]*Det2_34_23 - m[A12]*Det2_34_03 01421 + m[A13]*Det2_34_02; 01422 G4double Det3_134_024 = m[A10]*Det2_34_24 - m[A12]*Det2_34_04 01423 + m[A14]*Det2_34_02; 01424 G4double Det3_134_034 = m[A10]*Det2_34_34 - m[A13]*Det2_34_04 01425 + m[A14]*Det2_34_03; 01426 G4double Det3_134_123 = m[A11]*Det2_34_23 - m[A12]*Det2_34_13 01427 + m[A13]*Det2_34_12; 01428 G4double Det3_134_124 = m[A11]*Det2_34_24 - m[A12]*Det2_34_14 01429 + m[A14]*Det2_34_12; 01430 G4double Det3_134_134 = m[A11]*Det2_34_34 - m[A13]*Det2_34_14 01431 + m[A14]*Det2_34_13; 01432 G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02 01433 + m[A22]*Det2_34_01; 01434 G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03 01435 + m[A23]*Det2_34_01; 01436 G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04 01437 + m[A24]*Det2_34_01; 01438 G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03 01439 + m[A23]*Det2_34_02; 01440 G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04 01441 + m[A24]*Det2_34_02; 01442 G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04 01443 + m[A24]*Det2_34_03; 01444 G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13 01445 + m[A23]*Det2_34_12; 01446 G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14 01447 + m[A24]*Det2_34_12; 01448 G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14 01449 + m[A24]*Det2_34_13; 01450 G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24 01451 + m[A24]*Det2_34_23; 01452 01453 // Find all NECESSARY 4x4 dets: (15 of them) 01454 01455 G4double Det4_0123_0123 = m[A00]*Det3_123_123 - m[A01]*Det3_123_023 01456 + m[A02]*Det3_123_013 - m[A03]*Det3_123_012; 01457 G4double Det4_0124_0123 = m[A00]*Det3_124_123 - m[A01]*Det3_124_023 01458 + m[A02]*Det3_124_013 - m[A03]*Det3_124_012; 01459 G4double Det4_0124_0124 = m[A00]*Det3_124_124 - m[A01]*Det3_124_024 01460 + m[A02]*Det3_124_014 - m[A04]*Det3_124_012; 01461 G4double Det4_0134_0123 = m[A00]*Det3_134_123 - m[A01]*Det3_134_023 01462 + m[A02]*Det3_134_013 - m[A03]*Det3_134_012; 01463 G4double Det4_0134_0124 = m[A00]*Det3_134_124 - m[A01]*Det3_134_024 01464 + m[A02]*Det3_134_014 - m[A04]*Det3_134_012; 01465 G4double Det4_0134_0134 = m[A00]*Det3_134_134 - m[A01]*Det3_134_034 01466 + m[A03]*Det3_134_014 - m[A04]*Det3_134_013; 01467 G4double Det4_0234_0123 = m[A00]*Det3_234_123 - m[A01]*Det3_234_023 01468 + m[A02]*Det3_234_013 - m[A03]*Det3_234_012; 01469 G4double Det4_0234_0124 = m[A00]*Det3_234_124 - m[A01]*Det3_234_024 01470 + m[A02]*Det3_234_014 - m[A04]*Det3_234_012; 01471 G4double Det4_0234_0134 = m[A00]*Det3_234_134 - m[A01]*Det3_234_034 01472 + m[A03]*Det3_234_014 - m[A04]*Det3_234_013; 01473 G4double Det4_0234_0234 = m[A00]*Det3_234_234 - m[A02]*Det3_234_034 01474 + m[A03]*Det3_234_024 - m[A04]*Det3_234_023; 01475 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023 01476 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012; 01477 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024 01478 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012; 01479 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034 01480 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013; 01481 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034 01482 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023; 01483 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134 01484 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123; 01485 01486 // Find the 5x5 det: 01487 01488 G4double det = m[A00]*Det4_1234_1234 01489 - m[A01]*Det4_1234_0234 01490 + m[A02]*Det4_1234_0134 01491 - m[A03]*Det4_1234_0124 01492 + m[A04]*Det4_1234_0123; 01493 01494 if ( det == 0 ) 01495 { 01496 ifail = 1; 01497 return; 01498 } 01499 01500 G4double oneOverDet = 1.0/det; 01501 G4double mn1OverDet = - oneOverDet; 01502 01503 m[A00] = Det4_1234_1234 * oneOverDet; 01504 m[A01] = Det4_1234_0234 * mn1OverDet; 01505 m[A02] = Det4_1234_0134 * oneOverDet; 01506 m[A03] = Det4_1234_0124 * mn1OverDet; 01507 m[A04] = Det4_1234_0123 * oneOverDet; 01508 01509 m[A11] = Det4_0234_0234 * oneOverDet; 01510 m[A12] = Det4_0234_0134 * mn1OverDet; 01511 m[A13] = Det4_0234_0124 * oneOverDet; 01512 m[A14] = Det4_0234_0123 * mn1OverDet; 01513 01514 m[A22] = Det4_0134_0134 * oneOverDet; 01515 m[A23] = Det4_0134_0124 * mn1OverDet; 01516 m[A24] = Det4_0134_0123 * oneOverDet; 01517 01518 m[A33] = Det4_0124_0124 * oneOverDet; 01519 m[A34] = Det4_0124_0123 * mn1OverDet; 01520 01521 m[A44] = Det4_0123_0123 * oneOverDet; 01522 01523 return; 01524 }
void G4ErrorSymMatrix::invertHaywood6 | ( | G4int & | ifail | ) |
Definition at line 1526 of file G4ErrorSymMatrix.cc.
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, and A55.
01527 { 01528 ifail = 0; 01529 01530 // Find all NECESSARY 2x2 dets: (39 of them) 01531 01532 G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40]; 01533 G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40]; 01534 G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40]; 01535 G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40]; 01536 G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41]; 01537 G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41]; 01538 G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41]; 01539 G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42]; 01540 G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42]; 01541 G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43]; 01542 G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50]; 01543 G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50]; 01544 G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50]; 01545 G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50]; 01546 G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50]; 01547 G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51]; 01548 G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51]; 01549 G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51]; 01550 G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51]; 01551 G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52]; 01552 G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52]; 01553 G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52]; 01554 G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53]; 01555 G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53]; 01556 G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50]; 01557 G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50]; 01558 G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50]; 01559 G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50]; 01560 G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50]; 01561 G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51]; 01562 G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51]; 01563 G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51]; 01564 G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51]; 01565 G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52]; 01566 G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52]; 01567 G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52]; 01568 G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53]; 01569 G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53]; 01570 G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54]; 01571 01572 // Find all NECESSARY 3x3 dets: (65 of them) 01573 01574 G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02 01575 + m[A22]*Det2_34_01; 01576 G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03 01577 + m[A23]*Det2_34_01; 01578 G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04 01579 + m[A24]*Det2_34_01; 01580 G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03 01581 + m[A23]*Det2_34_02; 01582 G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04 01583 + m[A24]*Det2_34_02; 01584 G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04 01585 + m[A24]*Det2_34_03; 01586 G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13 01587 + m[A23]*Det2_34_12; 01588 G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14 01589 + m[A24]*Det2_34_12; 01590 G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14 01591 + m[A24]*Det2_34_13; 01592 G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24 01593 + m[A24]*Det2_34_23; 01594 G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02 01595 + m[A22]*Det2_35_01; 01596 G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03 01597 + m[A23]*Det2_35_01; 01598 G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04 01599 + m[A24]*Det2_35_01; 01600 G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05 01601 + m[A25]*Det2_35_01; 01602 G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03 01603 + m[A23]*Det2_35_02; 01604 G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04 01605 + m[A24]*Det2_35_02; 01606 G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05 01607 + m[A25]*Det2_35_02; 01608 G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04 01609 + m[A24]*Det2_35_03; 01610 G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05 01611 + m[A25]*Det2_35_03; 01612 G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13 01613 + m[A23]*Det2_35_12; 01614 G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14 01615 + m[A24]*Det2_35_12; 01616 G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15 01617 + m[A25]*Det2_35_12; 01618 G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14 01619 + m[A24]*Det2_35_13; 01620 G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15 01621 + m[A25]*Det2_35_13; 01622 G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24 01623 + m[A24]*Det2_35_23; 01624 G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25 01625 + m[A25]*Det2_35_23; 01626 G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02 01627 + m[A22]*Det2_45_01; 01628 G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03 01629 + m[A23]*Det2_45_01; 01630 G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04 01631 + m[A24]*Det2_45_01; 01632 G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05 01633 + m[A25]*Det2_45_01; 01634 G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03 01635 + m[A23]*Det2_45_02; 01636 G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04 01637 + m[A24]*Det2_45_02; 01638 G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05 01639 + m[A25]*Det2_45_02; 01640 G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04 01641 + m[A24]*Det2_45_03; 01642 G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05 01643 + m[A25]*Det2_45_03; 01644 G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05 01645 + m[A25]*Det2_45_04; 01646 G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13 01647 + m[A23]*Det2_45_12; 01648 G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14 01649 + m[A24]*Det2_45_12; 01650 G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15 01651 + m[A25]*Det2_45_12; 01652 G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14 01653 + m[A24]*Det2_45_13; 01654 G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15 01655 + m[A25]*Det2_45_13; 01656 G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15 01657 + m[A25]*Det2_45_14; 01658 G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24 01659 + m[A24]*Det2_45_23; 01660 G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25 01661 + m[A25]*Det2_45_23; 01662 G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25 01663 + m[A25]*Det2_45_24; 01664 G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02 01665 + m[A32]*Det2_45_01; 01666 G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03 01667 + m[A33]*Det2_45_01; 01668 G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04 01669 + m[A34]*Det2_45_01; 01670 G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05 01671 + m[A35]*Det2_45_01; 01672 G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03 01673 + m[A33]*Det2_45_02; 01674 G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04 01675 + m[A34]*Det2_45_02; 01676 G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05 01677 + m[A35]*Det2_45_02; 01678 G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04 01679 + m[A34]*Det2_45_03; 01680 G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05 01681 + m[A35]*Det2_45_03; 01682 G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05 01683 + m[A35]*Det2_45_04; 01684 G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13 01685 + m[A33]*Det2_45_12; 01686 G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14 01687 + m[A34]*Det2_45_12; 01688 G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15 01689 + m[A35]*Det2_45_12; 01690 G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14 01691 + m[A34]*Det2_45_13; 01692 G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15 01693 + m[A35]*Det2_45_13; 01694 G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15 01695 + m[A35]*Det2_45_14; 01696 G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24 01697 + m[A34]*Det2_45_23; 01698 G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25 01699 + m[A35]*Det2_45_23; 01700 G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25 01701 + m[A35]*Det2_45_24; 01702 G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35 01703 + m[A35]*Det2_45_34; 01704 01705 // Find all NECESSARY 4x4 dets: (55 of them) 01706 01707 G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023 01708 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012; 01709 G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024 01710 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012; 01711 G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034 01712 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013; 01713 G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034 01714 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023; 01715 G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134 01716 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123; 01717 G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023 01718 + m[A12]*Det3_235_013 - m[A13]*Det3_235_012; 01719 G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024 01720 + m[A12]*Det3_235_014 - m[A14]*Det3_235_012; 01721 G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025 01722 + m[A12]*Det3_235_015 - m[A15]*Det3_235_012; 01723 G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034 01724 + m[A13]*Det3_235_014 - m[A14]*Det3_235_013; 01725 G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035 01726 + m[A13]*Det3_235_015 - m[A15]*Det3_235_013; 01727 G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034 01728 + m[A13]*Det3_235_024 - m[A14]*Det3_235_023; 01729 G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035 01730 + m[A13]*Det3_235_025 - m[A15]*Det3_235_023; 01731 G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134 01732 + m[A13]*Det3_235_124 - m[A14]*Det3_235_123; 01733 G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135 01734 + m[A13]*Det3_235_125 - m[A15]*Det3_235_123; 01735 G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023 01736 + m[A12]*Det3_245_013 - m[A13]*Det3_245_012; 01737 G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024 01738 + m[A12]*Det3_245_014 - m[A14]*Det3_245_012; 01739 G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025 01740 + m[A12]*Det3_245_015 - m[A15]*Det3_245_012; 01741 G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034 01742 + m[A13]*Det3_245_014 - m[A14]*Det3_245_013; 01743 G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035 01744 + m[A13]*Det3_245_015 - m[A15]*Det3_245_013; 01745 G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045 01746 + m[A14]*Det3_245_015 - m[A15]*Det3_245_014; 01747 G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034 01748 + m[A13]*Det3_245_024 - m[A14]*Det3_245_023; 01749 G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035 01750 + m[A13]*Det3_245_025 - m[A15]*Det3_245_023; 01751 G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045 01752 + m[A14]*Det3_245_025 - m[A15]*Det3_245_024; 01753 G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134 01754 + m[A13]*Det3_245_124 - m[A14]*Det3_245_123; 01755 G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135 01756 + m[A13]*Det3_245_125 - m[A15]*Det3_245_123; 01757 G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145 01758 + m[A14]*Det3_245_125 - m[A15]*Det3_245_124; 01759 G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023 01760 + m[A12]*Det3_345_013 - m[A13]*Det3_345_012; 01761 G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024 01762 + m[A12]*Det3_345_014 - m[A14]*Det3_345_012; 01763 G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025 01764 + m[A12]*Det3_345_015 - m[A15]*Det3_345_012; 01765 G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034 01766 + m[A13]*Det3_345_014 - m[A14]*Det3_345_013; 01767 G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035 01768 + m[A13]*Det3_345_015 - m[A15]*Det3_345_013; 01769 G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045 01770 + m[A14]*Det3_345_015 - m[A15]*Det3_345_014; 01771 G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034 01772 + m[A13]*Det3_345_024 - m[A14]*Det3_345_023; 01773 G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035 01774 + m[A13]*Det3_345_025 - m[A15]*Det3_345_023; 01775 G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045 01776 + m[A14]*Det3_345_025 - m[A15]*Det3_345_024; 01777 G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045 01778 + m[A14]*Det3_345_035 - m[A15]*Det3_345_034; 01779 G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134 01780 + m[A13]*Det3_345_124 - m[A14]*Det3_345_123; 01781 G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135 01782 + m[A13]*Det3_345_125 - m[A15]*Det3_345_123; 01783 G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145 01784 + m[A14]*Det3_345_125 - m[A15]*Det3_345_124; 01785 G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145 01786 + m[A14]*Det3_345_135 - m[A15]*Det3_345_134; 01787 G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023 01788 + m[A22]*Det3_345_013 - m[A23]*Det3_345_012; 01789 G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024 01790 + m[A22]*Det3_345_014 - m[A24]*Det3_345_012; 01791 G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025 01792 + m[A22]*Det3_345_015 - m[A25]*Det3_345_012; 01793 G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034 01794 + m[A23]*Det3_345_014 - m[A24]*Det3_345_013; 01795 G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035 01796 + m[A23]*Det3_345_015 - m[A25]*Det3_345_013; 01797 G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045 01798 + m[A24]*Det3_345_015 - m[A25]*Det3_345_014; 01799 G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034 01800 + m[A23]*Det3_345_024 - m[A24]*Det3_345_023; 01801 G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035 01802 + m[A23]*Det3_345_025 - m[A25]*Det3_345_023; 01803 G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045 01804 + m[A24]*Det3_345_025 - m[A25]*Det3_345_024; 01805 G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045 01806 + m[A24]*Det3_345_035 - m[A25]*Det3_345_034; 01807 G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134 01808 + m[A23]*Det3_345_124 - m[A24]*Det3_345_123; 01809 G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135 01810 + m[A23]*Det3_345_125 - m[A25]*Det3_345_123; 01811 G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145 01812 + m[A24]*Det3_345_125 - m[A25]*Det3_345_124; 01813 G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145 01814 + m[A24]*Det3_345_135 - m[A25]*Det3_345_134; 01815 G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245 01816 + m[A24]*Det3_345_235 - m[A25]*Det3_345_234; 01817 01818 // Find all NECESSARY 5x5 dets: (19 of them) 01819 01820 G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234 01821 + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123; 01822 G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234 01823 + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123; 01824 G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235 01825 + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123; 01826 G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234 01827 + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123; 01828 G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235 01829 + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123; 01830 G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245 01831 + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124; 01832 G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234 01833 + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123; 01834 G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235 01835 + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123; 01836 G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245 01837 + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124; 01838 G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345 01839 + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134; 01840 G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234 01841 + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123; 01842 G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235 01843 + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123; 01844 G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245 01845 + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124; 01846 G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345 01847 + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134; 01848 G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345 01849 + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234; 01850 G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234 01851 + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123; 01852 G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235 01853 + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123; 01854 G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245 01855 + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124; 01856 G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345 01857 + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134; 01858 G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345 01859 + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234; 01860 G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345 01861 + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234; 01862 01863 // Find the determinant 01864 01865 G4double det = m[A00]*Det5_12345_12345 01866 - m[A01]*Det5_12345_02345 01867 + m[A02]*Det5_12345_01345 01868 - m[A03]*Det5_12345_01245 01869 + m[A04]*Det5_12345_01235 01870 - m[A05]*Det5_12345_01234; 01871 01872 if ( det == 0 ) 01873 { 01874 ifail = 1; 01875 return; 01876 } 01877 01878 G4double oneOverDet = 1.0/det; 01879 G4double mn1OverDet = - oneOverDet; 01880 01881 m[A00] = Det5_12345_12345*oneOverDet; 01882 m[A01] = Det5_12345_02345*mn1OverDet; 01883 m[A02] = Det5_12345_01345*oneOverDet; 01884 m[A03] = Det5_12345_01245*mn1OverDet; 01885 m[A04] = Det5_12345_01235*oneOverDet; 01886 m[A05] = Det5_12345_01234*mn1OverDet; 01887 01888 m[A11] = Det5_02345_02345*oneOverDet; 01889 m[A12] = Det5_02345_01345*mn1OverDet; 01890 m[A13] = Det5_02345_01245*oneOverDet; 01891 m[A14] = Det5_02345_01235*mn1OverDet; 01892 m[A15] = Det5_02345_01234*oneOverDet; 01893 01894 m[A22] = Det5_01345_01345*oneOverDet; 01895 m[A23] = Det5_01345_01245*mn1OverDet; 01896 m[A24] = Det5_01345_01235*oneOverDet; 01897 m[A25] = Det5_01345_01234*mn1OverDet; 01898 01899 m[A33] = Det5_01245_01245*oneOverDet; 01900 m[A34] = Det5_01245_01235*mn1OverDet; 01901 m[A35] = Det5_01245_01234*oneOverDet; 01902 01903 m[A44] = Det5_01235_01235*oneOverDet; 01904 m[A45] = Det5_01235_01234*mn1OverDet; 01905 01906 m[A55] = Det5_01234_01234*oneOverDet; 01907 01908 return; 01909 }
G4int G4ErrorSymMatrix::num_col | ( | ) | const [inline] |
Definition at line 43 of file G4ErrorSymMatrix.icc.
Referenced by operator *(), operator+(), operator+=(), G4ErrorMatrix::operator+=(), operator-(), operator-=(), G4ErrorMatrix::operator-=(), operator<<(), and similarity().
G4int G4ErrorSymMatrix::num_row | ( | ) | const [inline] |
Definition at line 38 of file G4ErrorSymMatrix.icc.
Referenced by apply(), dsum(), operator *(), operator+(), operator+=(), G4ErrorMatrix::operator+=(), operator-(), operator-=(), G4ErrorMatrix::operator-=(), operator<<(), similarity(), and sub().
G4int G4ErrorSymMatrix::num_size | ( | ) | const [inline, protected] |
G4ErrorSymMatrix & G4ErrorSymMatrix::operator *= | ( | G4double | t | ) |
Definition at line 472 of file G4ErrorSymMatrix.cc.
References SIMPLE_UOP.
00473 { 00474 SIMPLE_UOP(*=) 00475 return (*this); 00476 }
G4ErrorSymMatrix & G4ErrorSymMatrix::operator+= | ( | const G4ErrorSymMatrix & | m2 | ) |
Definition at line 427 of file G4ErrorSymMatrix.cc.
References CHK_DIM_2, num_col(), num_row(), and SIMPLE_BOP.
00428 { 00429 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=); 00430 SIMPLE_BOP(+=) 00431 return (*this); 00432 }
G4ErrorSymMatrix G4ErrorSymMatrix::operator- | ( | ) | const |
Definition at line 197 of file G4ErrorSymMatrix.cc.
References m, and num_size().
00198 { 00199 G4ErrorSymMatrix mat2(nrow); 00200 G4ErrorMatrixConstIter a=m.begin(); 00201 G4ErrorMatrixIter b=mat2.m.begin(); 00202 G4ErrorMatrixConstIter e=m.begin()+num_size(); 00203 for(;a<e; a++, b++) { (*b) = -(*a); } 00204 return mat2; 00205 }
G4ErrorSymMatrix & G4ErrorSymMatrix::operator-= | ( | const G4ErrorSymMatrix & | m2 | ) |
Definition at line 459 of file G4ErrorSymMatrix.cc.
References CHK_DIM_2, num_col(), num_row(), and SIMPLE_BOP.
00460 { 00461 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=); 00462 SIMPLE_BOP(-=) 00463 return (*this); 00464 }
G4ErrorSymMatrix & G4ErrorSymMatrix::operator/= | ( | G4double | t | ) |
Definition at line 466 of file G4ErrorSymMatrix.cc.
References SIMPLE_UOP.
00467 { 00468 SIMPLE_UOP(/=) 00469 return (*this); 00470 }
G4ErrorSymMatrix & G4ErrorSymMatrix::operator= | ( | const G4ErrorSymMatrix & | m2 | ) |
Definition at line 509 of file G4ErrorSymMatrix.cc.
00510 { 00511 if (&mat1 == this) { return *this; } 00512 if(mat1.nrow != nrow) 00513 { 00514 nrow = mat1.nrow; 00515 size = mat1.size; 00516 m.resize(size); 00517 } 00518 m = mat1.m; 00519 return (*this); 00520 }
G4ErrorSymMatrix::G4ErrorSymMatrix_row_const G4ErrorSymMatrix::operator[] | ( | G4int | ) | const [inline] |
Definition at line 93 of file G4ErrorSymMatrix.icc.
00094 { 00095 const G4ErrorSymMatrix_row_const b(*this,r); 00096 return b; 00097 }
G4ErrorSymMatrix::G4ErrorSymMatrix_row G4ErrorSymMatrix::operator[] | ( | G4int | ) | [inline] |
Definition at line 86 of file G4ErrorSymMatrix.icc.
00087 { 00088 G4ErrorSymMatrix_row b(*this,r); 00089 return b; 00090 }
G4ErrorSymMatrix G4ErrorSymMatrix::similarity | ( | const G4ErrorSymMatrix & | m1 | ) | const |
Definition at line 620 of file G4ErrorSymMatrix.cc.
References m, CLHEP::detail::n, num_col(), and num_row().
00621 { 00622 G4ErrorSymMatrix mret(mat1.num_row()); 00623 G4ErrorMatrix temp = mat1*(*this); 00624 G4int n = mat1.num_col(); 00625 G4ErrorMatrixIter mr = mret.m.begin(); 00626 G4ErrorMatrixIter tempr1 = temp.m.begin(); 00627 for(G4int r=1;r<=mret.num_row();r++) 00628 { 00629 G4ErrorMatrixConstIter m1c1 = mat1.m.begin(); 00630 G4int c; 00631 for(c=1;c<=r;c++) 00632 { 00633 G4double tmp = 0.0; 00634 G4ErrorMatrixIter tempri = tempr1; 00635 G4ErrorMatrixConstIter m1ci = m1c1; 00636 G4int i; 00637 for(i=1;i<c;i++) 00638 { 00639 tmp+=(*(tempri++))*(*(m1ci++)); 00640 } 00641 for(i=c;i<=mat1.num_col();i++) 00642 { 00643 tmp+=(*(tempri++))*(*(m1ci)); 00644 m1ci += i; 00645 } 00646 *(mr++) = tmp; 00647 m1c1 += c; 00648 } 00649 tempr1 += n; 00650 } 00651 return mret; 00652 }
G4ErrorSymMatrix G4ErrorSymMatrix::similarity | ( | const G4ErrorMatrix & | m1 | ) | const |
Definition at line 589 of file G4ErrorSymMatrix.cc.
References G4ErrorMatrix::m, CLHEP::detail::n, G4ErrorMatrix::num_col(), and G4ErrorMatrix::num_row().
Referenced by G4ErrorSurfaceTrajState::BuildErrorMatrix(), G4ErrorFreeTrajState::G4ErrorFreeTrajState(), and G4ErrorFreeTrajState::PropagateError().
00590 { 00591 G4ErrorSymMatrix mret(mat1.num_row()); 00592 G4ErrorMatrix temp = mat1*(*this); 00593 00594 // If mat1*(*this) has correct dimensions, then so will the mat1.T multiplication. 00595 // So there is no need to check dimensions again. 00596 00597 G4int n = mat1.num_col(); 00598 G4ErrorMatrixIter mr = mret.m.begin(); 00599 G4ErrorMatrixIter tempr1 = temp.m.begin(); 00600 for(G4int r=1;r<=mret.num_row();r++) 00601 { 00602 G4ErrorMatrixConstIter m1c1 = mat1.m.begin(); 00603 for(G4int c=1;c<=r;c++) 00604 { 00605 G4double tmp = 0.0; 00606 G4ErrorMatrixIter tempri = tempr1; 00607 G4ErrorMatrixConstIter m1ci = m1c1; 00608 for(G4int i=1;i<=mat1.num_col();i++) 00609 { 00610 tmp+=(*(tempri++))*(*(m1ci++)); 00611 } 00612 *(mr++) = tmp; 00613 m1c1 += n; 00614 } 00615 tempr1 += n; 00616 } 00617 return mret; 00618 }
G4ErrorSymMatrix G4ErrorSymMatrix::similarityT | ( | const G4ErrorMatrix & | m1 | ) | const |
Definition at line 654 of file G4ErrorSymMatrix.cc.
References CLHEP::detail::n, and G4ErrorMatrix::num_col().
00655 { 00656 G4ErrorSymMatrix mret(mat1.num_col()); 00657 G4ErrorMatrix temp = (*this)*mat1; 00658 G4int n = mat1.num_col(); 00659 G4ErrorMatrixIter mrc = mret.m.begin(); 00660 G4ErrorMatrixIter temp1r = temp.m.begin(); 00661 for(G4int r=1;r<=mret.num_row();r++) 00662 { 00663 G4ErrorMatrixConstIter m11c = mat1.m.begin(); 00664 for(G4int c=1;c<=r;c++) 00665 { 00666 G4double tmp = 0.0; 00667 G4ErrorMatrixIter tempir = temp1r; 00668 G4ErrorMatrixConstIter m1ic = m11c; 00669 for(G4int i=1;i<=mat1.num_row();i++) 00670 { 00671 tmp+=(*(tempir))*(*(m1ic)); 00672 tempir += n; 00673 m1ic += n; 00674 } 00675 *(mrc++) = tmp; 00676 m11c++; 00677 } 00678 temp1r++; 00679 } 00680 return mret; 00681 }
G4ErrorSymMatrix G4ErrorSymMatrix::sub | ( | G4int | min_row, | |
G4int | max_row | |||
) |
Definition at line 143 of file G4ErrorSymMatrix.cc.
References G4ErrorMatrix::error(), m, and num_row().
00144 { 00145 G4ErrorSymMatrix mret(max_row-min_row+1); 00146 if(max_row > num_row()) 00147 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); } 00148 G4ErrorMatrixIter a = mret.m.begin(); 00149 G4ErrorMatrixIter b1 = m.begin() + (min_row+2)*(min_row-1)/2; 00150 for(G4int irow=1; irow<=mret.num_row(); irow++) 00151 { 00152 G4ErrorMatrixIter b = b1; 00153 for(G4int icol=1; icol<=irow; icol++) 00154 { 00155 *(a++) = *(b++); 00156 } 00157 b1 += irow+min_row-1; 00158 } 00159 return mret; 00160 }
void G4ErrorSymMatrix::sub | ( | G4int | row, | |
const G4ErrorSymMatrix & | m1 | |||
) |
Definition at line 162 of file G4ErrorSymMatrix.cc.
References G4ErrorMatrix::error(), m, and num_row().
00163 { 00164 if(row <1 || row+mat1.num_row()-1 > num_row() ) 00165 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); } 00166 G4ErrorMatrixConstIter a = mat1.m.begin(); 00167 G4ErrorMatrixIter b1 = m.begin() + (row+2)*(row-1)/2; 00168 for(G4int irow=1; irow<=mat1.num_row(); irow++) 00169 { 00170 G4ErrorMatrixIter b = b1; 00171 for(G4int icol=1; icol<=irow; icol++) 00172 { 00173 *(b++) = *(a++); 00174 } 00175 b1 += irow+row-1; 00176 } 00177 }
G4ErrorSymMatrix G4ErrorSymMatrix::sub | ( | G4int | min_row, | |
G4int | max_row | |||
) | const |
Definition at line 124 of file G4ErrorSymMatrix.cc.
References G4ErrorMatrix::error(), m, and num_row().
00125 { 00126 G4ErrorSymMatrix mret(max_row-min_row+1); 00127 if(max_row > num_row()) 00128 { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); } 00129 G4ErrorMatrixIter a = mret.m.begin(); 00130 G4ErrorMatrixConstIter b1 = m.begin() + (min_row+2)*(min_row-1)/2; 00131 for(G4int irow=1; irow<=mret.num_row(); irow++) 00132 { 00133 G4ErrorMatrixConstIter b = b1; 00134 for(G4int icol=1; icol<=irow; icol++) 00135 { 00136 *(a++) = *(b++); 00137 } 00138 b1 += irow+min_row-1; 00139 } 00140 return mret; 00141 }
G4ErrorSymMatrix G4ErrorSymMatrix::T | ( | ) | const [inline] |
Definition at line 80 of file G4ErrorSymMatrix.icc.
References G4ErrorSymMatrix().
Referenced by G4ErrorFreeTrajState::PropagateError().
00081 { 00082 return G4ErrorSymMatrix(*this); 00083 }
G4double G4ErrorSymMatrix::trace | ( | ) | const |
Definition at line 819 of file G4ErrorSymMatrix.cc.
00820 { 00821 G4double t = 0.0; 00822 for (G4int i=0; i<nrow; i++) 00823 { t += *(m.begin() + (i+3)*i/2); } 00824 return t; 00825 }
G4double condition | ( | const G4ErrorSymMatrix & | m | ) | [friend] |
void diag_step | ( | G4ErrorSymMatrix * | t, | |
G4ErrorMatrix * | u, | |||
G4int | begin, | |||
G4int | end | |||
) | [friend] |
void diag_step | ( | G4ErrorSymMatrix * | t, | |
G4int | begin, | |||
G4int | end | |||
) | [friend] |
G4ErrorMatrix diagonalize | ( | G4ErrorSymMatrix * | s | ) | [friend] |
friend class G4ErrorMatrix [friend] |
Definition at line 192 of file G4ErrorSymMatrix.hh.
friend class G4ErrorSymMatrix_row [friend] |
Definition at line 190 of file G4ErrorSymMatrix.hh.
friend class G4ErrorSymMatrix_row_const [friend] |
Definition at line 191 of file G4ErrorSymMatrix.hh.
void house_with_update2 | ( | G4ErrorSymMatrix * | a, | |
G4ErrorMatrix * | v, | |||
G4int | row, | |||
G4int | col | |||
) | [friend] |
G4ErrorMatrix operator * | ( | const G4ErrorMatrix & | m1, | |
const G4ErrorSymMatrix & | m2 | |||
) | [friend] |
Definition at line 287 of file G4ErrorSymMatrix.cc.
00288 { 00289 G4ErrorMatrix mret(mat1.num_row(),mat2.num_col()); 00290 CHK_DIM_1(mat1.num_col(),mat2.num_row(),*); 00291 G4ErrorMatrixConstIter mit1, mit2, sp,snp; //mit2=0 00292 G4double temp; 00293 G4ErrorMatrixIter mir=mret.m.begin(); 00294 for(mit1=mat1.m.begin(); 00295 mit1<mat1.m.begin()+mat1.num_row()*mat1.num_col(); 00296 mit1 = mit2) 00297 { 00298 snp=mat2.m.begin(); 00299 for(int step=1;step<=mat2.num_row();++step) 00300 { 00301 mit2=mit1; 00302 sp=snp; 00303 snp+=step; 00304 temp=0; 00305 while(sp<snp) 00306 temp+=*(sp++)*(*(mit2++)); 00307 if( step<mat2.num_row() ) { // only if we aren't on the last row 00308 sp+=step-1; 00309 for(int stept=step+1;stept<=mat2.num_row();stept++) 00310 { 00311 temp+=*sp*(*(mit2++)); 00312 if(stept<mat2.num_row()) sp+=stept; 00313 } 00314 } // if(step 00315 *(mir++)=temp; 00316 } // for(step 00317 } // for(mit1 00318 return mret; 00319 }
G4ErrorMatrix operator * | ( | const G4ErrorSymMatrix & | m1, | |
const G4ErrorMatrix & | m2 | |||
) | [friend] |
Definition at line 321 of file G4ErrorSymMatrix.cc.
00322 { 00323 G4ErrorMatrix mret(mat1.num_row(),mat2.num_col()); 00324 CHK_DIM_1(mat1.num_col(),mat2.num_row(),*); 00325 G4int step,stept; 00326 G4ErrorMatrixConstIter mit1,mit2,sp,snp; 00327 G4double temp; 00328 G4ErrorMatrixIter mir=mret.m.begin(); 00329 for(step=1,snp=mat1.m.begin();step<=mat1.num_row();snp+=step++) 00330 { 00331 for(mit1=mat2.m.begin();mit1<mat2.m.begin()+mat2.num_col();mit1++) 00332 { 00333 mit2=mit1; 00334 sp=snp; 00335 temp=0; 00336 while(sp<snp+step) 00337 { 00338 temp+=*mit2*(*(sp++)); 00339 mit2+=mat2.num_col(); 00340 } 00341 sp+=step-1; 00342 for(stept=step+1;stept<=mat1.num_row();stept++) 00343 { 00344 temp+=*mit2*(*sp); 00345 mit2+=mat2.num_col(); 00346 sp+=stept; 00347 } 00348 *(mir++)=temp; 00349 } 00350 } 00351 return mret; 00352 }
G4ErrorMatrix operator * | ( | const G4ErrorSymMatrix & | m1, | |
const G4ErrorSymMatrix & | m2 | |||
) | [friend] |
Definition at line 354 of file G4ErrorSymMatrix.cc.
00355 { 00356 G4ErrorMatrix mret(mat1.num_row(),mat1.num_row()); 00357 CHK_DIM_1(mat1.num_col(),mat2.num_row(),*); 00358 G4int step1,stept1,step2,stept2; 00359 G4ErrorMatrixConstIter snp1,sp1,snp2,sp2; 00360 G4double temp; 00361 G4ErrorMatrixIter mr = mret.m.begin(); 00362 for(step1=1,snp1=mat1.m.begin();step1<=mat1.num_row();snp1+=step1++) 00363 { 00364 for(step2=1,snp2=mat2.m.begin();step2<=mat2.num_row();) 00365 { 00366 sp1=snp1; 00367 sp2=snp2; 00368 snp2+=step2; 00369 temp=0; 00370 if(step1<step2) 00371 { 00372 while(sp1<snp1+step1) 00373 { temp+=(*(sp1++))*(*(sp2++)); } 00374 sp1+=step1-1; 00375 for(stept1=step1+1;stept1!=step2+1;sp1+=stept1++) 00376 { temp+=(*sp1)*(*(sp2++)); } 00377 sp2+=step2-1; 00378 for(stept2=++step2;stept2<=mat2.num_row();sp1+=stept1++,sp2+=stept2++) 00379 { temp+=(*sp1)*(*sp2); } 00380 } 00381 else 00382 { 00383 while(sp2<snp2) 00384 { temp+=(*(sp1++))*(*(sp2++)); } 00385 sp2+=step2-1; 00386 for(stept2=++step2;stept2!=step1+1;sp2+=stept2++) 00387 { temp+=(*(sp1++))*(*sp2); } 00388 sp1+=step1-1; 00389 for(stept1=step1+1;stept1<=mat1.num_row();sp1+=stept1++,sp2+=stept2++) 00390 { temp+=(*sp1)*(*sp2); } 00391 } 00392 *(mr++)=temp; 00393 } 00394 } 00395 return mret; 00396 }
G4ErrorSymMatrix operator+ | ( | const G4ErrorSymMatrix & | m1, | |
const G4ErrorSymMatrix & | m2 | |||
) | [friend] |
Definition at line 223 of file G4ErrorSymMatrix.cc.
00225 { 00226 G4ErrorSymMatrix mret(mat1.nrow); 00227 CHK_DIM_1(mat1.nrow, mat2.nrow,+); 00228 SIMPLE_TOP(+) 00229 return mret; 00230 }
G4ErrorSymMatrix operator- | ( | const G4ErrorSymMatrix & | m1, | |
const G4ErrorSymMatrix & | m2 | |||
) | [friend] |
Definition at line 252 of file G4ErrorSymMatrix.cc.
00254 { 00255 G4ErrorSymMatrix mret(mat1.num_row()); 00256 CHK_DIM_1(mat1.num_row(),mat2.num_row(),-); 00257 SIMPLE_TOP(-) 00258 return mret; 00259 }
void tridiagonal | ( | G4ErrorSymMatrix * | a, | |
G4ErrorMatrix * | hsm | |||
) | [friend] |