G4ErrorSymMatrix Class Reference

#include <G4ErrorSymMatrix.hh>


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 G4doubleoperator() (G4int row, G4int col) const
G4doubleoperator() (G4int row, G4int col)
const G4doublefast (G4int row, G4int col) const
G4doublefast (G4int row, G4int col)
void assign (const G4ErrorMatrix &m2)
void assign (const G4ErrorSymMatrix &m2)
G4ErrorSymMatrixoperator *= (G4double t)
G4ErrorSymMatrixoperator/= (G4double t)
G4ErrorSymMatrixoperator+= (const G4ErrorSymMatrix &m2)
G4ErrorSymMatrixoperator-= (const G4ErrorSymMatrix &m2)
G4ErrorSymMatrixoperator= (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)

Data Structures

class  G4ErrorSymMatrix_row
class  G4ErrorSymMatrix_row_const


Detailed Description

Definition at line 43 of file G4ErrorSymMatrix.hh.


Constructor & Destructor Documentation

G4ErrorSymMatrix::G4ErrorSymMatrix (  )  [inline]

Definition at line 33 of file G4ErrorSymMatrix.icc.

Referenced by T().

00034   : m(0), nrow(0), size(0)
00035 {
00036 }

G4ErrorSymMatrix::G4ErrorSymMatrix ( G4int  p  )  [explicit]

Definition at line 71 of file G4ErrorSymMatrix.cc.

00072    : m(p*(p+1)/2), nrow(p)
00073 {
00074    size = nrow * (nrow+1) / 2;
00075    m.assign(size,0);
00076 }

G4ErrorSymMatrix::G4ErrorSymMatrix ( G4int  p,
G4int   
)

Definition at line 78 of file G4ErrorSymMatrix.cc.

References G4ErrorMatrix::error().

00079    : m(p*(p+1)/2), nrow(p)
00080 {
00081    size = nrow * (nrow+1) / 2;
00082 
00083    m.assign(size,0);
00084    switch(init)
00085    {
00086    case 0:
00087       break;
00088       
00089    case 1:
00090       {
00091          G4ErrorMatrixIter a = m.begin();
00092          for(G4int i=1;i<=nrow;i++)
00093          {
00094             *a = 1.0;
00095             a += (i+1);
00096          }
00097          break;
00098       }
00099    default:
00100       G4ErrorMatrix::error("G4ErrorSymMatrix: initialization must be 0 or 1.");
00101    }
00102 }

G4ErrorSymMatrix::G4ErrorSymMatrix ( const G4ErrorSymMatrix m1  ) 

Definition at line 112 of file G4ErrorSymMatrix.cc.

References m.

00113    : m(mat1.size), nrow(mat1.nrow), size(mat1.size)
00114 {
00115    m = mat1.m;
00116 }

G4ErrorSymMatrix::~G4ErrorSymMatrix (  )  [virtual]

Definition at line 108 of file G4ErrorSymMatrix.cc.

00109 {
00110 }


Member Function Documentation

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

Definition at line 553 of file G4ErrorSymMatrix.cc.

References m, and num_row().

00554 {
00555   G4ErrorSymMatrix mret(num_row());
00556   G4ErrorMatrixConstIter a = m.begin();
00557   G4ErrorMatrixIter b = mret.m.begin();
00558   for(G4int ir=1;ir<=num_row();ir++)
00559   {
00560     for(G4int ic=1;ic<=ir;ic++)
00561     {
00562       *(b++) = (*f)(*(a++), ir, ic);
00563     }
00564   }
00565   return mret;
00566 }

void G4ErrorSymMatrix::assign ( const G4ErrorSymMatrix m2  )  [inline]

Definition at line 75 of file G4ErrorSymMatrix.icc.

00076 {
00077   (*this)=mat;
00078 }

void G4ErrorSymMatrix::assign ( const G4ErrorMatrix m2  ) 

Definition at line 568 of file G4ErrorSymMatrix.cc.

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

00569 {
00570    if(mat1.nrow != nrow)
00571    {
00572       nrow = mat1.nrow;
00573       size = nrow * (nrow+1) / 2;
00574       m.resize(size);
00575    }
00576    G4ErrorMatrixConstIter a = mat1.m.begin();
00577    G4ErrorMatrixIter b = m.begin();
00578    for(G4int r=1;r<=nrow;r++)
00579    {
00580       G4ErrorMatrixConstIter d = a;
00581       for(G4int c=1;c<=r;c++)
00582       {
00583          *(b++) = *(d++);
00584       }
00585       a += nrow;
00586    }
00587 }

G4double G4ErrorSymMatrix::determinant (  )  const

Definition at line 799 of file G4ErrorSymMatrix.cc.

References G4ErrorMatrix::dfact_matrix().

00800 {
00801   static const G4int max_array = 20;
00802 
00803   // ir must point to an array which is ***1 longer than*** nrow
00804 
00805   static std::vector<G4int> ir_vec (max_array+1); 
00806   if (ir_vec.size() <= static_cast<unsigned int>(nrow))
00807   {
00808     ir_vec.resize(nrow+1);
00809   }
00810   G4int * ir = &ir_vec[0];   
00811 
00812   G4double det;
00813   G4ErrorMatrix mt(*this);
00814   G4int i = mt.dfact_matrix(det, ir);
00815   if(i==0) { return det; }
00816   return 0.0;
00817 }

G4double & G4ErrorSymMatrix::fast ( G4int  row,
G4int  col 
) [inline]

Definition at line 53 of file G4ErrorSymMatrix.icc.

00054 {
00055   return *(m.begin()+(row*(row-1))/2+(col-1));
00056 }

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

Definition at line 58 of file G4ErrorSymMatrix.icc.

Referenced by operator()().

00059 {
00060   return *(m.begin()+(row*(row-1))/2+(col-1));
00061 }

G4ErrorSymMatrix G4ErrorSymMatrix::inverse ( G4int ifail  )  const [inline]

Definition at line 135 of file G4ErrorSymMatrix.icc.

References invert().

00136 {
00137   G4ErrorSymMatrix mTmp(*this);
00138   mTmp.invert(ifail);
00139   return mTmp;
00140 }

void G4ErrorSymMatrix::invert ( G4int ifail  ) 

Definition at line 683 of file G4ErrorSymMatrix.cc.

References invertBunchKaufman().

Referenced by inverse().

00684 {  
00685   ifail = 0;
00686 
00687   switch(nrow)
00688   {
00689   case 3:
00690     {
00691       G4double det, temp;
00692       G4double t1, t2, t3;
00693       G4double c11,c12,c13,c22,c23,c33;
00694       c11 = (*(m.begin()+2)) * (*(m.begin()+5))
00695           - (*(m.begin()+4)) * (*(m.begin()+4));
00696       c12 = (*(m.begin()+4)) * (*(m.begin()+3))
00697           - (*(m.begin()+1)) * (*(m.begin()+5));
00698       c13 = (*(m.begin()+1)) * (*(m.begin()+4))
00699           - (*(m.begin()+2)) * (*(m.begin()+3));
00700       c22 = (*(m.begin()+5)) * (*m.begin())
00701           - (*(m.begin()+3)) * (*(m.begin()+3));
00702       c23 = (*(m.begin()+3)) * (*(m.begin()+1))
00703           - (*(m.begin()+4)) * (*m.begin());
00704       c33 = (*m.begin()) * (*(m.begin()+2))
00705           - (*(m.begin()+1)) * (*(m.begin()+1));
00706       t1 = std::fabs(*m.begin());
00707       t2 = std::fabs(*(m.begin()+1));
00708       t3 = std::fabs(*(m.begin()+3));
00709       if (t1 >= t2)
00710       {
00711         if (t3 >= t1)
00712         {
00713           temp = *(m.begin()+3);
00714           det = c23*c12-c22*c13;
00715         }
00716         else
00717         {
00718           temp = *m.begin();
00719           det = c22*c33-c23*c23;
00720         }
00721       }
00722       else if (t3 >= t2)
00723       {
00724         temp = *(m.begin()+3);
00725         det = c23*c12-c22*c13;
00726       }
00727       else
00728       {
00729         temp = *(m.begin()+1);
00730         det = c13*c23-c12*c33;
00731       }
00732       if (det==0)
00733       {
00734         ifail = 1;
00735         return;
00736       }
00737       {
00738         G4double ss = temp/det;
00739         G4ErrorMatrixIter mq = m.begin();
00740         *(mq++) = ss*c11;
00741         *(mq++) = ss*c12;
00742         *(mq++) = ss*c22;
00743         *(mq++) = ss*c13;
00744         *(mq++) = ss*c23;
00745         *(mq) = ss*c33;
00746       }
00747     }
00748     break;
00749  case 2:
00750     {
00751       G4double det, temp, ss;
00752       det = (*m.begin())*(*(m.begin()+2)) - (*(m.begin()+1))*(*(m.begin()+1));
00753       if (det==0)
00754       {
00755         ifail = 1;
00756         return;
00757       }
00758       ss = 1.0/det;
00759       *(m.begin()+1) *= -ss;
00760       temp = ss*(*(m.begin()+2));
00761       *(m.begin()+2) = ss*(*m.begin());
00762       *m.begin() = temp;
00763       break;
00764     }
00765  case 1:
00766     {
00767       if ((*m.begin())==0)
00768       {
00769         ifail = 1;
00770         return;
00771       }
00772       *m.begin() = 1.0/(*m.begin());
00773       break;
00774     }
00775  case 5:
00776     {
00777       invert5(ifail);
00778       return;
00779     }
00780  case 6:
00781     {
00782       invert6(ifail);
00783       return;
00784     }
00785  case 4:
00786     {
00787       invert4(ifail);
00788       return;
00789     }
00790  default:
00791     {
00792       invertBunchKaufman(ifail);
00793       return;
00794     }
00795   }
00796   return; // inversion successful
00797 }

void G4ErrorSymMatrix::invertBunchKaufman ( G4int ifail  ) 

Definition at line 827 of file G4ErrorSymMatrix.cc.

References G4InuclParticleNames::alpha, DBL_EPSILON, G4cerr, G4endl, and G4InuclParticleNames::lambda.

Referenced by invert().

00828 {
00829   // Bunch-Kaufman diagonal pivoting method
00830   // It is decribed in J.R. Bunch, L. Kaufman (1977). 
00831   // "Some Stable Methods for Calculating Inertia and Solving Symmetric 
00832   // Linear Systems", Math. Comp. 31, p. 162-179. or in Gene H. Golub, 
00833   // Charles F. van Loan, "Matrix Computations" (the second edition 
00834   // has a bug.) and implemented in "lapack"
00835   // Mario Stanke, 09/97
00836 
00837   G4int i, j, k, ss;
00838   G4int pivrow;
00839 
00840   // Establish the two working-space arrays needed:  x and piv are
00841   // used as pointers to arrays of doubles and ints respectively, each
00842   // of length nrow.  We do not want to reallocate each time through
00843   // unless the size needs to grow.  We do not want to leak memory, even
00844   // by having a new without a delete that is only done once.
00845   
00846   static const G4int max_array = 25;
00847   static std::vector<G4double> xvec (max_array);
00848   static std::vector<G4int>    pivv (max_array);
00849   typedef std::vector<G4int>::iterator pivIter; 
00850   if (xvec.size() < static_cast<unsigned int>(nrow)) xvec.resize(nrow);
00851   if (pivv.size() < static_cast<unsigned int>(nrow)) pivv.resize(nrow);
00852     // Note - resize should do  nothing if the size is already larger than nrow,
00853     //        but on VC++ there are indications that it does so we check.
00854     // Note - the data elements in a vector are guaranteed to be contiguous,
00855     //        so x[i] and piv[i] are optimally fast.
00856   G4ErrorMatrixIter   x   = xvec.begin();
00857     // x[i] is used as helper storage, needs to have at least size nrow.
00858   pivIter piv = pivv.begin();
00859     // piv[i] is used to store details of exchanges
00860       
00861   G4double temp1, temp2;
00862   G4ErrorMatrixIter ip, mjj, iq;
00863   G4double lambda, sigma;
00864   const G4double alpha = .6404; // = (1+sqrt(17))/8
00865   const G4double epsilon = 32*DBL_EPSILON;
00866     // whenever a sum of two doubles is below or equal to epsilon
00867     // it is set to zero.
00868     // this constant could be set to zero but then the algorithm
00869     // doesn't neccessarily detect that a matrix is singular
00870   
00871   for (i = 0; i < nrow; i++)
00872   {
00873     piv[i] = i+1;
00874   }
00875       
00876   ifail = 0;
00877       
00878   // compute the factorization P*A*P^T = L * D * L^T 
00879   // L is unit lower triangular, D is direct sum of 1x1 and 2x2 matrices
00880   // L and D^-1 are stored in A = *this, P is stored in piv[]
00881         
00882   for (j=1; j < nrow; j+=ss)  // main loop over columns
00883   {
00884           mjj = m.begin() + j*(j-1)/2 + j-1;
00885           lambda = 0;           // compute lambda = max of A(j+1:n,j)
00886           pivrow = j+1;
00887           ip = m.begin() + (j+1)*j/2 + j-1;
00888           for (i=j+1; i <= nrow ; ip += i++)
00889           {
00890             if (std::fabs(*ip) > lambda)
00891               {
00892                 lambda = std::fabs(*ip);
00893                 pivrow = i;
00894               }
00895           }
00896           if (lambda == 0 )
00897             {
00898               if (*mjj == 0)
00899                 {
00900                   ifail = 1;
00901                   return;
00902                 }
00903               ss=1;
00904               *mjj = 1./ *mjj;
00905             }
00906           else
00907             {
00908               if (std::fabs(*mjj) >= lambda*alpha)
00909                 {
00910                   ss=1;
00911                   pivrow=j;
00912                 }
00913               else
00914                 {
00915                   sigma = 0;  // compute sigma = max A(pivrow, j:pivrow-1)
00916                   ip = m.begin() + pivrow*(pivrow-1)/2+j-1;
00917                   for (k=j; k < pivrow; k++)
00918                     {
00919                       if (std::fabs(*ip) > sigma)
00920                         sigma = std::fabs(*ip);
00921                       ip++;
00922                     }
00923                   if (sigma * std::fabs(*mjj) >= alpha * lambda * lambda)
00924                     {
00925                       ss=1;
00926                       pivrow = j;
00927                     }
00928                   else if (std::fabs(*(m.begin()+pivrow*(pivrow-1)/2+pivrow-1)) 
00929                                 >= alpha * sigma)
00930                     { ss=1; }
00931                   else
00932                     { ss=2; }
00933                 }
00934               if (pivrow == j)  // no permutation neccessary
00935                 {
00936                   piv[j-1] = pivrow;
00937                   if (*mjj == 0)
00938                     {
00939                       ifail=1;
00940                       return;
00941                     }
00942                   temp2 = *mjj = 1./ *mjj; // invert D(j,j)
00943                   
00944                   // update A(j+1:n, j+1,n)
00945                   for (i=j+1; i <= nrow; i++)
00946                     {
00947                       temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
00948                       ip = m.begin()+i*(i-1)/2+j;
00949                       for (k=j+1; k<=i; k++)
00950                         {
00951                           *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
00952                           if (std::fabs(*ip) <= epsilon)
00953                             { *ip=0; }
00954                           ip++;
00955                         }
00956                     }
00957                   // update L 
00958                   ip = m.begin() + (j+1)*j/2 + j-1; 
00959                   for (i=j+1; i <= nrow; ip += i++)
00960                   {
00961                     *ip *= temp2;
00962                   }
00963                 }
00964               else if (ss==1) // 1x1 pivot 
00965                 {
00966                   piv[j-1] = pivrow;
00967                   
00968                   // interchange rows and columns j and pivrow in
00969                   // submatrix (j:n,j:n)
00970                   ip = m.begin() + pivrow*(pivrow-1)/2 + j;
00971                   for (i=j+1; i < pivrow; i++, ip++)
00972                     {
00973                       temp1 = *(m.begin() + i*(i-1)/2 + j-1);
00974                       *(m.begin() + i*(i-1)/2 + j-1)= *ip;
00975                       *ip = temp1;
00976                     }
00977                   temp1 = *mjj;
00978                   *mjj = *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1);
00979                   *(m.begin()+pivrow*(pivrow-1)/2+pivrow-1) = temp1;
00980                   ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;
00981                   iq = ip + pivrow-j;
00982                   for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
00983                     {
00984                       temp1 = *iq;
00985                       *iq = *ip;
00986                       *ip = temp1;
00987                     }
00988                   
00989                   if (*mjj == 0)
00990                     {
00991                       ifail = 1;
00992                       return;
00993                     }
00994                   temp2 = *mjj = 1./ *mjj; // invert D(j,j)
00995                   
00996                   // update A(j+1:n, j+1:n)
00997                   for (i = j+1; i <= nrow; i++)
00998                     {
00999                       temp1 = *(m.begin() + i*(i-1)/2 + j-1) * temp2;
01000                       ip = m.begin()+i*(i-1)/2+j;
01001                       for (k=j+1; k<=i; k++)
01002                         {
01003                           *ip -= temp1 * *(m.begin() + k*(k-1)/2 + j-1);
01004                           if (std::fabs(*ip) <= epsilon)
01005                             { *ip=0; }
01006                           ip++;
01007                         }
01008                     }
01009                   // update L
01010                   ip = m.begin() + (j+1)*j/2 + j-1;
01011                   for (i=j+1; i<=nrow; ip += i++)
01012                   {
01013                     *ip *= temp2;
01014                   }
01015                 }
01016               else // ss=2, ie use a 2x2 pivot
01017                 {
01018                   piv[j-1] = -pivrow;
01019                   piv[j] = 0; // that means this is the second row of a 2x2 pivot
01020                   
01021                   if (j+1 != pivrow) 
01022                     {
01023                       // interchange rows and columns j+1 and pivrow in
01024                       // submatrix (j:n,j:n) 
01025                       ip = m.begin() + pivrow*(pivrow-1)/2 + j+1;
01026                       for (i=j+2; i < pivrow; i++, ip++)
01027                         {
01028                           temp1 = *(m.begin() + i*(i-1)/2 + j);
01029                           *(m.begin() + i*(i-1)/2 + j) = *ip;
01030                           *ip = temp1;
01031                         }
01032                       temp1 = *(mjj + j + 1);
01033                       *(mjj + j + 1) = 
01034                         *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
01035                       *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
01036                       temp1 = *(mjj + j);
01037                       *(mjj + j) = *(m.begin() + pivrow*(pivrow-1)/2 + j-1);
01038                       *(m.begin() + pivrow*(pivrow-1)/2 + j-1) = temp1;
01039                       ip = m.begin() + (pivrow+1)*pivrow/2 + j;
01040                       iq = ip + pivrow-(j+1);
01041                       for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
01042                         {
01043                           temp1 = *iq;
01044                           *iq = *ip;
01045                           *ip = temp1;
01046                         }
01047                     } 
01048                   // invert D(j:j+1,j:j+1)
01049                   temp2 = *mjj * *(mjj + j + 1) - *(mjj + j) * *(mjj + j); 
01050                   if (temp2 == 0)
01051                   {
01052                     G4cerr
01053                       << "G4ErrorSymMatrix::bunch_invert: error in pivot choice"
01054                       << G4endl;
01055                   }
01056                   temp2 = 1. / temp2;
01057 
01058                   // this quotient is guaranteed to exist by the choice 
01059                   // of the pivot
01060 
01061                   temp1 = *mjj;
01062                   *mjj = *(mjj + j + 1) * temp2;
01063                   *(mjj + j + 1) = temp1 * temp2;
01064                   *(mjj + j) = - *(mjj + j) * temp2;
01065                   
01066                   if (j < nrow-1) // otherwise do nothing
01067                     {
01068                       // update A(j+2:n, j+2:n)
01069                       for (i=j+2; i <= nrow ; i++)
01070                         {
01071                           ip = m.begin() + i*(i-1)/2 + j-1;
01072                           temp1 = *ip * *mjj + *(ip + 1) * *(mjj + j);
01073                           if (std::fabs(temp1 ) <= epsilon)
01074                             { temp1 = 0; }
01075                           temp2 = *ip * *(mjj + j) + *(ip + 1) * *(mjj + j + 1);
01076                           if (std::fabs(temp2 ) <= epsilon)
01077                             { temp2 = 0; }
01078                           for (k = j+2; k <= i ; k++)
01079                             {
01080                               ip = m.begin() + i*(i-1)/2 + k-1;
01081                               iq = m.begin() + k*(k-1)/2 + j-1;
01082                               *ip -= temp1 * *iq + temp2 * *(iq+1);
01083                               if (std::fabs(*ip) <= epsilon)
01084                                 { *ip = 0; }
01085                             }
01086                         }
01087                       // update L
01088                       for (i=j+2; i <= nrow ; i++)
01089                         {
01090                           ip = m.begin() + i*(i-1)/2 + j-1;
01091                           temp1 = *ip * *mjj + *(ip+1) * *(mjj + j);
01092                           if (std::fabs(temp1) <= epsilon)
01093                             { temp1 = 0; }
01094                           *(ip+1) = *ip * *(mjj + j) + *(ip+1) * *(mjj + j + 1);
01095                           if (std::fabs(*(ip+1)) <= epsilon)
01096                             { *(ip+1) = 0; }
01097                           *ip = temp1;
01098                         }
01099                     }
01100                 }
01101             }
01102   } // end of main loop over columns
01103 
01104   if (j == nrow) // the the last pivot is 1x1
01105   {
01106     mjj = m.begin() + j*(j-1)/2 + j-1;
01107     if (*mjj == 0)
01108     {
01109       ifail = 1;
01110       return;
01111     }
01112     else
01113     {
01114       *mjj = 1. / *mjj;
01115     }
01116   } // end of last pivot code
01117 
01118   // computing the inverse from the factorization
01119          
01120   for (j = nrow ; j >= 1 ; j -= ss) // loop over columns
01121   {
01122           mjj = m.begin() + j*(j-1)/2 + j-1;
01123           if (piv[j-1] > 0) // 1x1 pivot, compute column j of inverse
01124             {
01125               ss = 1; 
01126               if (j < nrow)
01127                 {
01128                   ip = m.begin() + (j+1)*j/2 + j-1;
01129                   for (i=0; i < nrow-j; ip += 1+j+i++)
01130                   {
01131                     x[i] = *ip;
01132                   }
01133                   for (i=j+1; i<=nrow ; i++)
01134                   {
01135                     temp2=0;
01136                     ip = m.begin() + i*(i-1)/2 + j;
01137                     for (k=0; k <= i-j-1; k++)
01138                       { temp2 += *ip++ * x[k]; }
01139                     for (ip += i-1; k < nrow-j; ip += 1+j+k++) 
01140                       { temp2 += *ip * x[k]; }
01141                     *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
01142                   }
01143                   temp2 = 0;
01144                   ip = m.begin() + (j+1)*j/2 + j-1;
01145                   for (k=0; k < nrow-j; ip += 1+j+k++)
01146                     { temp2 += x[k] * *ip; }
01147                   *mjj -= temp2;
01148                 }
01149             }
01150           else //2x2 pivot, compute columns j and j-1 of the inverse
01151             {
01152               if (piv[j-1] != 0)
01153                 { G4cerr << "error in piv" << piv[j-1] << G4endl; }
01154               ss=2; 
01155               if (j < nrow)
01156                 {
01157                   ip = m.begin() + (j+1)*j/2 + j-1;
01158                   for (i=0; i < nrow-j; ip += 1+j+i++)
01159                   {
01160                     x[i] = *ip;
01161                   }
01162                   for (i=j+1; i<=nrow ; i++)
01163                   {
01164                     temp2 = 0;
01165                     ip = m.begin() + i*(i-1)/2 + j;
01166                     for (k=0; k <= i-j-1; k++)
01167                       { temp2 += *ip++ * x[k]; }
01168                     for (ip += i-1; k < nrow-j; ip += 1+j+k++)
01169                       { temp2 += *ip * x[k]; }
01170                     *(m.begin()+ i*(i-1)/2 + j-1) = -temp2;
01171                   }    
01172                   temp2 = 0;
01173                   ip = m.begin() + (j+1)*j/2 + j-1;
01174                   for (k=0; k < nrow-j; ip += 1+j+k++)
01175                     { temp2 += x[k] * *ip; }
01176                   *mjj -= temp2;
01177                   temp2 = 0;
01178                   ip = m.begin() + (j+1)*j/2 + j-2;
01179                   for (i=j+1; i <= nrow; ip += i++)
01180                     { temp2 += *ip * *(ip+1); }
01181                   *(mjj-1) -= temp2;
01182                   ip = m.begin() + (j+1)*j/2 + j-2;
01183                   for (i=0; i < nrow-j; ip += 1+j+i++)
01184                   {
01185                     x[i] = *ip;
01186                   }
01187                   for (i=j+1; i <= nrow ; i++)
01188                   {
01189                     temp2 = 0;
01190                     ip = m.begin() + i*(i-1)/2 + j;
01191                     for (k=0; k <= i-j-1; k++)
01192                       { temp2 += *ip++ * x[k]; }
01193                     for (ip += i-1; k < nrow-j; ip += 1+j+k++)
01194                       { temp2 += *ip * x[k]; }
01195                     *(m.begin()+ i*(i-1)/2 + j-2)= -temp2;
01196                   }
01197                   temp2 = 0;
01198                   ip = m.begin() + (j+1)*j/2 + j-2;
01199                   for (k=0; k < nrow-j; ip += 1+j+k++)
01200                     { temp2 += x[k] * *ip; }
01201                   *(mjj-j) -= temp2;
01202                 }
01203             }  
01204           
01205           // interchange rows and columns j and piv[j-1] 
01206           // or rows and columns j and -piv[j-2]
01207           
01208           pivrow = (piv[j-1]==0)? -piv[j-2] : piv[j-1];
01209           ip = m.begin() + pivrow*(pivrow-1)/2 + j;
01210           for (i=j+1;i < pivrow; i++, ip++)
01211             {
01212               temp1 = *(m.begin() + i*(i-1)/2 + j-1);
01213               *(m.begin() + i*(i-1)/2 + j-1) = *ip;
01214               *ip = temp1;
01215             }
01216           temp1 = *mjj;
01217           *mjj = *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1);
01218           *(m.begin() + pivrow*(pivrow-1)/2 + pivrow-1) = temp1;
01219           if (ss==2)
01220             {
01221               temp1 = *(mjj-1);
01222               *(mjj-1) = *( m.begin() + pivrow*(pivrow-1)/2 + j-2);
01223               *( m.begin() + pivrow*(pivrow-1)/2 + j-2) = temp1;
01224             }
01225           
01226           ip = m.begin() + (pivrow+1)*pivrow/2 + j-1;  // &A(i,j)
01227           iq = ip + pivrow-j;
01228           for (i = pivrow+1; i <= nrow; ip += i, iq += i++)
01229             {
01230               temp1 = *iq;
01231               *iq = *ip;
01232               *ip = temp1;
01233             } 
01234   } // end of loop over columns (in computing inverse from factorization)
01235 
01236   return; // inversion successful
01237 }

void G4ErrorSymMatrix::invertCholesky5 ( G4int ifail  ) 

Definition at line 1911 of file G4ErrorSymMatrix.cc.

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

01912 {
01913 
01914   // Invert by 
01915   //
01916   // a) decomposing M = G*G^T with G lower triangular
01917   //    (if M is not positive definite this will fail, leaving this unchanged)
01918   // b) inverting G to form H
01919   // c) multiplying H^T * H to get M^-1.
01920   //
01921   // If the matrix is pos. def. it is inverted and 1 is returned.
01922   // If the matrix is not pos. def. it remains unaltered and 0 is returned.
01923 
01924   G4double h10;                           // below-diagonal elements of H
01925   G4double h20, h21;
01926   G4double h30, h31, h32;
01927   G4double h40, h41, h42, h43;
01928 
01929   G4double h00, h11, h22, h33, h44;       // 1/diagonal elements of G = 
01930                                           // diagonal elements of H
01931 
01932   G4double g10;                           // below-diagonal elements of G
01933   G4double g20, g21;
01934   G4double g30, g31, g32;
01935   G4double g40, g41, g42, g43;
01936 
01937   ifail = 1;  // We start by assuing we won't succeed...
01938 
01939   // Form G -- compute diagonal members of H directly rather than of G
01940   //-------
01941 
01942   // Scale first column by 1/sqrt(A00)
01943 
01944   h00 = m[A00]; 
01945   if (h00 <= 0) { return; }
01946   h00 = 1.0 / std::sqrt(h00);
01947 
01948   g10 = m[A10] * h00;
01949   g20 = m[A20] * h00;
01950   g30 = m[A30] * h00;
01951   g40 = m[A40] * h00;
01952 
01953   // Form G11 (actually, just h11)
01954 
01955   h11 = m[A11] - (g10 * g10);
01956   if (h11 <= 0) { return; }
01957   h11 = 1.0 / std::sqrt(h11);
01958 
01959   // Subtract inter-column column dot products from rest of column 1 and
01960   // scale to get column 1 of G 
01961 
01962   g21 = (m[A21] - (g10 * g20)) * h11;
01963   g31 = (m[A31] - (g10 * g30)) * h11;
01964   g41 = (m[A41] - (g10 * g40)) * h11;
01965 
01966   // Form G22 (actually, just h22)
01967 
01968   h22 = m[A22] - (g20 * g20) - (g21 * g21);
01969   if (h22 <= 0) { return; }
01970   h22 = 1.0 / std::sqrt(h22);
01971 
01972   // Subtract inter-column column dot products from rest of column 2 and
01973   // scale to get column 2 of G 
01974 
01975   g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
01976   g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
01977 
01978   // Form G33 (actually, just h33)
01979 
01980   h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
01981   if (h33 <= 0) { return; }
01982   h33 = 1.0 / std::sqrt(h33);
01983 
01984   // Subtract inter-column column dot product from A43 and scale to get G43
01985 
01986   g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
01987 
01988   // Finally form h44 - if this is possible inversion succeeds
01989 
01990   h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
01991   if (h44 <= 0) { return; }
01992   h44 = 1.0 / std::sqrt(h44);
01993 
01994   // Form H = 1/G -- diagonal members of H are already correct
01995   //-------------
01996 
01997   // The order here is dictated by speed considerations
01998 
01999   h43 = -h33 *  g43 * h44;
02000   h32 = -h22 *  g32 * h33;
02001   h42 = -h22 * (g32 * h43 + g42 * h44);
02002   h21 = -h11 *  g21 * h22;
02003   h31 = -h11 * (g21 * h32 + g31 * h33);
02004   h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
02005   h10 = -h00 *  g10 * h11;
02006   h20 = -h00 * (g10 * h21 + g20 * h22);
02007   h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
02008   h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
02009 
02010   // Change this to its inverse = H^T*H
02011   //------------------------------------
02012 
02013   m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40;
02014   m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41;
02015   m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41;
02016   m[A02] = h20 * h22 + h30 * h32 + h40 * h42;
02017   m[A12] = h21 * h22 + h31 * h32 + h41 * h42;
02018   m[A22] = h22 * h22 + h32 * h32 + h42 * h42;
02019   m[A03] = h30 * h33 + h40 * h43;
02020   m[A13] = h31 * h33 + h41 * h43;
02021   m[A23] = h32 * h33 + h42 * h43;
02022   m[A33] = h33 * h33 + h43 * h43;
02023   m[A04] = h40 * h44;
02024   m[A14] = h41 * h44;
02025   m[A24] = h42 * h44;
02026   m[A34] = h43 * h44;
02027   m[A44] = h44 * h44;
02028 
02029   ifail = 0;
02030   return;
02031 }

void G4ErrorSymMatrix::invertCholesky6 ( G4int ifail  ) 

Definition at line 2033 of file G4ErrorSymMatrix.cc.

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

02034 {
02035   // Invert by 
02036   //
02037   // a) decomposing M = G*G^T with G lower triangular
02038   //    (if M is not positive definite this will fail, leaving this unchanged)
02039   // b) inverting G to form H
02040   // c) multiplying H^T * H to get M^-1.
02041   //
02042   // If the matrix is pos. def. it is inverted and 1 is returned.
02043   // If the matrix is not pos. def. it remains unaltered and 0 is returned.
02044 
02045   G4double h10;                           // below-diagonal elements of H
02046   G4double h20, h21;
02047   G4double h30, h31, h32;
02048   G4double h40, h41, h42, h43;
02049   G4double h50, h51, h52, h53, h54;
02050 
02051   G4double h00, h11, h22, h33, h44, h55;  // 1/diagonal elements of G = 
02052                                           // diagonal elements of H
02053 
02054   G4double g10;                           // below-diagonal elements of G
02055   G4double g20, g21;
02056   G4double g30, g31, g32;
02057   G4double g40, g41, g42, g43;
02058   G4double g50, g51, g52, g53, g54;
02059 
02060   ifail = 1;  // We start by assuing we won't succeed...
02061 
02062   // Form G -- compute diagonal members of H directly rather than of G
02063   //-------
02064 
02065   // Scale first column by 1/sqrt(A00)
02066 
02067   h00 = m[A00]; 
02068   if (h00 <= 0) { return; }
02069   h00 = 1.0 / std::sqrt(h00);
02070 
02071   g10 = m[A10] * h00;
02072   g20 = m[A20] * h00;
02073   g30 = m[A30] * h00;
02074   g40 = m[A40] * h00;
02075   g50 = m[A50] * h00;
02076 
02077   // Form G11 (actually, just h11)
02078 
02079   h11 = m[A11] - (g10 * g10);
02080   if (h11 <= 0) { return; }
02081   h11 = 1.0 / std::sqrt(h11);
02082 
02083   // Subtract inter-column column dot products from rest of column 1 and
02084   // scale to get column 1 of G 
02085 
02086   g21 = (m[A21] - (g10 * g20)) * h11;
02087   g31 = (m[A31] - (g10 * g30)) * h11;
02088   g41 = (m[A41] - (g10 * g40)) * h11;
02089   g51 = (m[A51] - (g10 * g50)) * h11;
02090 
02091   // Form G22 (actually, just h22)
02092 
02093   h22 = m[A22] - (g20 * g20) - (g21 * g21);
02094   if (h22 <= 0) { return; }
02095   h22 = 1.0 / std::sqrt(h22);
02096 
02097   // Subtract inter-column column dot products from rest of column 2 and
02098   // scale to get column 2 of G 
02099 
02100   g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
02101   g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
02102   g52 = (m[A52] - (g20 * g50) - (g21 * g51)) * h22;
02103 
02104   // Form G33 (actually, just h33)
02105 
02106   h33 = m[A33] - (g30 * g30) - (g31 * g31) - (g32 * g32);
02107   if (h33 <= 0) { return; }
02108   h33 = 1.0 / std::sqrt(h33);
02109 
02110   // Subtract inter-column column dot products from rest of column 3 and
02111   // scale to get column 3 of G 
02112 
02113   g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
02114   g53 = (m[A53] - (g30 * g50) - (g31 * g51) - (g32 * g52)) * h33;
02115 
02116   // Form G44 (actually, just h44)
02117 
02118   h44 = m[A44] - (g40 * g40) - (g41 * g41) - (g42 * g42) - (g43 * g43);
02119   if (h44 <= 0) { return; }
02120   h44 = 1.0 / std::sqrt(h44);
02121 
02122   // Subtract inter-column column dot product from M54 and scale to get G54
02123 
02124   g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;
02125 
02126   // Finally form h55 - if this is possible inversion succeeds
02127 
02128   h55 = m[A55] - (g50*g50) - (g51*g51) - (g52*g52) - (g53*g53) - (g54*g54);
02129   if (h55 <= 0) { return; }
02130   h55 = 1.0 / std::sqrt(h55);
02131 
02132   // Form H = 1/G -- diagonal members of H are already correct
02133   //-------------
02134 
02135   // The order here is dictated by speed considerations
02136 
02137   h54 = -h44 *  g54 * h55;
02138   h43 = -h33 *  g43 * h44;
02139   h53 = -h33 * (g43 * h54 + g53 * h55);
02140   h32 = -h22 *  g32 * h33;
02141   h42 = -h22 * (g32 * h43 + g42 * h44);
02142   h52 = -h22 * (g32 * h53 + g42 * h54 + g52 * h55);
02143   h21 = -h11 *  g21 * h22;
02144   h31 = -h11 * (g21 * h32 + g31 * h33);
02145   h41 = -h11 * (g21 * h42 + g31 * h43 + g41 * h44);
02146   h51 = -h11 * (g21 * h52 + g31 * h53 + g41 * h54 + g51 * h55);
02147   h10 = -h00 *  g10 * h11;
02148   h20 = -h00 * (g10 * h21 + g20 * h22);
02149   h30 = -h00 * (g10 * h31 + g20 * h32 + g30 * h33);
02150   h40 = -h00 * (g10 * h41 + g20 * h42 + g30 * h43 + g40 * h44);
02151   h50 = -h00 * (g10 * h51 + g20 * h52 + g30 * h53 + g40 * h54 + g50 * h55);
02152 
02153   // Change this to its inverse = H^T*H
02154   //------------------------------------
02155 
02156   m[A00] = h00 * h00 + h10 * h10 + h20 * h20 + h30 * h30 + h40 * h40 + h50*h50;
02157   m[A01] = h10 * h11 + h20 * h21 + h30 * h31 + h40 * h41 + h50 * h51;
02158   m[A11] = h11 * h11 + h21 * h21 + h31 * h31 + h41 * h41 + h51 * h51;
02159   m[A02] = h20 * h22 + h30 * h32 + h40 * h42 + h50 * h52;
02160   m[A12] = h21 * h22 + h31 * h32 + h41 * h42 + h51 * h52;
02161   m[A22] = h22 * h22 + h32 * h32 + h42 * h42 + h52 * h52;
02162   m[A03] = h30 * h33 + h40 * h43 + h50 * h53;
02163   m[A13] = h31 * h33 + h41 * h43 + h51 * h53;
02164   m[A23] = h32 * h33 + h42 * h43 + h52 * h53;
02165   m[A33] = h33 * h33 + h43 * h43 + h53 * h53;
02166   m[A04] = h40 * h44 + h50 * h54;
02167   m[A14] = h41 * h44 + h51 * h54;
02168   m[A24] = h42 * h44 + h52 * h54;
02169   m[A34] = h43 * h44 + h53 * h54;
02170   m[A44] = h44 * h44 + h54 * h54;
02171   m[A05] = h50 * h55;
02172   m[A15] = h51 * h55;
02173   m[A25] = h52 * h55;
02174   m[A35] = h53 * h55;
02175   m[A45] = h54 * h55;
02176   m[A55] = h55 * h55;
02177 
02178   ifail = 0;
02179   return;
02180 }

void G4ErrorSymMatrix::invertHaywood4 ( G4int ifail  ) 

Definition at line 2260 of file G4ErrorSymMatrix.cc.

02261 {
02262   invert4(ifail); // For the 4x4 case, the method we use for invert is already
02263                   // the Haywood method.
02264 }

void G4ErrorSymMatrix::invertHaywood5 ( G4int ifail  ) 

Definition at line 1358 of file G4ErrorSymMatrix.cc.

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

01359 {
01360   ifail = 0;
01361 
01362   // Find all NECESSARY 2x2 dets:  (25 of them)
01363 
01364   G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
01365   G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
01366   G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
01367   G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
01368   G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
01369   G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
01370   G4double Det2_24_01 = m[A20]*m[A41] - m[A21]*m[A40];
01371   G4double Det2_24_02 = m[A20]*m[A42] - m[A22]*m[A40];
01372   G4double Det2_24_03 = m[A20]*m[A43] - m[A23]*m[A40];
01373   G4double Det2_24_04 = m[A20]*m[A44] - m[A24]*m[A40];
01374   G4double Det2_24_12 = m[A21]*m[A42] - m[A22]*m[A41];
01375   G4double Det2_24_13 = m[A21]*m[A43] - m[A23]*m[A41];
01376   G4double Det2_24_14 = m[A21]*m[A44] - m[A24]*m[A41];
01377   G4double Det2_24_23 = m[A22]*m[A43] - m[A23]*m[A42];
01378   G4double Det2_24_24 = m[A22]*m[A44] - m[A24]*m[A42];
01379   G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
01380   G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
01381   G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
01382   G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
01383   G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
01384   G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
01385   G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
01386   G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
01387   G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
01388   G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
01389 
01390   // Find all NECESSARY 3x3 dets:   (30 of them)
01391 
01392   G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02 
01393                                 + m[A12]*Det2_23_01;
01394   G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03 
01395                                 + m[A13]*Det2_23_01;
01396   G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03 
01397                                 + m[A13]*Det2_23_02;
01398   G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13 
01399                                 + m[A13]*Det2_23_12;
01400   G4double Det3_124_012 = m[A10]*Det2_24_12 - m[A11]*Det2_24_02 
01401                                 + m[A12]*Det2_24_01;
01402   G4double Det3_124_013 = m[A10]*Det2_24_13 - m[A11]*Det2_24_03 
01403                                 + m[A13]*Det2_24_01;
01404   G4double Det3_124_014 = m[A10]*Det2_24_14 - m[A11]*Det2_24_04 
01405                                 + m[A14]*Det2_24_01;
01406   G4double Det3_124_023 = m[A10]*Det2_24_23 - m[A12]*Det2_24_03 
01407                                 + m[A13]*Det2_24_02;
01408   G4double Det3_124_024 = m[A10]*Det2_24_24 - m[A12]*Det2_24_04 
01409                                 + m[A14]*Det2_24_02;
01410   G4double Det3_124_123 = m[A11]*Det2_24_23 - m[A12]*Det2_24_13 
01411                                 + m[A13]*Det2_24_12;
01412   G4double Det3_124_124 = m[A11]*Det2_24_24 - m[A12]*Det2_24_14 
01413                                 + m[A14]*Det2_24_12;
01414   G4double Det3_134_012 = m[A10]*Det2_34_12 - m[A11]*Det2_34_02 
01415                                 + m[A12]*Det2_34_01;
01416   G4double Det3_134_013 = m[A10]*Det2_34_13 - m[A11]*Det2_34_03 
01417                                 + m[A13]*Det2_34_01;
01418   G4double Det3_134_014 = m[A10]*Det2_34_14 - m[A11]*Det2_34_04 
01419                                 + m[A14]*Det2_34_01;
01420   G4double Det3_134_023 = m[A10]*Det2_34_23 - m[A12]*Det2_34_03 
01421                                 + m[A13]*Det2_34_02;
01422   G4double Det3_134_024 = m[A10]*Det2_34_24 - m[A12]*Det2_34_04 
01423                                 + m[A14]*Det2_34_02;
01424   G4double Det3_134_034 = m[A10]*Det2_34_34 - m[A13]*Det2_34_04 
01425                                 + m[A14]*Det2_34_03;
01426   G4double Det3_134_123 = m[A11]*Det2_34_23 - m[A12]*Det2_34_13 
01427                                 + m[A13]*Det2_34_12;
01428   G4double Det3_134_124 = m[A11]*Det2_34_24 - m[A12]*Det2_34_14 
01429                                 + m[A14]*Det2_34_12;
01430   G4double Det3_134_134 = m[A11]*Det2_34_34 - m[A13]*Det2_34_14 
01431                                 + m[A14]*Det2_34_13;
01432   G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02 
01433                                 + m[A22]*Det2_34_01;
01434   G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03 
01435                                 + m[A23]*Det2_34_01;
01436   G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04 
01437                                 + m[A24]*Det2_34_01;
01438   G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03 
01439                                 + m[A23]*Det2_34_02;
01440   G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04 
01441                                 + m[A24]*Det2_34_02;
01442   G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04 
01443                                 + m[A24]*Det2_34_03;
01444   G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13 
01445                                 + m[A23]*Det2_34_12;
01446   G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14 
01447                                 + m[A24]*Det2_34_12;
01448   G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14 
01449                                 + m[A24]*Det2_34_13;
01450   G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24 
01451                                 + m[A24]*Det2_34_23;
01452 
01453   // Find all NECESSARY 4x4 dets:   (15 of them)
01454 
01455   G4double Det4_0123_0123 = m[A00]*Det3_123_123 - m[A01]*Det3_123_023 
01456                                 + m[A02]*Det3_123_013 - m[A03]*Det3_123_012;
01457   G4double Det4_0124_0123 = m[A00]*Det3_124_123 - m[A01]*Det3_124_023 
01458                                 + m[A02]*Det3_124_013 - m[A03]*Det3_124_012;
01459   G4double Det4_0124_0124 = m[A00]*Det3_124_124 - m[A01]*Det3_124_024 
01460                                 + m[A02]*Det3_124_014 - m[A04]*Det3_124_012;
01461   G4double Det4_0134_0123 = m[A00]*Det3_134_123 - m[A01]*Det3_134_023 
01462                                 + m[A02]*Det3_134_013 - m[A03]*Det3_134_012;
01463   G4double Det4_0134_0124 = m[A00]*Det3_134_124 - m[A01]*Det3_134_024 
01464                                 + m[A02]*Det3_134_014 - m[A04]*Det3_134_012;
01465   G4double Det4_0134_0134 = m[A00]*Det3_134_134 - m[A01]*Det3_134_034 
01466                                 + m[A03]*Det3_134_014 - m[A04]*Det3_134_013;
01467   G4double Det4_0234_0123 = m[A00]*Det3_234_123 - m[A01]*Det3_234_023 
01468                                 + m[A02]*Det3_234_013 - m[A03]*Det3_234_012;
01469   G4double Det4_0234_0124 = m[A00]*Det3_234_124 - m[A01]*Det3_234_024 
01470                                 + m[A02]*Det3_234_014 - m[A04]*Det3_234_012;
01471   G4double Det4_0234_0134 = m[A00]*Det3_234_134 - m[A01]*Det3_234_034 
01472                                 + m[A03]*Det3_234_014 - m[A04]*Det3_234_013;
01473   G4double Det4_0234_0234 = m[A00]*Det3_234_234 - m[A02]*Det3_234_034 
01474                                 + m[A03]*Det3_234_024 - m[A04]*Det3_234_023;
01475   G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023 
01476                                 + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
01477   G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024 
01478                                 + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
01479   G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034 
01480                                 + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
01481   G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034 
01482                                 + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
01483   G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134 
01484                                 + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
01485 
01486   // Find the 5x5 det:
01487 
01488   G4double det =  m[A00]*Det4_1234_1234 
01489                 - m[A01]*Det4_1234_0234 
01490                 + m[A02]*Det4_1234_0134 
01491                 - m[A03]*Det4_1234_0124 
01492                 + m[A04]*Det4_1234_0123;
01493 
01494   if ( det == 0 )
01495   {  
01496     ifail = 1;
01497     return;
01498   } 
01499 
01500   G4double oneOverDet = 1.0/det;
01501   G4double mn1OverDet = - oneOverDet;
01502 
01503   m[A00] =  Det4_1234_1234 * oneOverDet;
01504   m[A01] =  Det4_1234_0234 * mn1OverDet;
01505   m[A02] =  Det4_1234_0134 * oneOverDet;
01506   m[A03] =  Det4_1234_0124 * mn1OverDet;
01507   m[A04] =  Det4_1234_0123 * oneOverDet;
01508 
01509   m[A11] =  Det4_0234_0234 * oneOverDet;
01510   m[A12] =  Det4_0234_0134 * mn1OverDet;
01511   m[A13] =  Det4_0234_0124 * oneOverDet;
01512   m[A14] =  Det4_0234_0123 * mn1OverDet;
01513 
01514   m[A22] =  Det4_0134_0134 * oneOverDet;
01515   m[A23] =  Det4_0134_0124 * mn1OverDet;
01516   m[A24] =  Det4_0134_0123 * oneOverDet;
01517 
01518   m[A33] =  Det4_0124_0124 * oneOverDet;
01519   m[A34] =  Det4_0124_0123 * mn1OverDet;
01520 
01521   m[A44] =  Det4_0123_0123 * oneOverDet;
01522 
01523   return;
01524 }

void G4ErrorSymMatrix::invertHaywood6 ( G4int ifail  ) 

Definition at line 1526 of file G4ErrorSymMatrix.cc.

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

01527 {
01528   ifail = 0;
01529 
01530   // Find all NECESSARY 2x2 dets:  (39 of them)
01531 
01532   G4double Det2_34_01 = m[A30]*m[A41] - m[A31]*m[A40];
01533   G4double Det2_34_02 = m[A30]*m[A42] - m[A32]*m[A40];
01534   G4double Det2_34_03 = m[A30]*m[A43] - m[A33]*m[A40];
01535   G4double Det2_34_04 = m[A30]*m[A44] - m[A34]*m[A40];
01536   G4double Det2_34_12 = m[A31]*m[A42] - m[A32]*m[A41];
01537   G4double Det2_34_13 = m[A31]*m[A43] - m[A33]*m[A41];
01538   G4double Det2_34_14 = m[A31]*m[A44] - m[A34]*m[A41];
01539   G4double Det2_34_23 = m[A32]*m[A43] - m[A33]*m[A42];
01540   G4double Det2_34_24 = m[A32]*m[A44] - m[A34]*m[A42];
01541   G4double Det2_34_34 = m[A33]*m[A44] - m[A34]*m[A43];
01542   G4double Det2_35_01 = m[A30]*m[A51] - m[A31]*m[A50];
01543   G4double Det2_35_02 = m[A30]*m[A52] - m[A32]*m[A50];
01544   G4double Det2_35_03 = m[A30]*m[A53] - m[A33]*m[A50];
01545   G4double Det2_35_04 = m[A30]*m[A54] - m[A34]*m[A50];
01546   G4double Det2_35_05 = m[A30]*m[A55] - m[A35]*m[A50];
01547   G4double Det2_35_12 = m[A31]*m[A52] - m[A32]*m[A51];
01548   G4double Det2_35_13 = m[A31]*m[A53] - m[A33]*m[A51];
01549   G4double Det2_35_14 = m[A31]*m[A54] - m[A34]*m[A51];
01550   G4double Det2_35_15 = m[A31]*m[A55] - m[A35]*m[A51];
01551   G4double Det2_35_23 = m[A32]*m[A53] - m[A33]*m[A52];
01552   G4double Det2_35_24 = m[A32]*m[A54] - m[A34]*m[A52];
01553   G4double Det2_35_25 = m[A32]*m[A55] - m[A35]*m[A52];
01554   G4double Det2_35_34 = m[A33]*m[A54] - m[A34]*m[A53];
01555   G4double Det2_35_35 = m[A33]*m[A55] - m[A35]*m[A53];
01556   G4double Det2_45_01 = m[A40]*m[A51] - m[A41]*m[A50];
01557   G4double Det2_45_02 = m[A40]*m[A52] - m[A42]*m[A50];
01558   G4double Det2_45_03 = m[A40]*m[A53] - m[A43]*m[A50];
01559   G4double Det2_45_04 = m[A40]*m[A54] - m[A44]*m[A50];
01560   G4double Det2_45_05 = m[A40]*m[A55] - m[A45]*m[A50];
01561   G4double Det2_45_12 = m[A41]*m[A52] - m[A42]*m[A51];
01562   G4double Det2_45_13 = m[A41]*m[A53] - m[A43]*m[A51];
01563   G4double Det2_45_14 = m[A41]*m[A54] - m[A44]*m[A51];
01564   G4double Det2_45_15 = m[A41]*m[A55] - m[A45]*m[A51];
01565   G4double Det2_45_23 = m[A42]*m[A53] - m[A43]*m[A52];
01566   G4double Det2_45_24 = m[A42]*m[A54] - m[A44]*m[A52];
01567   G4double Det2_45_25 = m[A42]*m[A55] - m[A45]*m[A52];
01568   G4double Det2_45_34 = m[A43]*m[A54] - m[A44]*m[A53];
01569   G4double Det2_45_35 = m[A43]*m[A55] - m[A45]*m[A53];
01570   G4double Det2_45_45 = m[A44]*m[A55] - m[A45]*m[A54];
01571 
01572   // Find all NECESSARY 3x3 dets:  (65 of them)
01573 
01574   G4double Det3_234_012 = m[A20]*Det2_34_12 - m[A21]*Det2_34_02 
01575                                                 + m[A22]*Det2_34_01;
01576   G4double Det3_234_013 = m[A20]*Det2_34_13 - m[A21]*Det2_34_03 
01577                                                 + m[A23]*Det2_34_01;
01578   G4double Det3_234_014 = m[A20]*Det2_34_14 - m[A21]*Det2_34_04 
01579                                                 + m[A24]*Det2_34_01;
01580   G4double Det3_234_023 = m[A20]*Det2_34_23 - m[A22]*Det2_34_03 
01581                                                 + m[A23]*Det2_34_02;
01582   G4double Det3_234_024 = m[A20]*Det2_34_24 - m[A22]*Det2_34_04 
01583                                                 + m[A24]*Det2_34_02;
01584   G4double Det3_234_034 = m[A20]*Det2_34_34 - m[A23]*Det2_34_04 
01585                                                 + m[A24]*Det2_34_03;
01586   G4double Det3_234_123 = m[A21]*Det2_34_23 - m[A22]*Det2_34_13 
01587                                                 + m[A23]*Det2_34_12;
01588   G4double Det3_234_124 = m[A21]*Det2_34_24 - m[A22]*Det2_34_14 
01589                                                 + m[A24]*Det2_34_12;
01590   G4double Det3_234_134 = m[A21]*Det2_34_34 - m[A23]*Det2_34_14 
01591                                                 + m[A24]*Det2_34_13;
01592   G4double Det3_234_234 = m[A22]*Det2_34_34 - m[A23]*Det2_34_24 
01593                                                 + m[A24]*Det2_34_23;
01594   G4double Det3_235_012 = m[A20]*Det2_35_12 - m[A21]*Det2_35_02 
01595                                                 + m[A22]*Det2_35_01;
01596   G4double Det3_235_013 = m[A20]*Det2_35_13 - m[A21]*Det2_35_03 
01597                                                 + m[A23]*Det2_35_01;
01598   G4double Det3_235_014 = m[A20]*Det2_35_14 - m[A21]*Det2_35_04 
01599                                                 + m[A24]*Det2_35_01;
01600   G4double Det3_235_015 = m[A20]*Det2_35_15 - m[A21]*Det2_35_05 
01601                                                 + m[A25]*Det2_35_01;
01602   G4double Det3_235_023 = m[A20]*Det2_35_23 - m[A22]*Det2_35_03 
01603                                                 + m[A23]*Det2_35_02;
01604   G4double Det3_235_024 = m[A20]*Det2_35_24 - m[A22]*Det2_35_04 
01605                                                 + m[A24]*Det2_35_02;
01606   G4double Det3_235_025 = m[A20]*Det2_35_25 - m[A22]*Det2_35_05 
01607                                                 + m[A25]*Det2_35_02;
01608   G4double Det3_235_034 = m[A20]*Det2_35_34 - m[A23]*Det2_35_04 
01609                                                 + m[A24]*Det2_35_03;
01610   G4double Det3_235_035 = m[A20]*Det2_35_35 - m[A23]*Det2_35_05 
01611                                                 + m[A25]*Det2_35_03;
01612   G4double Det3_235_123 = m[A21]*Det2_35_23 - m[A22]*Det2_35_13 
01613                                                 + m[A23]*Det2_35_12;
01614   G4double Det3_235_124 = m[A21]*Det2_35_24 - m[A22]*Det2_35_14 
01615                                                 + m[A24]*Det2_35_12;
01616   G4double Det3_235_125 = m[A21]*Det2_35_25 - m[A22]*Det2_35_15 
01617                                                 + m[A25]*Det2_35_12;
01618   G4double Det3_235_134 = m[A21]*Det2_35_34 - m[A23]*Det2_35_14 
01619                                                 + m[A24]*Det2_35_13;
01620   G4double Det3_235_135 = m[A21]*Det2_35_35 - m[A23]*Det2_35_15 
01621                                                 + m[A25]*Det2_35_13;
01622   G4double Det3_235_234 = m[A22]*Det2_35_34 - m[A23]*Det2_35_24 
01623                                                 + m[A24]*Det2_35_23;
01624   G4double Det3_235_235 = m[A22]*Det2_35_35 - m[A23]*Det2_35_25 
01625                                                 + m[A25]*Det2_35_23;
01626   G4double Det3_245_012 = m[A20]*Det2_45_12 - m[A21]*Det2_45_02 
01627                                                 + m[A22]*Det2_45_01;
01628   G4double Det3_245_013 = m[A20]*Det2_45_13 - m[A21]*Det2_45_03 
01629                                                 + m[A23]*Det2_45_01;
01630   G4double Det3_245_014 = m[A20]*Det2_45_14 - m[A21]*Det2_45_04 
01631                                                 + m[A24]*Det2_45_01;
01632   G4double Det3_245_015 = m[A20]*Det2_45_15 - m[A21]*Det2_45_05 
01633                                                 + m[A25]*Det2_45_01;
01634   G4double Det3_245_023 = m[A20]*Det2_45_23 - m[A22]*Det2_45_03 
01635                                                 + m[A23]*Det2_45_02;
01636   G4double Det3_245_024 = m[A20]*Det2_45_24 - m[A22]*Det2_45_04 
01637                                                 + m[A24]*Det2_45_02;
01638   G4double Det3_245_025 = m[A20]*Det2_45_25 - m[A22]*Det2_45_05 
01639                                                 + m[A25]*Det2_45_02;
01640   G4double Det3_245_034 = m[A20]*Det2_45_34 - m[A23]*Det2_45_04 
01641                                                 + m[A24]*Det2_45_03;
01642   G4double Det3_245_035 = m[A20]*Det2_45_35 - m[A23]*Det2_45_05 
01643                                                 + m[A25]*Det2_45_03;
01644   G4double Det3_245_045 = m[A20]*Det2_45_45 - m[A24]*Det2_45_05 
01645                                                 + m[A25]*Det2_45_04;
01646   G4double Det3_245_123 = m[A21]*Det2_45_23 - m[A22]*Det2_45_13 
01647                                                 + m[A23]*Det2_45_12;
01648   G4double Det3_245_124 = m[A21]*Det2_45_24 - m[A22]*Det2_45_14 
01649                                                 + m[A24]*Det2_45_12;
01650   G4double Det3_245_125 = m[A21]*Det2_45_25 - m[A22]*Det2_45_15 
01651                                                 + m[A25]*Det2_45_12;
01652   G4double Det3_245_134 = m[A21]*Det2_45_34 - m[A23]*Det2_45_14 
01653                                                 + m[A24]*Det2_45_13;
01654   G4double Det3_245_135 = m[A21]*Det2_45_35 - m[A23]*Det2_45_15 
01655                                                 + m[A25]*Det2_45_13;
01656   G4double Det3_245_145 = m[A21]*Det2_45_45 - m[A24]*Det2_45_15 
01657                                                 + m[A25]*Det2_45_14;
01658   G4double Det3_245_234 = m[A22]*Det2_45_34 - m[A23]*Det2_45_24 
01659                                                 + m[A24]*Det2_45_23;
01660   G4double Det3_245_235 = m[A22]*Det2_45_35 - m[A23]*Det2_45_25 
01661                                                 + m[A25]*Det2_45_23;
01662   G4double Det3_245_245 = m[A22]*Det2_45_45 - m[A24]*Det2_45_25 
01663                                                 + m[A25]*Det2_45_24;
01664   G4double Det3_345_012 = m[A30]*Det2_45_12 - m[A31]*Det2_45_02 
01665                                                 + m[A32]*Det2_45_01;
01666   G4double Det3_345_013 = m[A30]*Det2_45_13 - m[A31]*Det2_45_03 
01667                                                 + m[A33]*Det2_45_01;
01668   G4double Det3_345_014 = m[A30]*Det2_45_14 - m[A31]*Det2_45_04 
01669                                                 + m[A34]*Det2_45_01;
01670   G4double Det3_345_015 = m[A30]*Det2_45_15 - m[A31]*Det2_45_05 
01671                                                 + m[A35]*Det2_45_01;
01672   G4double Det3_345_023 = m[A30]*Det2_45_23 - m[A32]*Det2_45_03 
01673                                                 + m[A33]*Det2_45_02;
01674   G4double Det3_345_024 = m[A30]*Det2_45_24 - m[A32]*Det2_45_04 
01675                                                 + m[A34]*Det2_45_02;
01676   G4double Det3_345_025 = m[A30]*Det2_45_25 - m[A32]*Det2_45_05 
01677                                                 + m[A35]*Det2_45_02;
01678   G4double Det3_345_034 = m[A30]*Det2_45_34 - m[A33]*Det2_45_04 
01679                                                 + m[A34]*Det2_45_03;
01680   G4double Det3_345_035 = m[A30]*Det2_45_35 - m[A33]*Det2_45_05 
01681                                                 + m[A35]*Det2_45_03;
01682   G4double Det3_345_045 = m[A30]*Det2_45_45 - m[A34]*Det2_45_05 
01683                                                 + m[A35]*Det2_45_04;
01684   G4double Det3_345_123 = m[A31]*Det2_45_23 - m[A32]*Det2_45_13 
01685                                                 + m[A33]*Det2_45_12;
01686   G4double Det3_345_124 = m[A31]*Det2_45_24 - m[A32]*Det2_45_14 
01687                                                 + m[A34]*Det2_45_12;
01688   G4double Det3_345_125 = m[A31]*Det2_45_25 - m[A32]*Det2_45_15 
01689                                                 + m[A35]*Det2_45_12;
01690   G4double Det3_345_134 = m[A31]*Det2_45_34 - m[A33]*Det2_45_14 
01691                                                 + m[A34]*Det2_45_13;
01692   G4double Det3_345_135 = m[A31]*Det2_45_35 - m[A33]*Det2_45_15 
01693                                                 + m[A35]*Det2_45_13;
01694   G4double Det3_345_145 = m[A31]*Det2_45_45 - m[A34]*Det2_45_15 
01695                                                 + m[A35]*Det2_45_14;
01696   G4double Det3_345_234 = m[A32]*Det2_45_34 - m[A33]*Det2_45_24 
01697                                                 + m[A34]*Det2_45_23;
01698   G4double Det3_345_235 = m[A32]*Det2_45_35 - m[A33]*Det2_45_25 
01699                                                 + m[A35]*Det2_45_23;
01700   G4double Det3_345_245 = m[A32]*Det2_45_45 - m[A34]*Det2_45_25 
01701                                                 + m[A35]*Det2_45_24;
01702   G4double Det3_345_345 = m[A33]*Det2_45_45 - m[A34]*Det2_45_35 
01703                                                 + m[A35]*Det2_45_34;
01704 
01705   // Find all NECESSARY 4x4 dets:  (55 of them)
01706 
01707   G4double Det4_1234_0123 = m[A10]*Det3_234_123 - m[A11]*Det3_234_023 
01708                         + m[A12]*Det3_234_013 - m[A13]*Det3_234_012;
01709   G4double Det4_1234_0124 = m[A10]*Det3_234_124 - m[A11]*Det3_234_024 
01710                         + m[A12]*Det3_234_014 - m[A14]*Det3_234_012;
01711   G4double Det4_1234_0134 = m[A10]*Det3_234_134 - m[A11]*Det3_234_034 
01712                         + m[A13]*Det3_234_014 - m[A14]*Det3_234_013;
01713   G4double Det4_1234_0234 = m[A10]*Det3_234_234 - m[A12]*Det3_234_034 
01714                         + m[A13]*Det3_234_024 - m[A14]*Det3_234_023;
01715   G4double Det4_1234_1234 = m[A11]*Det3_234_234 - m[A12]*Det3_234_134 
01716                         + m[A13]*Det3_234_124 - m[A14]*Det3_234_123;
01717   G4double Det4_1235_0123 = m[A10]*Det3_235_123 - m[A11]*Det3_235_023 
01718                         + m[A12]*Det3_235_013 - m[A13]*Det3_235_012;
01719   G4double Det4_1235_0124 = m[A10]*Det3_235_124 - m[A11]*Det3_235_024 
01720                         + m[A12]*Det3_235_014 - m[A14]*Det3_235_012;
01721   G4double Det4_1235_0125 = m[A10]*Det3_235_125 - m[A11]*Det3_235_025 
01722                         + m[A12]*Det3_235_015 - m[A15]*Det3_235_012;
01723   G4double Det4_1235_0134 = m[A10]*Det3_235_134 - m[A11]*Det3_235_034 
01724                         + m[A13]*Det3_235_014 - m[A14]*Det3_235_013;
01725   G4double Det4_1235_0135 = m[A10]*Det3_235_135 - m[A11]*Det3_235_035 
01726                         + m[A13]*Det3_235_015 - m[A15]*Det3_235_013;
01727   G4double Det4_1235_0234 = m[A10]*Det3_235_234 - m[A12]*Det3_235_034 
01728                         + m[A13]*Det3_235_024 - m[A14]*Det3_235_023;
01729   G4double Det4_1235_0235 = m[A10]*Det3_235_235 - m[A12]*Det3_235_035 
01730                         + m[A13]*Det3_235_025 - m[A15]*Det3_235_023;
01731   G4double Det4_1235_1234 = m[A11]*Det3_235_234 - m[A12]*Det3_235_134 
01732                         + m[A13]*Det3_235_124 - m[A14]*Det3_235_123;
01733   G4double Det4_1235_1235 = m[A11]*Det3_235_235 - m[A12]*Det3_235_135 
01734                         + m[A13]*Det3_235_125 - m[A15]*Det3_235_123;
01735   G4double Det4_1245_0123 = m[A10]*Det3_245_123 - m[A11]*Det3_245_023 
01736                         + m[A12]*Det3_245_013 - m[A13]*Det3_245_012;
01737   G4double Det4_1245_0124 = m[A10]*Det3_245_124 - m[A11]*Det3_245_024 
01738                         + m[A12]*Det3_245_014 - m[A14]*Det3_245_012;
01739   G4double Det4_1245_0125 = m[A10]*Det3_245_125 - m[A11]*Det3_245_025 
01740                         + m[A12]*Det3_245_015 - m[A15]*Det3_245_012;
01741   G4double Det4_1245_0134 = m[A10]*Det3_245_134 - m[A11]*Det3_245_034 
01742                         + m[A13]*Det3_245_014 - m[A14]*Det3_245_013;
01743   G4double Det4_1245_0135 = m[A10]*Det3_245_135 - m[A11]*Det3_245_035 
01744                         + m[A13]*Det3_245_015 - m[A15]*Det3_245_013;
01745   G4double Det4_1245_0145 = m[A10]*Det3_245_145 - m[A11]*Det3_245_045 
01746                         + m[A14]*Det3_245_015 - m[A15]*Det3_245_014;
01747   G4double Det4_1245_0234 = m[A10]*Det3_245_234 - m[A12]*Det3_245_034 
01748                         + m[A13]*Det3_245_024 - m[A14]*Det3_245_023;
01749   G4double Det4_1245_0235 = m[A10]*Det3_245_235 - m[A12]*Det3_245_035 
01750                         + m[A13]*Det3_245_025 - m[A15]*Det3_245_023;
01751   G4double Det4_1245_0245 = m[A10]*Det3_245_245 - m[A12]*Det3_245_045 
01752                         + m[A14]*Det3_245_025 - m[A15]*Det3_245_024;
01753   G4double Det4_1245_1234 = m[A11]*Det3_245_234 - m[A12]*Det3_245_134 
01754                         + m[A13]*Det3_245_124 - m[A14]*Det3_245_123;
01755   G4double Det4_1245_1235 = m[A11]*Det3_245_235 - m[A12]*Det3_245_135 
01756                         + m[A13]*Det3_245_125 - m[A15]*Det3_245_123;
01757   G4double Det4_1245_1245 = m[A11]*Det3_245_245 - m[A12]*Det3_245_145 
01758                         + m[A14]*Det3_245_125 - m[A15]*Det3_245_124;
01759   G4double Det4_1345_0123 = m[A10]*Det3_345_123 - m[A11]*Det3_345_023 
01760                         + m[A12]*Det3_345_013 - m[A13]*Det3_345_012;
01761   G4double Det4_1345_0124 = m[A10]*Det3_345_124 - m[A11]*Det3_345_024 
01762                         + m[A12]*Det3_345_014 - m[A14]*Det3_345_012;
01763   G4double Det4_1345_0125 = m[A10]*Det3_345_125 - m[A11]*Det3_345_025 
01764                         + m[A12]*Det3_345_015 - m[A15]*Det3_345_012;
01765   G4double Det4_1345_0134 = m[A10]*Det3_345_134 - m[A11]*Det3_345_034 
01766                         + m[A13]*Det3_345_014 - m[A14]*Det3_345_013;
01767   G4double Det4_1345_0135 = m[A10]*Det3_345_135 - m[A11]*Det3_345_035 
01768                         + m[A13]*Det3_345_015 - m[A15]*Det3_345_013;
01769   G4double Det4_1345_0145 = m[A10]*Det3_345_145 - m[A11]*Det3_345_045 
01770                         + m[A14]*Det3_345_015 - m[A15]*Det3_345_014;
01771   G4double Det4_1345_0234 = m[A10]*Det3_345_234 - m[A12]*Det3_345_034 
01772                         + m[A13]*Det3_345_024 - m[A14]*Det3_345_023;
01773   G4double Det4_1345_0235 = m[A10]*Det3_345_235 - m[A12]*Det3_345_035 
01774                         + m[A13]*Det3_345_025 - m[A15]*Det3_345_023;
01775   G4double Det4_1345_0245 = m[A10]*Det3_345_245 - m[A12]*Det3_345_045 
01776                         + m[A14]*Det3_345_025 - m[A15]*Det3_345_024;
01777   G4double Det4_1345_0345 = m[A10]*Det3_345_345 - m[A13]*Det3_345_045 
01778                         + m[A14]*Det3_345_035 - m[A15]*Det3_345_034;
01779   G4double Det4_1345_1234 = m[A11]*Det3_345_234 - m[A12]*Det3_345_134 
01780                         + m[A13]*Det3_345_124 - m[A14]*Det3_345_123;
01781   G4double Det4_1345_1235 = m[A11]*Det3_345_235 - m[A12]*Det3_345_135 
01782                         + m[A13]*Det3_345_125 - m[A15]*Det3_345_123;
01783   G4double Det4_1345_1245 = m[A11]*Det3_345_245 - m[A12]*Det3_345_145 
01784                         + m[A14]*Det3_345_125 - m[A15]*Det3_345_124;
01785   G4double Det4_1345_1345 = m[A11]*Det3_345_345 - m[A13]*Det3_345_145 
01786                         + m[A14]*Det3_345_135 - m[A15]*Det3_345_134;
01787   G4double Det4_2345_0123 = m[A20]*Det3_345_123 - m[A21]*Det3_345_023 
01788                         + m[A22]*Det3_345_013 - m[A23]*Det3_345_012;
01789   G4double Det4_2345_0124 = m[A20]*Det3_345_124 - m[A21]*Det3_345_024 
01790                         + m[A22]*Det3_345_014 - m[A24]*Det3_345_012;
01791   G4double Det4_2345_0125 = m[A20]*Det3_345_125 - m[A21]*Det3_345_025 
01792                         + m[A22]*Det3_345_015 - m[A25]*Det3_345_012;
01793   G4double Det4_2345_0134 = m[A20]*Det3_345_134 - m[A21]*Det3_345_034 
01794                         + m[A23]*Det3_345_014 - m[A24]*Det3_345_013;
01795   G4double Det4_2345_0135 = m[A20]*Det3_345_135 - m[A21]*Det3_345_035 
01796                         + m[A23]*Det3_345_015 - m[A25]*Det3_345_013;
01797   G4double Det4_2345_0145 = m[A20]*Det3_345_145 - m[A21]*Det3_345_045 
01798                         + m[A24]*Det3_345_015 - m[A25]*Det3_345_014;
01799   G4double Det4_2345_0234 = m[A20]*Det3_345_234 - m[A22]*Det3_345_034 
01800                         + m[A23]*Det3_345_024 - m[A24]*Det3_345_023;
01801   G4double Det4_2345_0235 = m[A20]*Det3_345_235 - m[A22]*Det3_345_035 
01802                         + m[A23]*Det3_345_025 - m[A25]*Det3_345_023;
01803   G4double Det4_2345_0245 = m[A20]*Det3_345_245 - m[A22]*Det3_345_045 
01804                         + m[A24]*Det3_345_025 - m[A25]*Det3_345_024;
01805   G4double Det4_2345_0345 = m[A20]*Det3_345_345 - m[A23]*Det3_345_045 
01806                         + m[A24]*Det3_345_035 - m[A25]*Det3_345_034;
01807   G4double Det4_2345_1234 = m[A21]*Det3_345_234 - m[A22]*Det3_345_134 
01808                         + m[A23]*Det3_345_124 - m[A24]*Det3_345_123;
01809   G4double Det4_2345_1235 = m[A21]*Det3_345_235 - m[A22]*Det3_345_135 
01810                         + m[A23]*Det3_345_125 - m[A25]*Det3_345_123;
01811   G4double Det4_2345_1245 = m[A21]*Det3_345_245 - m[A22]*Det3_345_145 
01812                         + m[A24]*Det3_345_125 - m[A25]*Det3_345_124;
01813   G4double Det4_2345_1345 = m[A21]*Det3_345_345 - m[A23]*Det3_345_145 
01814                         + m[A24]*Det3_345_135 - m[A25]*Det3_345_134;
01815   G4double Det4_2345_2345 = m[A22]*Det3_345_345 - m[A23]*Det3_345_245 
01816                         + m[A24]*Det3_345_235 - m[A25]*Det3_345_234;
01817 
01818   // Find all NECESSARY 5x5 dets:  (19 of them)
01819 
01820   G4double Det5_01234_01234 = m[A00]*Det4_1234_1234 - m[A01]*Det4_1234_0234 
01821     + m[A02]*Det4_1234_0134 - m[A03]*Det4_1234_0124 + m[A04]*Det4_1234_0123;
01822   G4double Det5_01235_01234 = m[A00]*Det4_1235_1234 - m[A01]*Det4_1235_0234 
01823     + m[A02]*Det4_1235_0134 - m[A03]*Det4_1235_0124 + m[A04]*Det4_1235_0123;
01824   G4double Det5_01235_01235 = m[A00]*Det4_1235_1235 - m[A01]*Det4_1235_0235 
01825     + m[A02]*Det4_1235_0135 - m[A03]*Det4_1235_0125 + m[A05]*Det4_1235_0123;
01826   G4double Det5_01245_01234 = m[A00]*Det4_1245_1234 - m[A01]*Det4_1245_0234 
01827     + m[A02]*Det4_1245_0134 - m[A03]*Det4_1245_0124 + m[A04]*Det4_1245_0123;
01828   G4double Det5_01245_01235 = m[A00]*Det4_1245_1235 - m[A01]*Det4_1245_0235 
01829     + m[A02]*Det4_1245_0135 - m[A03]*Det4_1245_0125 + m[A05]*Det4_1245_0123;
01830   G4double Det5_01245_01245 = m[A00]*Det4_1245_1245 - m[A01]*Det4_1245_0245 
01831     + m[A02]*Det4_1245_0145 - m[A04]*Det4_1245_0125 + m[A05]*Det4_1245_0124;
01832   G4double Det5_01345_01234 = m[A00]*Det4_1345_1234 - m[A01]*Det4_1345_0234 
01833     + m[A02]*Det4_1345_0134 - m[A03]*Det4_1345_0124 + m[A04]*Det4_1345_0123;
01834   G4double Det5_01345_01235 = m[A00]*Det4_1345_1235 - m[A01]*Det4_1345_0235 
01835     + m[A02]*Det4_1345_0135 - m[A03]*Det4_1345_0125 + m[A05]*Det4_1345_0123;
01836   G4double Det5_01345_01245 = m[A00]*Det4_1345_1245 - m[A01]*Det4_1345_0245 
01837     + m[A02]*Det4_1345_0145 - m[A04]*Det4_1345_0125 + m[A05]*Det4_1345_0124;
01838   G4double Det5_01345_01345 = m[A00]*Det4_1345_1345 - m[A01]*Det4_1345_0345 
01839     + m[A03]*Det4_1345_0145 - m[A04]*Det4_1345_0135 + m[A05]*Det4_1345_0134;
01840   G4double Det5_02345_01234 = m[A00]*Det4_2345_1234 - m[A01]*Det4_2345_0234 
01841     + m[A02]*Det4_2345_0134 - m[A03]*Det4_2345_0124 + m[A04]*Det4_2345_0123;
01842   G4double Det5_02345_01235 = m[A00]*Det4_2345_1235 - m[A01]*Det4_2345_0235 
01843     + m[A02]*Det4_2345_0135 - m[A03]*Det4_2345_0125 + m[A05]*Det4_2345_0123;
01844   G4double Det5_02345_01245 = m[A00]*Det4_2345_1245 - m[A01]*Det4_2345_0245 
01845     + m[A02]*Det4_2345_0145 - m[A04]*Det4_2345_0125 + m[A05]*Det4_2345_0124;
01846   G4double Det5_02345_01345 = m[A00]*Det4_2345_1345 - m[A01]*Det4_2345_0345 
01847     + m[A03]*Det4_2345_0145 - m[A04]*Det4_2345_0135 + m[A05]*Det4_2345_0134;
01848   G4double Det5_02345_02345 = m[A00]*Det4_2345_2345 - m[A02]*Det4_2345_0345 
01849     + m[A03]*Det4_2345_0245 - m[A04]*Det4_2345_0235 + m[A05]*Det4_2345_0234;
01850   G4double Det5_12345_01234 = m[A10]*Det4_2345_1234 - m[A11]*Det4_2345_0234 
01851     + m[A12]*Det4_2345_0134 - m[A13]*Det4_2345_0124 + m[A14]*Det4_2345_0123;
01852   G4double Det5_12345_01235 = m[A10]*Det4_2345_1235 - m[A11]*Det4_2345_0235 
01853     + m[A12]*Det4_2345_0135 - m[A13]*Det4_2345_0125 + m[A15]*Det4_2345_0123;
01854   G4double Det5_12345_01245 = m[A10]*Det4_2345_1245 - m[A11]*Det4_2345_0245 
01855     + m[A12]*Det4_2345_0145 - m[A14]*Det4_2345_0125 + m[A15]*Det4_2345_0124;
01856   G4double Det5_12345_01345 = m[A10]*Det4_2345_1345 - m[A11]*Det4_2345_0345 
01857     + m[A13]*Det4_2345_0145 - m[A14]*Det4_2345_0135 + m[A15]*Det4_2345_0134;
01858   G4double Det5_12345_02345 = m[A10]*Det4_2345_2345 - m[A12]*Det4_2345_0345 
01859     + m[A13]*Det4_2345_0245 - m[A14]*Det4_2345_0235 + m[A15]*Det4_2345_0234;
01860   G4double Det5_12345_12345 = m[A11]*Det4_2345_2345 - m[A12]*Det4_2345_1345 
01861     + m[A13]*Det4_2345_1245 - m[A14]*Det4_2345_1235 + m[A15]*Det4_2345_1234;
01862 
01863   // Find the determinant 
01864 
01865   G4double det =  m[A00]*Det5_12345_12345 
01866                 - m[A01]*Det5_12345_02345 
01867                 + m[A02]*Det5_12345_01345 
01868                 - m[A03]*Det5_12345_01245 
01869                 + m[A04]*Det5_12345_01235 
01870                 - m[A05]*Det5_12345_01234;
01871 
01872   if ( det == 0 )
01873   {  
01874     ifail = 1;
01875     return;
01876   } 
01877 
01878   G4double oneOverDet = 1.0/det;
01879   G4double mn1OverDet = - oneOverDet;
01880 
01881   m[A00] =  Det5_12345_12345*oneOverDet;
01882   m[A01] =  Det5_12345_02345*mn1OverDet;
01883   m[A02] =  Det5_12345_01345*oneOverDet;
01884   m[A03] =  Det5_12345_01245*mn1OverDet;
01885   m[A04] =  Det5_12345_01235*oneOverDet;
01886   m[A05] =  Det5_12345_01234*mn1OverDet;
01887 
01888   m[A11] =  Det5_02345_02345*oneOverDet;
01889   m[A12] =  Det5_02345_01345*mn1OverDet;
01890   m[A13] =  Det5_02345_01245*oneOverDet;
01891   m[A14] =  Det5_02345_01235*mn1OverDet;
01892   m[A15] =  Det5_02345_01234*oneOverDet;
01893 
01894   m[A22] =  Det5_01345_01345*oneOverDet;
01895   m[A23] =  Det5_01345_01245*mn1OverDet;
01896   m[A24] =  Det5_01345_01235*oneOverDet;
01897   m[A25] =  Det5_01345_01234*mn1OverDet;
01898 
01899   m[A33] =  Det5_01245_01245*oneOverDet;
01900   m[A34] =  Det5_01245_01235*mn1OverDet;
01901   m[A35] =  Det5_01245_01234*oneOverDet;
01902 
01903   m[A44] =  Det5_01235_01235*oneOverDet;
01904   m[A45] =  Det5_01235_01234*mn1OverDet;
01905 
01906   m[A55] =  Det5_01234_01234*oneOverDet;
01907 
01908   return;
01909 }

G4int G4ErrorSymMatrix::num_col (  )  const [inline]

Definition at line 43 of file G4ErrorSymMatrix.icc.

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

00044 {
00045   return nrow;
00046 }

G4int G4ErrorSymMatrix::num_row (  )  const [inline]

Definition at line 38 of file G4ErrorSymMatrix.icc.

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

00039 {
00040   return nrow;
00041 }

G4int G4ErrorSymMatrix::num_size (  )  const [inline, protected]

Definition at line 48 of file G4ErrorSymMatrix.icc.

Referenced by operator-().

00049 {
00050   return size;
00051 }

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

Definition at line 472 of file G4ErrorSymMatrix.cc.

References SIMPLE_UOP.

00473 {
00474   SIMPLE_UOP(*=)
00475   return (*this);
00476 }

G4double & G4ErrorSymMatrix::operator() ( G4int  row,
G4int  col 
) [inline]

Definition at line 64 of file G4ErrorSymMatrix.icc.

References fast().

00065 {
00066   return (row>=col? fast(row,col) : fast(col,row));
00067 }

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

Definition at line 70 of file G4ErrorSymMatrix.icc.

References fast().

00071 {
00072   return (row>=col? fast(row,col) : fast(col,row));
00073 }

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

Definition at line 427 of file G4ErrorSymMatrix.cc.

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

00428 {
00429   CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=);
00430   SIMPLE_BOP(+=)
00431   return (*this);
00432 }

G4ErrorSymMatrix G4ErrorSymMatrix::operator- (  )  const

Definition at line 197 of file G4ErrorSymMatrix.cc.

References m, and num_size().

00198 {
00199   G4ErrorSymMatrix mat2(nrow);
00200   G4ErrorMatrixConstIter a=m.begin();
00201   G4ErrorMatrixIter b=mat2.m.begin();
00202   G4ErrorMatrixConstIter e=m.begin()+num_size();
00203   for(;a<e; a++, b++) { (*b) = -(*a); }
00204   return mat2;
00205 }

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

Definition at line 459 of file G4ErrorSymMatrix.cc.

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

00460 {
00461   CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
00462   SIMPLE_BOP(-=)
00463   return (*this);
00464 }

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

Definition at line 466 of file G4ErrorSymMatrix.cc.

References SIMPLE_UOP.

00467 {
00468   SIMPLE_UOP(/=)
00469   return (*this);
00470 }

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

Definition at line 509 of file G4ErrorSymMatrix.cc.

References m, nrow, and size.

00510 {
00511    if (&mat1 == this) { return *this; }
00512    if(mat1.nrow != nrow)
00513    {
00514       nrow = mat1.nrow;
00515       size = mat1.size;
00516       m.resize(size);
00517    }
00518    m = mat1.m;
00519    return (*this);
00520 }

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

Definition at line 93 of file G4ErrorSymMatrix.icc.

00094 {
00095   const G4ErrorSymMatrix_row_const b(*this,r);
00096   return b;
00097 }

G4ErrorSymMatrix::G4ErrorSymMatrix_row G4ErrorSymMatrix::operator[] ( G4int   )  [inline]

Definition at line 86 of file G4ErrorSymMatrix.icc.

00087 {
00088   G4ErrorSymMatrix_row b(*this,r);
00089   return b;
00090 }

G4ErrorSymMatrix G4ErrorSymMatrix::similarity ( const G4ErrorSymMatrix m1  )  const

Definition at line 620 of file G4ErrorSymMatrix.cc.

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

00621 {
00622   G4ErrorSymMatrix mret(mat1.num_row());
00623   G4ErrorMatrix temp = mat1*(*this);
00624   G4int n = mat1.num_col();
00625   G4ErrorMatrixIter mr = mret.m.begin();
00626   G4ErrorMatrixIter tempr1 = temp.m.begin();
00627   for(G4int r=1;r<=mret.num_row();r++)
00628   {
00629     G4ErrorMatrixConstIter m1c1 = mat1.m.begin();
00630     G4int c;
00631     for(c=1;c<=r;c++)
00632     {
00633       G4double tmp = 0.0;
00634       G4ErrorMatrixIter tempri = tempr1;
00635       G4ErrorMatrixConstIter m1ci = m1c1;
00636       G4int i;
00637       for(i=1;i<c;i++)
00638       {
00639         tmp+=(*(tempri++))*(*(m1ci++));
00640       }
00641       for(i=c;i<=mat1.num_col();i++)
00642       {
00643         tmp+=(*(tempri++))*(*(m1ci));
00644         m1ci += i;
00645       }
00646       *(mr++) = tmp;
00647       m1c1 += c;
00648     }
00649     tempr1 += n;
00650   }
00651   return mret;
00652 }

G4ErrorSymMatrix G4ErrorSymMatrix::similarity ( const G4ErrorMatrix m1  )  const

Definition at line 589 of file G4ErrorSymMatrix.cc.

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

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

00590 {
00591   G4ErrorSymMatrix mret(mat1.num_row());
00592   G4ErrorMatrix temp = mat1*(*this);
00593 
00594 // If mat1*(*this) has correct dimensions, then so will the mat1.T multiplication.
00595 // So there is no need to check dimensions again.
00596 
00597   G4int n = mat1.num_col();
00598   G4ErrorMatrixIter mr = mret.m.begin();
00599   G4ErrorMatrixIter tempr1 = temp.m.begin();
00600   for(G4int r=1;r<=mret.num_row();r++)
00601   {
00602     G4ErrorMatrixConstIter m1c1 = mat1.m.begin();
00603     for(G4int c=1;c<=r;c++)
00604     {
00605       G4double tmp = 0.0;
00606       G4ErrorMatrixIter tempri = tempr1;
00607       G4ErrorMatrixConstIter m1ci = m1c1;
00608       for(G4int i=1;i<=mat1.num_col();i++)
00609       {
00610         tmp+=(*(tempri++))*(*(m1ci++));
00611       }
00612       *(mr++) = tmp;
00613       m1c1 += n;
00614     }
00615     tempr1 += n;
00616   }
00617   return mret;
00618 }

G4ErrorSymMatrix G4ErrorSymMatrix::similarityT ( const G4ErrorMatrix m1  )  const

Definition at line 654 of file G4ErrorSymMatrix.cc.

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

00655 {
00656   G4ErrorSymMatrix mret(mat1.num_col());
00657   G4ErrorMatrix temp = (*this)*mat1;
00658   G4int n = mat1.num_col();
00659   G4ErrorMatrixIter mrc = mret.m.begin();
00660   G4ErrorMatrixIter temp1r = temp.m.begin();
00661   for(G4int r=1;r<=mret.num_row();r++)
00662   {
00663     G4ErrorMatrixConstIter m11c = mat1.m.begin();
00664     for(G4int c=1;c<=r;c++)
00665     {
00666       G4double tmp = 0.0;
00667       G4ErrorMatrixIter tempir = temp1r;
00668       G4ErrorMatrixConstIter m1ic = m11c;
00669       for(G4int i=1;i<=mat1.num_row();i++)
00670       {
00671         tmp+=(*(tempir))*(*(m1ic));
00672         tempir += n;
00673         m1ic += n;
00674       }
00675       *(mrc++) = tmp;
00676       m11c++;
00677     }
00678     temp1r++;
00679   }
00680   return mret;
00681 }

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

Definition at line 143 of file G4ErrorSymMatrix.cc.

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

00144 {
00145   G4ErrorSymMatrix mret(max_row-min_row+1);
00146   if(max_row > num_row())
00147     { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
00148   G4ErrorMatrixIter a = mret.m.begin();
00149   G4ErrorMatrixIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
00150   for(G4int irow=1; irow<=mret.num_row(); irow++)
00151   {
00152     G4ErrorMatrixIter b = b1;
00153     for(G4int icol=1; icol<=irow; icol++)
00154     {
00155       *(a++) = *(b++);
00156     }
00157     b1 += irow+min_row-1;
00158   }
00159   return mret;
00160 }

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

Definition at line 162 of file G4ErrorSymMatrix.cc.

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

00163 {
00164   if(row <1 || row+mat1.num_row()-1 > num_row() )
00165     { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
00166   G4ErrorMatrixConstIter a = mat1.m.begin();
00167   G4ErrorMatrixIter b1 = m.begin() + (row+2)*(row-1)/2;
00168   for(G4int irow=1; irow<=mat1.num_row(); irow++)
00169   {
00170     G4ErrorMatrixIter b = b1;
00171     for(G4int icol=1; icol<=irow; icol++)
00172     {
00173       *(b++) = *(a++);
00174     }
00175     b1 += irow+row-1;
00176   }
00177 }

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

Definition at line 124 of file G4ErrorSymMatrix.cc.

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

00125 {
00126   G4ErrorSymMatrix mret(max_row-min_row+1);
00127   if(max_row > num_row())
00128     { G4ErrorMatrix::error("G4ErrorSymMatrix::sub: Index out of range"); }
00129   G4ErrorMatrixIter a = mret.m.begin();
00130   G4ErrorMatrixConstIter b1 = m.begin() + (min_row+2)*(min_row-1)/2;
00131   for(G4int irow=1; irow<=mret.num_row(); irow++)
00132   {
00133     G4ErrorMatrixConstIter b = b1;
00134     for(G4int icol=1; icol<=irow; icol++)
00135     {
00136       *(a++) = *(b++);
00137     }
00138     b1 += irow+min_row-1;
00139   }
00140   return mret;
00141 }

G4ErrorSymMatrix G4ErrorSymMatrix::T (  )  const [inline]

Definition at line 80 of file G4ErrorSymMatrix.icc.

References G4ErrorSymMatrix().

Referenced by G4ErrorFreeTrajState::PropagateError().

00081 {
00082   return G4ErrorSymMatrix(*this);
00083 }

G4double G4ErrorSymMatrix::trace (  )  const

Definition at line 819 of file G4ErrorSymMatrix.cc.

00820 {
00821    G4double t = 0.0;
00822    for (G4int i=0; i<nrow; i++) 
00823      { t += *(m.begin() + (i+3)*i/2); }
00824    return t;
00825 }


Friends And Related Function Documentation

G4double condition ( const G4ErrorSymMatrix m  )  [friend]

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

void diag_step ( G4ErrorSymMatrix t,
G4int  begin,
G4int  end 
) [friend]

G4ErrorMatrix diagonalize ( G4ErrorSymMatrix s  )  [friend]

friend class G4ErrorMatrix [friend]

Definition at line 192 of file G4ErrorSymMatrix.hh.

friend class G4ErrorSymMatrix_row [friend]

Definition at line 190 of file G4ErrorSymMatrix.hh.

friend class G4ErrorSymMatrix_row_const [friend]

Definition at line 191 of file G4ErrorSymMatrix.hh.

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

G4ErrorMatrix operator * ( const G4ErrorMatrix m1,
const G4ErrorSymMatrix m2 
) [friend]

Definition at line 287 of file G4ErrorSymMatrix.cc.

00288 {
00289   G4ErrorMatrix mret(mat1.num_row(),mat2.num_col());
00290   CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
00291   G4ErrorMatrixConstIter mit1, mit2, sp,snp; //mit2=0
00292   G4double temp;
00293   G4ErrorMatrixIter mir=mret.m.begin();
00294   for(mit1=mat1.m.begin();
00295       mit1<mat1.m.begin()+mat1.num_row()*mat1.num_col();
00296       mit1 = mit2)
00297     {
00298       snp=mat2.m.begin();
00299       for(int step=1;step<=mat2.num_row();++step)
00300         {
00301             mit2=mit1;
00302             sp=snp;
00303             snp+=step;
00304             temp=0;
00305             while(sp<snp)
00306             temp+=*(sp++)*(*(mit2++));
00307             if( step<mat2.num_row() ) { // only if we aren't on the last row
00308                 sp+=step-1;
00309                 for(int stept=step+1;stept<=mat2.num_row();stept++)
00310                 {
00311                    temp+=*sp*(*(mit2++));
00312                    if(stept<mat2.num_row()) sp+=stept;
00313                 }
00314             }   // if(step
00315             *(mir++)=temp;
00316         }       // for(step
00317     }   // for(mit1
00318   return mret;
00319 }

G4ErrorMatrix operator * ( const G4ErrorSymMatrix m1,
const G4ErrorMatrix m2 
) [friend]

Definition at line 321 of file G4ErrorSymMatrix.cc.

00322 {
00323   G4ErrorMatrix mret(mat1.num_row(),mat2.num_col());
00324   CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
00325   G4int step,stept;
00326   G4ErrorMatrixConstIter mit1,mit2,sp,snp;
00327   G4double temp;
00328   G4ErrorMatrixIter mir=mret.m.begin();
00329   for(step=1,snp=mat1.m.begin();step<=mat1.num_row();snp+=step++)
00330   {
00331     for(mit1=mat2.m.begin();mit1<mat2.m.begin()+mat2.num_col();mit1++)
00332     {
00333       mit2=mit1;
00334       sp=snp;
00335       temp=0;
00336       while(sp<snp+step)
00337       {
00338         temp+=*mit2*(*(sp++));
00339         mit2+=mat2.num_col();
00340       }
00341       sp+=step-1;
00342       for(stept=step+1;stept<=mat1.num_row();stept++)
00343       {
00344         temp+=*mit2*(*sp);
00345         mit2+=mat2.num_col();
00346         sp+=stept;
00347       }
00348       *(mir++)=temp;
00349     }
00350   }
00351   return mret;
00352 }

G4ErrorMatrix operator * ( const G4ErrorSymMatrix m1,
const G4ErrorSymMatrix m2 
) [friend]

Definition at line 354 of file G4ErrorSymMatrix.cc.

00355 {
00356   G4ErrorMatrix mret(mat1.num_row(),mat1.num_row());
00357   CHK_DIM_1(mat1.num_col(),mat2.num_row(),*);
00358   G4int step1,stept1,step2,stept2;
00359   G4ErrorMatrixConstIter snp1,sp1,snp2,sp2;
00360   G4double temp;
00361   G4ErrorMatrixIter mr = mret.m.begin();
00362   for(step1=1,snp1=mat1.m.begin();step1<=mat1.num_row();snp1+=step1++)
00363   {
00364     for(step2=1,snp2=mat2.m.begin();step2<=mat2.num_row();)
00365     {
00366       sp1=snp1;
00367       sp2=snp2;
00368       snp2+=step2;
00369       temp=0;
00370       if(step1<step2)
00371       {
00372         while(sp1<snp1+step1)
00373           { temp+=(*(sp1++))*(*(sp2++)); }
00374         sp1+=step1-1;
00375         for(stept1=step1+1;stept1!=step2+1;sp1+=stept1++)
00376           { temp+=(*sp1)*(*(sp2++)); }
00377         sp2+=step2-1;
00378         for(stept2=++step2;stept2<=mat2.num_row();sp1+=stept1++,sp2+=stept2++)
00379           { temp+=(*sp1)*(*sp2); }
00380       }
00381       else
00382       {
00383         while(sp2<snp2)
00384           { temp+=(*(sp1++))*(*(sp2++)); }
00385         sp2+=step2-1;
00386         for(stept2=++step2;stept2!=step1+1;sp2+=stept2++)
00387           { temp+=(*(sp1++))*(*sp2); }
00388         sp1+=step1-1;
00389         for(stept1=step1+1;stept1<=mat1.num_row();sp1+=stept1++,sp2+=stept2++)
00390           { temp+=(*sp1)*(*sp2); }
00391       }
00392       *(mr++)=temp;
00393     }
00394   }
00395   return mret;
00396 }

G4ErrorSymMatrix operator+ ( const G4ErrorSymMatrix m1,
const G4ErrorSymMatrix m2 
) [friend]

Definition at line 223 of file G4ErrorSymMatrix.cc.

00225 {
00226   G4ErrorSymMatrix mret(mat1.nrow);
00227   CHK_DIM_1(mat1.nrow, mat2.nrow,+);
00228   SIMPLE_TOP(+)
00229   return mret;
00230 }

G4ErrorSymMatrix operator- ( const G4ErrorSymMatrix m1,
const G4ErrorSymMatrix m2 
) [friend]

Definition at line 252 of file G4ErrorSymMatrix.cc.

00254 {
00255   G4ErrorSymMatrix mret(mat1.num_row());
00256   CHK_DIM_1(mat1.num_row(),mat2.num_row(),-);
00257   SIMPLE_TOP(-)
00258   return mret;
00259 }

void tridiagonal ( G4ErrorSymMatrix a,
G4ErrorMatrix hsm 
) [friend]


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:55 2013 for Geant4 by  doxygen 1.4.7