Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Data Structures | Public Member Functions | Protected Member Functions | Friends
G4ErrorSymMatrix Class Reference

#include <G4ErrorSymMatrix.hh>

Data Structures

class  G4ErrorSymMatrix_row
 
class  G4ErrorSymMatrix_row_const
 

Public Member Functions

 G4ErrorSymMatrix ()
 
 G4ErrorSymMatrix (G4int p)
 
 G4ErrorSymMatrix (G4int p, G4int)
 
 G4ErrorSymMatrix (const G4ErrorSymMatrix &m1)
 
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- () 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::G4ErrorSymMatrix ( )
inline
G4ErrorSymMatrix::G4ErrorSymMatrix ( G4int  p)
explicit

Definition at line 71 of file G4ErrorSymMatrix.cc.

72  : m(p*(p+1)/2), nrow(p)
73 {
74  size = nrow * (nrow+1) / 2;
75  m.assign(size,0);
76 }
p
const char * p
Definition: xmltok.h:285
G4ErrorSymMatrix::G4ErrorSymMatrix ( G4int  p,
G4int  init 
)

Definition at line 78 of file G4ErrorSymMatrix.cc.

References test::a, and G4ErrorMatrix::error().

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

Definition at line 112 of file G4ErrorSymMatrix.cc.

113  : m(mat1.size), nrow(mat1.nrow), size(mat1.size)
114 {
115  m = mat1.m;
116 }
G4ErrorSymMatrix::~G4ErrorSymMatrix ( )
virtual

Definition at line 108 of file G4ErrorSymMatrix.cc.

109 {
110 }

Member Function Documentation

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

Definition at line 553 of file G4ErrorSymMatrix.cc.

References test::a, test::b, and num_row().

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

Definition at line 568 of file G4ErrorSymMatrix.cc.

References test::a, test::b, and test::c.

569 {
570  if(mat1.nrow != nrow)
571  {
572  nrow = mat1.nrow;
573  size = nrow * (nrow+1) / 2;
574  m.resize(size);
575  }
576  G4ErrorMatrixConstIter a = mat1.m.begin();
577  G4ErrorMatrixIter b = m.begin();
578  for(G4int r=1;r<=nrow;r++)
579  {
580  G4ErrorMatrixConstIter d = a;
581  for(G4int c=1;c<=r;c++)
582  {
583  *(b++) = *(d++);
584  }
585  a += nrow;
586  }
587 }
G4int
int G4int
Definition: G4Types.hh:78
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
void G4ErrorSymMatrix::assign ( const G4ErrorSymMatrix &  m2)
G4double G4ErrorSymMatrix::determinant ( ) const

Definition at line 799 of file G4ErrorSymMatrix.cc.

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

Definition at line 683 of file G4ErrorSymMatrix.cc.

References invertBunchKaufman(), and plottest35::t1.

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

Definition at line 827 of file G4ErrorSymMatrix.cc.

References DBL_EPSILON, G4cerr, G4endl, G4ThreadLocal, G4InuclParticleNames::lambda, and test::x.

Referenced by invert().

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

Definition at line 1913 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.

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

Definition at line 2035 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.

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

Definition at line 2262 of file G4ErrorSymMatrix.cc.

2263 {
2264  invert4(ifail); // For the 4x4 case, the method we use for invert is already
2265  // the Haywood method.
2266 }
void G4ErrorSymMatrix::invertHaywood5 ( G4int &  ifail)

Definition at line 1360 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.

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

Definition at line 1528 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.

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

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

G4int G4ErrorSymMatrix::num_row ( ) const
inline

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

G4int G4ErrorSymMatrix::num_size ( ) const
inlineprotected

Referenced by operator-().

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

Definition at line 472 of file G4ErrorSymMatrix.cc.

References SIMPLE_UOP.

473 {
474  SIMPLE_UOP(*=)
475  return (*this);
476 }
SIMPLE_UOP
#define SIMPLE_UOP(OPER)
Definition: G4ErrorSymMatrix.cc:41
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.

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

Definition at line 197 of file G4ErrorSymMatrix.cc.

References test::a, test::b, and num_size().

198 {
199  G4ErrorSymMatrix mat2(nrow);
200  G4ErrorMatrixConstIter a=m.begin();
201  G4ErrorMatrixIter b=mat2.m.begin();
202  G4ErrorMatrixConstIter e=m.begin()+num_size();
203  for(;a<e; a++, b++) { (*b) = -(*a); }
204  return mat2;
205 }
G4ErrorSymMatrix::num_size
G4int num_size() const
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
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.

460 {
461  CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
462  SIMPLE_BOP(-=)
463  return (*this);
464 }
SIMPLE_BOP
#define SIMPLE_BOP(OPER)
Definition: G4ErrorSymMatrix.cc:46
G4ErrorSymMatrix::num_col
G4int num_col() const
G4ErrorSymMatrix::num_row
G4int num_row() const
CHK_DIM_2
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition: G4ErrorSymMatrix.cc:59
G4ErrorSymMatrix & G4ErrorSymMatrix::operator/= ( G4double  t)

Definition at line 466 of file G4ErrorSymMatrix.cc.

References SIMPLE_UOP.

467 {
468  SIMPLE_UOP(/=)
469  return (*this);
470 }
SIMPLE_UOP
#define SIMPLE_UOP(OPER)
Definition: G4ErrorSymMatrix.cc:41
G4ErrorSymMatrix & G4ErrorSymMatrix::operator= ( const G4ErrorSymMatrix &  m2)

Definition at line 509 of file G4ErrorSymMatrix.cc.

510 {
511  if (&mat1 == this) { return *this; }
512  if(mat1.nrow != nrow)
513  {
514  nrow = mat1.nrow;
515  size = mat1.size;
516  m.resize(size);
517  }
518  m = mat1.m;
519  return (*this);
520 }
G4ErrorSymMatrix_row G4ErrorSymMatrix::operator[] ( G4int  )
inline
G4ErrorSymMatrix_row_const G4ErrorSymMatrix::operator[] ( G4int  ) const
inline
G4ErrorSymMatrix G4ErrorSymMatrix::similarity ( const G4ErrorMatrix &  m1) const

Definition at line 589 of file G4ErrorSymMatrix.cc.

References test::c, n, G4ErrorMatrix::num_col(), and G4ErrorMatrix::num_row().

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

590 {
591  G4ErrorSymMatrix mret(mat1.num_row());
592  G4ErrorMatrix temp = mat1*(*this);
593 
594 // If mat1*(*this) has correct dimensions, then so will the mat1.T multiplication.
595 // So there is no need to check dimensions again.
596 
597  G4int n = mat1.num_col();
598  G4ErrorMatrixIter mr = mret.m.begin();
599  G4ErrorMatrixIter tempr1 = temp.m.begin();
600  for(G4int r=1;r<=mret.num_row();r++)
601  {
602  G4ErrorMatrixConstIter m1c1 = mat1.m.begin();
603  for(G4int c=1;c<=r;c++)
604  {
605  G4double tmp = 0.0;
606  G4ErrorMatrixIter tempri = tempr1;
607  G4ErrorMatrixConstIter m1ci = m1c1;
608  for(G4int i=1;i<=mat1.num_col();i++)
609  {
610  tmp+=(*(tempri++))*(*(m1ci++));
611  }
612  *(mr++) = tmp;
613  m1c1 += n;
614  }
615  tempr1 += n;
616  }
617  return mret;
618 }
G4int
int G4int
Definition: G4Types.hh:78
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4double
double G4double
Definition: G4Types.hh:76
G4ErrorSymMatrix G4ErrorSymMatrix::similarity ( const G4ErrorSymMatrix &  m1) const

Definition at line 620 of file G4ErrorSymMatrix.cc.

References test::c, n, num_col(), and num_row().

621 {
622  G4ErrorSymMatrix mret(mat1.num_row());
623  G4ErrorMatrix temp = mat1*(*this);
624  G4int n = mat1.num_col();
625  G4ErrorMatrixIter mr = mret.m.begin();
626  G4ErrorMatrixIter tempr1 = temp.m.begin();
627  for(G4int r=1;r<=mret.num_row();r++)
628  {
629  G4ErrorMatrixConstIter m1c1 = mat1.m.begin();
630  G4int c;
631  for(c=1;c<=r;c++)
632  {
633  G4double tmp = 0.0;
634  G4ErrorMatrixIter tempri = tempr1;
635  G4ErrorMatrixConstIter m1ci = m1c1;
636  G4int i;
637  for(i=1;i<c;i++)
638  {
639  tmp+=(*(tempri++))*(*(m1ci++));
640  }
641  for(i=c;i<=mat1.num_col();i++)
642  {
643  tmp+=(*(tempri++))*(*(m1ci));
644  m1ci += i;
645  }
646  *(mr++) = tmp;
647  m1c1 += c;
648  }
649  tempr1 += n;
650  }
651  return mret;
652 }
G4int
int G4int
Definition: G4Types.hh:78
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4double
double G4double
Definition: G4Types.hh:76
G4ErrorSymMatrix G4ErrorSymMatrix::similarityT ( const G4ErrorMatrix &  m1) const

Definition at line 654 of file G4ErrorSymMatrix.cc.

References test::c, n, and G4ErrorMatrix::num_col().

655 {
656  G4ErrorSymMatrix mret(mat1.num_col());
657  G4ErrorMatrix temp = (*this)*mat1;
658  G4int n = mat1.num_col();
659  G4ErrorMatrixIter mrc = mret.m.begin();
660  G4ErrorMatrixIter temp1r = temp.m.begin();
661  for(G4int r=1;r<=mret.num_row();r++)
662  {
663  G4ErrorMatrixConstIter m11c = mat1.m.begin();
664  for(G4int c=1;c<=r;c++)
665  {
666  G4double tmp = 0.0;
667  G4ErrorMatrixIter tempir = temp1r;
668  G4ErrorMatrixConstIter m1ic = m11c;
669  for(G4int i=1;i<=mat1.num_row();i++)
670  {
671  tmp+=(*(tempir))*(*(m1ic));
672  tempir += n;
673  m1ic += n;
674  }
675  *(mrc++) = tmp;
676  m11c++;
677  }
678  temp1r++;
679  }
680  return mret;
681 }
G4int
int G4int
Definition: G4Types.hh:78
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4double
double G4double
Definition: G4Types.hh:76
G4ErrorSymMatrix G4ErrorSymMatrix::sub ( G4int  min_row,
G4int  max_row 
) const

Definition at line 124 of file G4ErrorSymMatrix.cc.

References test::a, test::b, G4ErrorMatrix::error(), and num_row().

Referenced by dsum().

125 {
126  G4ErrorSymMatrix mret(max_row-min_row+1);
127  if(max_row > num_row())
128  { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
129  G4ErrorMatrixIter a = mret.m.begin();
130  G4ErrorMatrixConstIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
131  for(G4int irow=1; irow<=mret.num_row(); irow++)
132  {
133  G4ErrorMatrixConstIter b = b1;
134  for(G4int icol=1; icol<=irow; icol++)
135  {
136  *(a++) = *(b++);
137  }
138  b1 += irow+min_row-1;
139  }
140  return mret;
141 }
G4int
int G4int
Definition: G4Types.hh:78
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
G4ErrorSymMatrix::num_row
G4int num_row() const
G4ErrorMatrix::error
static void error(const char *s)
Definition: G4ErrorMatrix.cc:805
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
void G4ErrorSymMatrix::sub ( G4int  row,
const G4ErrorSymMatrix &  m1 
)

Definition at line 162 of file G4ErrorSymMatrix.cc.

References test::a, test::b, G4ErrorMatrix::error(), and num_row().

163 {
164  if(row <1 || row+mat1.num_row()-1 > num_row() )
165  { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
166  G4ErrorMatrixConstIter a = mat1.m.begin();
167  G4ErrorMatrixIter b1 = m.begin() + (row+2)*(row-1)/2;
168  for(G4int irow=1; irow<=mat1.num_row(); irow++)
169  {
170  G4ErrorMatrixIter b = b1;
171  for(G4int icol=1; icol<=irow; icol++)
172  {
173  *(b++) = *(a++);
174  }
175  b1 += irow+row-1;
176  }
177 }
G4int
int G4int
Definition: G4Types.hh:78
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
G4ErrorSymMatrix::num_row
G4int num_row() const
G4ErrorMatrix::error
static void error(const char *s)
Definition: G4ErrorMatrix.cc:805
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4ErrorSymMatrix G4ErrorSymMatrix::sub ( G4int  min_row,
G4int  max_row 
)

Definition at line 143 of file G4ErrorSymMatrix.cc.

References test::a, test::b, G4ErrorMatrix::error(), and num_row().

144 {
145  G4ErrorSymMatrix mret(max_row-min_row+1);
146  if(max_row > num_row())
147  { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
148  G4ErrorMatrixIter a = mret.m.begin();
149  G4ErrorMatrixIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
150  for(G4int irow=1; irow<=mret.num_row(); irow++)
151  {
152  G4ErrorMatrixIter b = b1;
153  for(G4int icol=1; icol<=irow; icol++)
154  {
155  *(a++) = *(b++);
156  }
157  b1 += irow+min_row-1;
158  }
159  return mret;
160 }
G4int
int G4int
Definition: G4Types.hh:78
G4ErrorSymMatrix::num_row
G4int num_row() const
G4ErrorMatrix::error
static void error(const char *s)
Definition: G4ErrorMatrix.cc:805
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4ErrorSymMatrix G4ErrorSymMatrix::T ( ) const

Referenced by G4ErrorFreeTrajState::PropagateError().

G4double G4ErrorSymMatrix::trace ( ) const

Definition at line 819 of file G4ErrorSymMatrix.cc.

820 {
821  G4double t = 0.0;
822  for (G4int i=0; i<nrow; i++)
823  { t += *(m.begin() + (i+3)*i/2); }
824  return t;
825 }
G4int
int G4int
Definition: G4Types.hh:78
G4double
double G4double
Definition: G4Types.hh:76

Friends And Related Function Documentation

G4double condition ( const G4ErrorSymMatrix &  m)
friend
void diag_step ( G4ErrorSymMatrix *  t,
G4int  begin,
G4int  end 
)
friend
void diag_step ( G4ErrorSymMatrix *  t,
G4ErrorMatrix *  u,
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 G4ErrorSymMatrix &  m1,
const G4ErrorSymMatrix &  m2 
)
friend

Definition at line 354 of file G4ErrorSymMatrix.cc.

355 {
356  G4ErrorMatrix mret(mat1.num_row(),mat1.num_row());
357  CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
358  G4int step1,stept1,step2,stept2;
359  G4ErrorMatrixConstIter snp1,sp1,snp2,sp2;
360  G4double temp;
361  G4ErrorMatrixIter mr = mret.m.begin();
362  for(step1=1,snp1=mat1.m.begin();step1<=mat1.num_row();snp1+=step1++)
363  {
364  for(step2=1,snp2=mat2.m.begin();step2<=mat2.num_row();)
365  {
366  sp1=snp1;
367  sp2=snp2;
368  snp2+=step2;
369  temp=0;
370  if(step1<step2)
371  {
372  while(sp1<snp1+step1)
373  { temp+=(*(sp1++))*(*(sp2++)); }
374  sp1+=step1-1;
375  for(stept1=step1+1;stept1!=step2+1;sp1+=stept1++)
376  { temp+=(*sp1)*(*(sp2++)); }
377  sp2+=step2-1;
378  for(stept2=++step2;stept2<=mat2.num_row();sp1+=stept1++,sp2+=stept2++)
379  { temp+=(*sp1)*(*sp2); }
380  }
381  else
382  {
383  while(sp2<snp2)
384  { temp+=(*(sp1++))*(*(sp2++)); }
385  sp2+=step2-1;
386  for(stept2=++step2;stept2!=step1+1;sp2+=stept2++)
387  { temp+=(*(sp1++))*(*sp2); }
388  sp1+=step1-1;
389  for(stept1=step1+1;stept1<=mat1.num_row();sp1+=stept1++,sp2+=stept2++)
390  { temp+=(*sp1)*(*sp2); }
391  }
392  *(mr++)=temp;
393  }
394  }
395  return mret;
396 }
G4int
int G4int
Definition: G4Types.hh:78
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: G4ErrorSymMatrix.cc:64
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4double
double G4double
Definition: G4Types.hh:76
G4ErrorMatrix operator* ( const G4ErrorSymMatrix &  m1,
const G4ErrorMatrix &  m2 
)
friend

Definition at line 321 of file G4ErrorSymMatrix.cc.

322 {
323  G4ErrorMatrix mret(mat1.num_row(),mat2.num_col());
324  CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
325  G4int step,stept;
326  G4ErrorMatrixConstIter mit1,mit2,sp,snp;
327  G4double temp;
328  G4ErrorMatrixIter mir=mret.m.begin();
329  for(step=1,snp=mat1.m.begin();step<=mat1.num_row();snp+=step++)
330  {
331  for(mit1=mat2.m.begin();mit1<mat2.m.begin()+mat2.num_col();mit1++)
332  {
333  mit2=mit1;
334  sp=snp;
335  temp=0;
336  while(sp<snp+step)
337  {
338  temp+=*mit2*(*(sp++));
339  mit2+=mat2.num_col();
340  }
341  sp+=step-1;
342  for(stept=step+1;stept<=mat1.num_row();stept++)
343  {
344  temp+=*mit2*(*sp);
345  mit2+=mat2.num_col();
346  sp+=stept;
347  }
348  *(mir++)=temp;
349  }
350  }
351  return mret;
352 }
G4int
int G4int
Definition: G4Types.hh:78
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: G4ErrorSymMatrix.cc:64
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4double
double G4double
Definition: G4Types.hh:76
G4ErrorMatrix operator* ( const G4ErrorMatrix &  m1,
const G4ErrorSymMatrix &  m2 
)
friend

Definition at line 287 of file G4ErrorSymMatrix.cc.

288 {
289  G4ErrorMatrix mret(mat1.num_row(),mat2.num_col());
290  CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
291  G4ErrorMatrixConstIter mit1, mit2, sp,snp; //mit2=0
292  G4double temp;
293  G4ErrorMatrixIter mir=mret.m.begin();
294  for(mit1=mat1.m.begin();
295  mit1<mat1.m.begin()+mat1.num_row()*mat1.num_col();
296  mit1 = mit2)
297  {
298  snp=mat2.m.begin();
299  for(int step=1;step<=mat2.num_row();++step)
300  {
301  mit2=mit1;
302  sp=snp;
303  snp+=step;
304  temp=0;
305  while(sp<snp)
306  temp+=*(sp++)*(*(mit2++));
307  if( step<mat2.num_row() ) { // only if we aren't on the last row
308  sp+=step-1;
309  for(int stept=step+1;stept<=mat2.num_row();stept++)
310  {
311  temp+=*sp*(*(mit2++));
312  if(stept<mat2.num_row()) sp+=stept;
313  }
314  } // if(step
315  *(mir++)=temp;
316  } // for(step
317  } // for(mit1
318  return mret;
319 }
G4ErrorMatrixConstIter
std::vector< G4double >::const_iterator G4ErrorMatrixConstIter
Definition: G4ErrorMatrix.hh:44
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: G4ErrorSymMatrix.cc:64
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4double
double G4double
Definition: G4Types.hh:76
G4ErrorSymMatrix operator+ ( const G4ErrorSymMatrix &  m1,
const G4ErrorSymMatrix &  m2 
)
friend

Definition at line 223 of file G4ErrorSymMatrix.cc.

225 {
226  G4ErrorSymMatrix mret(mat1.nrow);
227  CHK_DIM_1(mat1.nrow, mat2.nrow,+);
228  SIMPLE_TOP(+)
229  return mret;
230 }
SIMPLE_TOP
#define SIMPLE_TOP(OPER)
Definition: G4ErrorSymMatrix.cc:52
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: G4ErrorSymMatrix.cc:64
G4ErrorSymMatrix operator- ( const G4ErrorSymMatrix &  m1,
const G4ErrorSymMatrix &  m2 
)
friend

Definition at line 252 of file G4ErrorSymMatrix.cc.

254 {
255  G4ErrorSymMatrix mret(mat1.num_row());
256  CHK_DIM_1(mat1.num_row(),mat2.num_row(),-);
257  SIMPLE_TOP(-)
258  return mret;
259 }
SIMPLE_TOP
#define SIMPLE_TOP(OPER)
Definition: G4ErrorSymMatrix.cc:52
CHK_DIM_1
#define CHK_DIM_1(c1, r2, fun)
Definition: G4ErrorSymMatrix.cc:64
void tridiagonal ( G4ErrorSymMatrix *  a,
G4ErrorMatrix *  hsm 
)
friend
The documentation for this class was generated from the following files:
Generated on Wed Apr 30 2014 15:56:24 for Geant4.10 by   doxygen 1.8.7