00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "globals.hh"
00033 #include <iostream>
00034 #include <cmath>
00035
00036 #include "G4ErrorSymMatrix.hh"
00037 #include "G4ErrorMatrix.hh"
00038
00039
00040
00041 #define SIMPLE_UOP(OPER) \
00042 G4ErrorMatrixIter a=m.begin(); \
00043 G4ErrorMatrixIter e=m.begin()+num_size(); \
00044 for(;a<e; a++) (*a) OPER t;
00045
00046 #define SIMPLE_BOP(OPER) \
00047 G4ErrorMatrixIter a=m.begin(); \
00048 G4ErrorMatrixConstIter b=mat2.m.begin(); \
00049 G4ErrorMatrixConstIter e=m.begin()+num_size(); \
00050 for(;a<e; a++, b++) (*a) OPER (*b);
00051
00052 #define SIMPLE_TOP(OPER) \
00053 G4ErrorMatrixConstIter a=mat1.m.begin(); \
00054 G4ErrorMatrixConstIter b=mat2.m.begin(); \
00055 G4ErrorMatrixIter t=mret.m.begin(); \
00056 G4ErrorMatrixConstIter e=mat1.m.begin()+mat1.num_size(); \
00057 for( ;a<e; a++, b++, t++) (*t) = (*a) OPER (*b);
00058
00059 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
00060 if (r1!=r2 || c1!=c2) { \
00061 G4ErrorMatrix::error("Range error in Matrix function " #fun "(1)."); \
00062 }
00063
00064 #define CHK_DIM_1(c1,r2,fun) \
00065 if (c1!=r2) { \
00066 G4ErrorMatrix::error("Range error in Matrix function " #fun "(2)."); \
00067 }
00068
00069
00070
00071 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p)
00072 : m(p*(p+1)/2), nrow(p)
00073 {
00074 size = nrow * (nrow+1) / 2;
00075 m.assign(size,0);
00076 }
00077
00078 G4ErrorSymMatrix::G4ErrorSymMatrix(G4int p, G4int init)
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 }
00103
00104
00105
00106
00107
00108 G4ErrorSymMatrix::~G4ErrorSymMatrix()
00109 {
00110 }
00111
00112 G4ErrorSymMatrix::G4ErrorSymMatrix(const G4ErrorSymMatrix &mat1)
00113 : m(mat1.size), nrow(mat1.nrow), size(mat1.size)
00114 {
00115 m = mat1.m;
00116 }
00117
00118
00119
00120
00121
00122
00123
00124 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int min_row, G4int max_row) const
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 }
00142
00143 G4ErrorSymMatrix G4ErrorSymMatrix::sub(G4int min_row, G4int max_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 }
00161
00162 void G4ErrorSymMatrix::sub(G4int row,const G4ErrorSymMatrix &mat1)
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 }
00178
00179
00180
00181
00182
00183 G4ErrorSymMatrix dsum(const G4ErrorSymMatrix &mat1,
00184 const G4ErrorSymMatrix &mat2)
00185 {
00186 G4ErrorSymMatrix mret(mat1.num_row() + mat2.num_row(), 0);
00187 mret.sub(1,mat1);
00188 mret.sub(mat1.num_row()+1,mat2);
00189 return mret;
00190 }
00191
00192
00193
00194
00195
00196
00197 G4ErrorSymMatrix G4ErrorSymMatrix::operator- () const
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 }
00206
00207 G4ErrorMatrix operator+(const G4ErrorMatrix &mat1, const G4ErrorSymMatrix &mat2)
00208 {
00209 G4ErrorMatrix mret(mat1);
00210 CHK_DIM_2(mat1.num_row(),mat2.num_row(), mat1.num_col(),mat2.num_col(),+);
00211 mret += mat2;
00212 return mret;
00213 }
00214
00215 G4ErrorMatrix operator+(const G4ErrorSymMatrix &mat1, const G4ErrorMatrix &mat2)
00216 {
00217 G4ErrorMatrix mret(mat2);
00218 CHK_DIM_2(mat1.num_row(),mat2.num_row(),mat1.num_col(),mat2.num_col(),+);
00219 mret += mat1;
00220 return mret;
00221 }
00222
00223 G4ErrorSymMatrix operator+(const G4ErrorSymMatrix &mat1,
00224 const G4ErrorSymMatrix &mat2)
00225 {
00226 G4ErrorSymMatrix mret(mat1.nrow);
00227 CHK_DIM_1(mat1.nrow, mat2.nrow,+);
00228 SIMPLE_TOP(+)
00229 return mret;
00230 }
00231
00232
00233
00234
00235
00236 G4ErrorMatrix operator-(const G4ErrorMatrix &mat1, const G4ErrorSymMatrix &mat2)
00237 {
00238 G4ErrorMatrix mret(mat1);
00239 CHK_DIM_2(mat1.num_row(),mat2.num_row(),mat1.num_col(),mat2.num_col(),-);
00240 mret -= mat2;
00241 return mret;
00242 }
00243
00244 G4ErrorMatrix operator-(const G4ErrorSymMatrix &mat1, const G4ErrorMatrix &mat2)
00245 {
00246 G4ErrorMatrix mret(mat1);
00247 CHK_DIM_2(mat1.num_row(),mat2.num_row(),mat1.num_col(),mat2.num_col(),-);
00248 mret -= mat2;
00249 return mret;
00250 }
00251
00252 G4ErrorSymMatrix operator-(const G4ErrorSymMatrix &mat1,
00253 const G4ErrorSymMatrix &mat2)
00254 {
00255 G4ErrorSymMatrix mret(mat1.num_row());
00256 CHK_DIM_1(mat1.num_row(),mat2.num_row(),-);
00257 SIMPLE_TOP(-)
00258 return mret;
00259 }
00260
00261
00262
00263
00264
00265
00266 G4ErrorSymMatrix operator/(const G4ErrorSymMatrix &mat1,G4double t)
00267 {
00268 G4ErrorSymMatrix mret(mat1);
00269 mret /= t;
00270 return mret;
00271 }
00272
00273 G4ErrorSymMatrix operator*(const G4ErrorSymMatrix &mat1,G4double t)
00274 {
00275 G4ErrorSymMatrix mret(mat1);
00276 mret *= t;
00277 return mret;
00278 }
00279
00280 G4ErrorSymMatrix operator*(G4double t,const G4ErrorSymMatrix &mat1)
00281 {
00282 G4ErrorSymMatrix mret(mat1);
00283 mret *= t;
00284 return mret;
00285 }
00286
00287 G4ErrorMatrix operator*(const G4ErrorMatrix &mat1, const G4ErrorSymMatrix &mat2)
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;
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() ) {
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 }
00315 *(mir++)=temp;
00316 }
00317 }
00318 return mret;
00319 }
00320
00321 G4ErrorMatrix operator*(const G4ErrorSymMatrix &mat1, const G4ErrorMatrix &mat2)
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 }
00353
00354 G4ErrorMatrix operator*(const G4ErrorSymMatrix &mat1, const G4ErrorSymMatrix &mat2)
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 }
00397
00398
00399
00400
00401
00402 G4ErrorMatrix & G4ErrorMatrix::operator+=(const G4ErrorSymMatrix &mat2)
00403 {
00404 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=);
00405 G4int n = num_col();
00406 G4ErrorMatrixConstIter sjk = mat2.m.begin();
00407 G4ErrorMatrixIter m1j = m.begin();
00408 G4ErrorMatrixIter mj = m.begin();
00409
00410 for(G4int j=1;j<=num_row();j++)
00411 {
00412 G4ErrorMatrixIter mjk = mj;
00413 G4ErrorMatrixIter mkj = m1j;
00414 for(G4int k=1;k<=j;k++)
00415 {
00416 *(mjk++) += *sjk;
00417 if(j!=k) *mkj += *sjk;
00418 sjk++;
00419 mkj += n;
00420 }
00421 mj += n;
00422 m1j++;
00423 }
00424 return (*this);
00425 }
00426
00427 G4ErrorSymMatrix & G4ErrorSymMatrix::operator+=(const G4ErrorSymMatrix &mat2)
00428 {
00429 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),+=);
00430 SIMPLE_BOP(+=)
00431 return (*this);
00432 }
00433
00434 G4ErrorMatrix & G4ErrorMatrix::operator-=(const G4ErrorSymMatrix &mat2)
00435 {
00436 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
00437 G4int n = num_col();
00438 G4ErrorMatrixConstIter sjk = mat2.m.begin();
00439 G4ErrorMatrixIter m1j = m.begin();
00440 G4ErrorMatrixIter mj = m.begin();
00441
00442 for(G4int j=1;j<=num_row();j++)
00443 {
00444 G4ErrorMatrixIter mjk = mj;
00445 G4ErrorMatrixIter mkj = m1j;
00446 for(G4int k=1;k<=j;k++)
00447 {
00448 *(mjk++) -= *sjk;
00449 if(j!=k) *mkj -= *sjk;
00450 sjk++;
00451 mkj += n;
00452 }
00453 mj += n;
00454 m1j++;
00455 }
00456 return (*this);
00457 }
00458
00459 G4ErrorSymMatrix & G4ErrorSymMatrix::operator-=(const G4ErrorSymMatrix &mat2)
00460 {
00461 CHK_DIM_2(num_row(),mat2.num_row(),num_col(),mat2.num_col(),-=);
00462 SIMPLE_BOP(-=)
00463 return (*this);
00464 }
00465
00466 G4ErrorSymMatrix & G4ErrorSymMatrix::operator/=(G4double t)
00467 {
00468 SIMPLE_UOP(/=)
00469 return (*this);
00470 }
00471
00472 G4ErrorSymMatrix & G4ErrorSymMatrix::operator*=(G4double t)
00473 {
00474 SIMPLE_UOP(*=)
00475 return (*this);
00476 }
00477
00478 G4ErrorMatrix & G4ErrorMatrix::operator=(const G4ErrorSymMatrix &mat1)
00479 {
00480 if(mat1.nrow*mat1.nrow != size)
00481 {
00482 size = mat1.nrow * mat1.nrow;
00483 m.resize(size);
00484 }
00485 nrow = mat1.nrow;
00486 ncol = mat1.nrow;
00487 G4int n = ncol;
00488 G4ErrorMatrixConstIter sjk = mat1.m.begin();
00489 G4ErrorMatrixIter m1j = m.begin();
00490 G4ErrorMatrixIter mj = m.begin();
00491
00492 for(G4int j=1;j<=num_row();j++)
00493 {
00494 G4ErrorMatrixIter mjk = mj;
00495 G4ErrorMatrixIter mkj = m1j;
00496 for(G4int k=1;k<=j;k++)
00497 {
00498 *(mjk++) = *sjk;
00499 if(j!=k) *mkj = *sjk;
00500 sjk++;
00501 mkj += n;
00502 }
00503 mj += n;
00504 m1j++;
00505 }
00506 return (*this);
00507 }
00508
00509 G4ErrorSymMatrix & G4ErrorSymMatrix::operator=(const G4ErrorSymMatrix &mat1)
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 }
00521
00522
00523
00524 std::ostream& operator<<(std::ostream &os, const G4ErrorSymMatrix &q)
00525 {
00526 os << G4endl;
00527
00528
00529
00530
00531 G4int width;
00532 if(os.flags() & std::ios::fixed)
00533 {
00534 width = os.precision()+3;
00535 }
00536 else
00537 {
00538 width = os.precision()+7;
00539 }
00540 for(G4int irow = 1; irow<= q.num_row(); irow++)
00541 {
00542 for(G4int icol = 1; icol <= q.num_col(); icol++)
00543 {
00544 os.width(width);
00545 os << q(irow,icol) << " ";
00546 }
00547 os << G4endl;
00548 }
00549 return os;
00550 }
00551
00552 G4ErrorSymMatrix G4ErrorSymMatrix::
00553 apply(G4double (*f)(G4double, G4int, G4int)) const
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 }
00567
00568 void G4ErrorSymMatrix::assign (const G4ErrorMatrix &mat1)
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 }
00588
00589 G4ErrorSymMatrix G4ErrorSymMatrix::similarity(const G4ErrorMatrix &mat1) const
00590 {
00591 G4ErrorSymMatrix mret(mat1.num_row());
00592 G4ErrorMatrix temp = mat1*(*this);
00593
00594
00595
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 }
00619
00620 G4ErrorSymMatrix G4ErrorSymMatrix::similarity(const G4ErrorSymMatrix &mat1) const
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 }
00653
00654 G4ErrorSymMatrix G4ErrorSymMatrix::similarityT(const G4ErrorMatrix &mat1) const
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 }
00682
00683 void G4ErrorSymMatrix::invert(G4int &ifail)
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;
00797 }
00798
00799 G4double G4ErrorSymMatrix::determinant() const
00800 {
00801 static const G4int max_array = 20;
00802
00803
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 }
00818
00819 G4double G4ErrorSymMatrix::trace() const
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 }
00826
00827 void G4ErrorSymMatrix::invertBunchKaufman(G4int &ifail)
00828 {
00829
00830
00831
00832
00833
00834
00835
00836
00837 G4int i, j, k, ss;
00838 G4int pivrow;
00839
00840
00841
00842
00843
00844
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
00853
00854
00855
00856 G4ErrorMatrixIter x = xvec.begin();
00857
00858 pivIter piv = pivv.begin();
00859
00860
00861 G4double temp1, temp2;
00862 G4ErrorMatrixIter ip, mjj, iq;
00863 G4double lambda, sigma;
00864 const G4double alpha = .6404;
00865 const G4double epsilon = 32*DBL_EPSILON;
00866
00867
00868
00869
00870
00871 for (i = 0; i < nrow; i++)
00872 {
00873 piv[i] = i+1;
00874 }
00875
00876 ifail = 0;
00877
00878
00879
00880
00881
00882 for (j=1; j < nrow; j+=ss)
00883 {
00884 mjj = m.begin() + j*(j-1)/2 + j-1;
00885 lambda = 0;
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;
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)
00935 {
00936 piv[j-1] = pivrow;
00937 if (*mjj == 0)
00938 {
00939 ifail=1;
00940 return;
00941 }
00942 temp2 = *mjj = 1./ *mjj;
00943
00944
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
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)
00965 {
00966 piv[j-1] = pivrow;
00967
00968
00969
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;
00995
00996
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
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
01017 {
01018 piv[j-1] = -pivrow;
01019 piv[j] = 0;
01020
01021 if (j+1 != pivrow)
01022 {
01023
01024
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
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
01059
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)
01067 {
01068
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
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 }
01103
01104 if (j == nrow)
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 }
01117
01118
01119
01120 for (j = nrow ; j >= 1 ; j -= ss)
01121 {
01122 mjj = m.begin() + j*(j-1)/2 + j-1;
01123 if (piv[j-1] > 0)
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
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
01206
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;
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 }
01235
01236 return;
01237 }
01238
01239 G4double G4ErrorSymMatrix::posDefFraction5x5 = 1.0;
01240 G4double G4ErrorSymMatrix::posDefFraction6x6 = 1.0;
01241 G4double G4ErrorSymMatrix::adjustment5x5 = 0.0;
01242 G4double G4ErrorSymMatrix::adjustment6x6 = 0.0;
01243 const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_5x5 = .5;
01244 const G4double G4ErrorSymMatrix::CHOLESKY_THRESHOLD_6x6 = .2;
01245 const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_5x5 = .005;
01246 const G4double G4ErrorSymMatrix::CHOLESKY_CREEP_6x6 = .002;
01247
01248
01249
01250
01251
01252 #define A00 0
01253 #define A01 1
01254 #define A02 3
01255 #define A03 6
01256 #define A04 10
01257 #define A05 15
01258
01259 #define A10 1
01260 #define A11 2
01261 #define A12 4
01262 #define A13 7
01263 #define A14 11
01264 #define A15 16
01265
01266 #define A20 3
01267 #define A21 4
01268 #define A22 5
01269 #define A23 8
01270 #define A24 12
01271 #define A25 17
01272
01273 #define A30 6
01274 #define A31 7
01275 #define A32 8
01276 #define A33 9
01277 #define A34 13
01278 #define A35 18
01279
01280 #define A40 10
01281 #define A41 11
01282 #define A42 12
01283 #define A43 13
01284 #define A44 14
01285 #define A45 19
01286
01287 #define A50 15
01288 #define A51 16
01289 #define A52 17
01290 #define A53 18
01291 #define A54 19
01292 #define A55 20
01293
01294 void G4ErrorSymMatrix::invert5(G4int & ifail)
01295 {
01296 if (posDefFraction5x5 >= CHOLESKY_THRESHOLD_5x5)
01297 {
01298 invertCholesky5(ifail);
01299 posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
01300 if (ifail!=0)
01301 {
01302 invertHaywood5(ifail);
01303 }
01304 }
01305 else
01306 {
01307 if (posDefFraction5x5 + adjustment5x5 >= CHOLESKY_THRESHOLD_5x5)
01308 {
01309 invertCholesky5(ifail);
01310 posDefFraction5x5 = .9*posDefFraction5x5 + .1*(1-ifail);
01311 if (ifail!=0)
01312 {
01313 invertHaywood5(ifail);
01314 adjustment5x5 = 0;
01315 }
01316 }
01317 else
01318 {
01319 invertHaywood5(ifail);
01320 adjustment5x5 += CHOLESKY_CREEP_5x5;
01321 }
01322 }
01323 return;
01324 }
01325
01326 void G4ErrorSymMatrix::invert6(G4int & ifail)
01327 {
01328 if (posDefFraction6x6 >= CHOLESKY_THRESHOLD_6x6)
01329 {
01330 invertCholesky6(ifail);
01331 posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
01332 if (ifail!=0)
01333 {
01334 invertHaywood6(ifail);
01335 }
01336 }
01337 else
01338 {
01339 if (posDefFraction6x6 + adjustment6x6 >= CHOLESKY_THRESHOLD_6x6)
01340 {
01341 invertCholesky6(ifail);
01342 posDefFraction6x6 = .9*posDefFraction6x6 + .1*(1-ifail);
01343 if (ifail!=0)
01344 {
01345 invertHaywood6(ifail);
01346 adjustment6x6 = 0;
01347 }
01348 }
01349 else
01350 {
01351 invertHaywood6(ifail);
01352 adjustment6x6 += CHOLESKY_CREEP_6x6;
01353 }
01354 }
01355 return;
01356 }
01357
01358 void G4ErrorSymMatrix::invertHaywood5 (G4int & ifail)
01359 {
01360 ifail = 0;
01361
01362
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
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
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
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 }
01525
01526 void G4ErrorSymMatrix::invertHaywood6 (G4int & ifail)
01527 {
01528 ifail = 0;
01529
01530
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
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
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
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
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 }
01910
01911 void G4ErrorSymMatrix::invertCholesky5 (G4int & ifail)
01912 {
01913
01914
01915
01916
01917
01918
01919
01920
01921
01922
01923
01924 G4double h10;
01925 G4double h20, h21;
01926 G4double h30, h31, h32;
01927 G4double h40, h41, h42, h43;
01928
01929 G4double h00, h11, h22, h33, h44;
01930
01931
01932 G4double g10;
01933 G4double g20, g21;
01934 G4double g30, g31, g32;
01935 G4double g40, g41, g42, g43;
01936
01937 ifail = 1;
01938
01939
01940
01941
01942
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
01954
01955 h11 = m[A11] - (g10 * g10);
01956 if (h11 <= 0) { return; }
01957 h11 = 1.0 / std::sqrt(h11);
01958
01959
01960
01961
01962 g21 = (m[A21] - (g10 * g20)) * h11;
01963 g31 = (m[A31] - (g10 * g30)) * h11;
01964 g41 = (m[A41] - (g10 * g40)) * h11;
01965
01966
01967
01968 h22 = m[A22] - (g20 * g20) - (g21 * g21);
01969 if (h22 <= 0) { return; }
01970 h22 = 1.0 / std::sqrt(h22);
01971
01972
01973
01974
01975 g32 = (m[A32] - (g20 * g30) - (g21 * g31)) * h22;
01976 g42 = (m[A42] - (g20 * g40) - (g21 * g41)) * h22;
01977
01978
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
01985
01986 g43 = (m[A43] - (g30 * g40) - (g31 * g41) - (g32 * g42)) * h33;
01987
01988
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
01995
01996
01997
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
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 }
02032
02033 void G4ErrorSymMatrix::invertCholesky6 (G4int & ifail)
02034 {
02035
02036
02037
02038
02039
02040
02041
02042
02043
02044
02045 G4double h10;
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;
02052
02053
02054 G4double g10;
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;
02061
02062
02063
02064
02065
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
02078
02079 h11 = m[A11] - (g10 * g10);
02080 if (h11 <= 0) { return; }
02081 h11 = 1.0 / std::sqrt(h11);
02082
02083
02084
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
02092
02093 h22 = m[A22] - (g20 * g20) - (g21 * g21);
02094 if (h22 <= 0) { return; }
02095 h22 = 1.0 / std::sqrt(h22);
02096
02097
02098
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
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
02111
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
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
02123
02124 g54 = (m[A54] - (g40 * g50) - (g41 * g51) - (g42 * g52) - (g43 * g53)) * h44;
02125
02126
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
02133
02134
02135
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
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 }
02181
02182 void G4ErrorSymMatrix::invert4 (G4int & ifail)
02183 {
02184 ifail = 0;
02185
02186
02187
02188 G4double Det2_12_01 = m[A10]*m[A21] - m[A11]*m[A20];
02189 G4double Det2_12_02 = m[A10]*m[A22] - m[A12]*m[A20];
02190 G4double Det2_12_12 = m[A11]*m[A22] - m[A12]*m[A21];
02191 G4double Det2_13_01 = m[A10]*m[A31] - m[A11]*m[A30];
02192 G4double Det2_13_02 = m[A10]*m[A32] - m[A12]*m[A30];
02193 G4double Det2_13_03 = m[A10]*m[A33] - m[A13]*m[A30];
02194 G4double Det2_13_12 = m[A11]*m[A32] - m[A12]*m[A31];
02195 G4double Det2_13_13 = m[A11]*m[A33] - m[A13]*m[A31];
02196 G4double Det2_23_01 = m[A20]*m[A31] - m[A21]*m[A30];
02197 G4double Det2_23_02 = m[A20]*m[A32] - m[A22]*m[A30];
02198 G4double Det2_23_03 = m[A20]*m[A33] - m[A23]*m[A30];
02199 G4double Det2_23_12 = m[A21]*m[A32] - m[A22]*m[A31];
02200 G4double Det2_23_13 = m[A21]*m[A33] - m[A23]*m[A31];
02201 G4double Det2_23_23 = m[A22]*m[A33] - m[A23]*m[A32];
02202
02203
02204
02205 G4double Det3_012_012 = m[A00]*Det2_12_12 - m[A01]*Det2_12_02
02206 + m[A02]*Det2_12_01;
02207 G4double Det3_013_012 = m[A00]*Det2_13_12 - m[A01]*Det2_13_02
02208 + m[A02]*Det2_13_01;
02209 G4double Det3_013_013 = m[A00]*Det2_13_13 - m[A01]*Det2_13_03
02210 + m[A03]*Det2_13_01;
02211 G4double Det3_023_012 = m[A00]*Det2_23_12 - m[A01]*Det2_23_02
02212 + m[A02]*Det2_23_01;
02213 G4double Det3_023_013 = m[A00]*Det2_23_13 - m[A01]*Det2_23_03
02214 + m[A03]*Det2_23_01;
02215 G4double Det3_023_023 = m[A00]*Det2_23_23 - m[A02]*Det2_23_03
02216 + m[A03]*Det2_23_02;
02217 G4double Det3_123_012 = m[A10]*Det2_23_12 - m[A11]*Det2_23_02
02218 + m[A12]*Det2_23_01;
02219 G4double Det3_123_013 = m[A10]*Det2_23_13 - m[A11]*Det2_23_03
02220 + m[A13]*Det2_23_01;
02221 G4double Det3_123_023 = m[A10]*Det2_23_23 - m[A12]*Det2_23_03
02222 + m[A13]*Det2_23_02;
02223 G4double Det3_123_123 = m[A11]*Det2_23_23 - m[A12]*Det2_23_13
02224 + m[A13]*Det2_23_12;
02225
02226
02227
02228 G4double det = m[A00]*Det3_123_123
02229 - m[A01]*Det3_123_023
02230 + m[A02]*Det3_123_013
02231 - m[A03]*Det3_123_012;
02232
02233 if ( det == 0 )
02234 {
02235 ifail = 1;
02236 return;
02237 }
02238
02239 G4double oneOverDet = 1.0/det;
02240 G4double mn1OverDet = - oneOverDet;
02241
02242 m[A00] = Det3_123_123 * oneOverDet;
02243 m[A01] = Det3_123_023 * mn1OverDet;
02244 m[A02] = Det3_123_013 * oneOverDet;
02245 m[A03] = Det3_123_012 * mn1OverDet;
02246
02247
02248 m[A11] = Det3_023_023 * oneOverDet;
02249 m[A12] = Det3_023_013 * mn1OverDet;
02250 m[A13] = Det3_023_012 * oneOverDet;
02251
02252 m[A22] = Det3_013_013 * oneOverDet;
02253 m[A23] = Det3_013_012 * mn1OverDet;
02254
02255 m[A33] = Det3_012_012 * oneOverDet;
02256
02257 return;
02258 }
02259
02260 void G4ErrorSymMatrix::invertHaywood4 (G4int & ifail)
02261 {
02262 invert4(ifail);
02263
02264 }
02265