Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Classes | Public Member Functions | Protected Member Functions | Friends | List of all members
G4ErrorSymMatrix Class Reference

#include <G4ErrorSymMatrix.hh>

Classes

class  G4ErrorSymMatrix_row
 
class  G4ErrorSymMatrix_row_const
 

Public Member Functions

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

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

Here is the call graph for this function:

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.

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

Here is the call graph for this function:

void G4ErrorSymMatrix::assign ( const G4ErrorMatrix &  m2)

Definition at line 568 of file G4ErrorSymMatrix.cc.

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.

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
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4ErrorSymMatrix::invertBunchKaufman ( G4int &  ifail)

Definition at line 827 of file G4ErrorSymMatrix.cc.

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  G4Exception("G4ErrorSymMatrix::bunch_invert()",
1055  "GEANT4e-Notification", JustWarning,
1056  "Error in pivot choice!");
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  {
1156  std::ostringstream message;
1157  message << "Error in pivot: " << piv[j-1];
1158  G4Exception("G4ErrorSymMatrix::invertBunchKaufman()",
1159  "GEANT4e-Notification", JustWarning, message);
1160  }
1161  ss=2;
1162  if (j < nrow)
1163  {
1164  ip = m.begin() + (j+1)*j/2 + j-1;
1165  for (i=0; i < nrow-j; ip += 1+j+i++)
1166  {
1167  x[i] = *ip;
1168  }
1169  for (i=j+1; i<=nrow ; i++)
1170  {
1171  temp2 = 0;
1172  ip = m.begin() + i*(i-1)/2 + j;
1173  for (k=0; k <= i-j-1; k++)
1174  { temp2 += *ip++ * x[k]; }
1175  for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1176  { temp2 += *ip * x[k]; }
1177  *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
1178  }
1179  temp2 = 0;
1180  ip = m.begin() + (j+1)*j/2 + j-1;
1181  for (k=0; k < nrow-j; ip += 1+j+k++)
1182  { temp2 += x[k] * *ip; }
1183  *mjj -= temp2;
1184  temp2 = 0;
1185  ip = m.begin() + (j+1)*j/2 + j-2;
1186  for (i=j+1; i <= nrow; ip += i++)
1187  { temp2 += *ip * *(ip+1); }
1188  *(mjj-1) -= temp2;
1189  ip = m.begin() + (j+1)*j/2 + j-2;
1190  for (i=0; i < nrow-j; ip += 1+j+i++)
1191  {
1192  x[i] = *ip;
1193  }
1194  for (i=j+1; i <= nrow ; i++)
1195  {
1196  temp2 = 0;
1197  ip = m.begin() + i*(i-1)/2 + j;
1198  for (k=0; k <= i-j-1; k++)
1199  { temp2 += *ip++ * x[k]; }
1200  for (ip += i-1; k < nrow-j; ip += 1+j+k++)
1201  { temp2 += *ip * x[k]; }
1202  *(m.begin()+ i*(i-1)/2 + j-2)= -temp2;
1203  }
1204  temp2 = 0;
1205  ip = m.begin() + (j+1)*j/2 + j-2;
1206  for (k=0; k < nrow-j; ip += 1+j+k++)
1207  { temp2 += x[k] * *ip; }
1208  *(mjj-j) -= temp2;
1209  }
1210  }
1211 
1212  // interchange rows and columns j and piv[j-1]
1213  // or rows and columns j and -piv[j-2]
1214 
1215  pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
1216  ip = m.begin() + pivrow*(pivrow-1)/2 + j;
1217  for (i=j+1;i < pivrow; i++, ip++)
1218  {
1219  temp1 = *(m.begin() + i*(i-1)/2 + j-1);
1220  *(m.begin() + i*(i-1)/2 + j-1) = *ip;
1221  *ip = temp1;
1222  }
1223  temp1 = *mjj;
1224  *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
1225  *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
1226  if (ss==2)
1227  {
1228  temp1 = *(mjj-1);
1229  *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2);
1230  *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1;
1231  }
1232 
1233  ip = m.begin() + (pivrow+1)*pivrow/2 + j-1; // &A(i,j)
1234  iq = ip + pivrow-j;
1235  for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
1236  {
1237  temp1 = *iq;
1238  *iq = *ip;
1239  *ip = temp1;
1240  }
1241  } // end of loop over columns (in computing inverse from factorization)
1242 
1243  return; // inversion successful
1244 }
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:89
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
DBL_EPSILON
#define DBL_EPSILON
Definition: templates.hh:87
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4double
double G4double
Definition: G4Types.hh:76
alpha
static const G4double alpha
Definition: G4LivermoreBremsstrahlungModel.cc:79
epsilon
double epsilon(double density, double temperature)
Definition: G4DNAElectronHoleRecombination.cc:96

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ErrorSymMatrix::invertCholesky5 ( G4int &  ifail)

Definition at line 1918 of file G4ErrorSymMatrix.cc.

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

Definition at line 2040 of file G4ErrorSymMatrix.cc.

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

Definition at line 2267 of file G4ErrorSymMatrix.cc.

2268 {
2269  invert4(ifail); // For the 4x4 case, the method we use for invert is already
2270  // the Haywood method.
2271 }
void G4ErrorSymMatrix::invertHaywood5 ( G4int &  ifail)

Definition at line 1365 of file G4ErrorSymMatrix.cc.

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

Definition at line 1533 of file G4ErrorSymMatrix.cc.

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

Here is the caller graph for this function:

G4int G4ErrorSymMatrix::num_row ( ) const
inline

Here is the caller graph for this function:

G4int G4ErrorSymMatrix::num_size ( ) const
inlineprotected

Here is the caller graph for this function:

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.

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.

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

Here is the call graph for this function:

G4ErrorSymMatrix G4ErrorSymMatrix::operator- ( ) const

Definition at line 197 of file G4ErrorSymMatrix.cc.

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

Here is the call graph for this function:

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

Definition at line 459 of file G4ErrorSymMatrix.cc.

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

Here is the call graph for this function:

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

Definition at line 466 of file G4ErrorSymMatrix.cc.

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.

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
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 620 of file G4ErrorSymMatrix.cc.

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
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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

Definition at line 654 of file G4ErrorSymMatrix.cc.

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
G4ErrorMatrixIter
std::vector< G4double >::iterator G4ErrorMatrixIter
Definition: G4ErrorMatrix.hh:41
G4double
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

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

Definition at line 124 of file G4ErrorSymMatrix.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

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

Definition at line 162 of file G4ErrorSymMatrix.cc.

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

Here is the call graph for this function:

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

Definition at line 143 of file G4ErrorSymMatrix.cc.

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

Here is the call graph for this function:

G4ErrorSymMatrix G4ErrorSymMatrix::T ( ) const

Here is the caller graph for this function:

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) // Loop checking, 06.08.2015, G.Cosmo
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) // Loop checking, 06.08.2015, G.Cosmo
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) // Loop checking, 06.08.2015, G.Cosmo
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) // Loop checking, 06.08.2015, G.Cosmo
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: