CLHEP VERSION Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // -*- C++ -*- 00002 // --------------------------------------------------------------------------- 00003 // 00004 // This file is a part of the CLHEP - a Class Library for High Energy Physics. 00005 // 00006 00007 #ifdef GNUPRAGMA 00008 #pragma implementation 00009 #endif 00010 00011 #include <string.h> 00012 #include <float.h> // for DBL_EPSILON 00013 #include <cmath> 00014 #include <stdlib.h> 00015 00016 #include "CLHEP/Matrix/defs.h" 00017 #include "CLHEP/Random/Random.h" 00018 #include "CLHEP/Matrix/Matrix.h" 00019 #include "CLHEP/Matrix/SymMatrix.h" 00020 #include "CLHEP/Matrix/DiagMatrix.h" 00021 #include "CLHEP/Matrix/Vector.h" 00022 00023 #ifdef HEP_DEBUG_INLINE 00024 #include "CLHEP/Matrix/Matrix.icc" 00025 #endif 00026 00027 namespace CLHEP { 00028 00029 // Simple operation for all elements 00030 00031 #define SIMPLE_UOP(OPER) \ 00032 mIter a=m.begin(); \ 00033 mIter e=m.end(); \ 00034 for(;a!=e; a++) (*a) OPER t; 00035 00036 #define SIMPLE_BOP(OPER) \ 00037 HepMatrix::mIter a=m.begin(); \ 00038 HepMatrix::mcIter b=m2.m.begin(); \ 00039 HepMatrix::mIter e=m.end(); \ 00040 for(;a!=e; a++, b++) (*a) OPER (*b); 00041 00042 #define SIMPLE_TOP(OPER) \ 00043 HepMatrix::mcIter a=m1.m.begin(); \ 00044 HepMatrix::mcIter b=m2.m.begin(); \ 00045 HepMatrix::mIter t=mret.m.begin(); \ 00046 HepMatrix::mcIter e=m1.m.end(); \ 00047 for(;a!=e; a++, b++, t++) (*t) = (*a) OPER (*b); 00048 00049 // Static functions. 00050 00051 #define CHK_DIM_2(r1,r2,c1,c2,fun) \ 00052 if (r1!=r2 || c1!=c2) { \ 00053 HepGenMatrix::error("Range error in Matrix function " #fun "(1)."); \ 00054 } 00055 00056 #define CHK_DIM_1(c1,r2,fun) \ 00057 if (c1!=r2) { \ 00058 HepGenMatrix::error("Range error in Matrix function " #fun "(2)."); \ 00059 } 00060 00061 // Constructors. (Default constructors are inlined and in .icc file) 00062 00063 HepMatrix::HepMatrix(int p,int q) 00064 : m(p*q), nrow(p), ncol(q) 00065 { 00066 size_ = nrow * ncol; 00067 } 00068 00069 HepMatrix::HepMatrix(int p,int q,int init) 00070 : m(p*q), nrow(p), ncol(q) 00071 { 00072 size_ = nrow * ncol; 00073 00074 if (size_ > 0) { 00075 switch(init) 00076 { 00077 case 0: 00078 break; 00079 00080 case 1: 00081 { 00082 if ( ncol == nrow ) { 00083 mIter a = m.begin(); 00084 for( int step=0; step < size_; step+=(ncol+1) ) *(a+step) = 1.0; 00085 } else { 00086 error("Invalid dimension in HepMatrix(int,int,1)."); 00087 } 00088 break; 00089 } 00090 default: 00091 error("Matrix: initialization must be either 0 or 1."); 00092 } 00093 } 00094 } 00095 00096 HepMatrix::HepMatrix(int p,int q, HepRandom &r) 00097 : m(p*q), nrow(p), ncol(q) 00098 { 00099 size_ = nrow * ncol; 00100 00101 mIter a = m.begin(); 00102 mIter b = m.end(); 00103 for(; a<b; a++) *a = r(); 00104 } 00105 // 00106 // Destructor 00107 // 00108 HepMatrix::~HepMatrix() { 00109 } 00110 00111 HepMatrix::HepMatrix(const HepMatrix &m1) 00112 : HepGenMatrix(m1), m(m1.size_), nrow(m1.nrow), ncol(m1.ncol), size_(m1.size_) 00113 { 00114 m = m1.m; 00115 00116 } 00117 00118 // trivial 00119 00120 int HepMatrix::num_row() const { return nrow;} 00121 00122 int HepMatrix::num_col() const { return ncol;} 00123 00124 int HepMatrix::num_size() const { return size_;} 00125 00126 // operator() 00127 00128 double & HepMatrix::operator()(int row, int col) 00129 { 00130 #ifdef MATRIX_BOUND_CHECK 00131 if(row<1 || row>num_row() || col<1 || col>num_col()) 00132 error("Range error in HepMatrix::operator()"); 00133 #endif 00134 return *(m.begin()+(row-1)*ncol+col-1); 00135 } 00136 00137 const double & HepMatrix::operator()(int row, int col) const 00138 { 00139 #ifdef MATRIX_BOUND_CHECK 00140 if(row<1 || row>num_row() || col<1 || col>num_col()) 00141 error("Range error in HepMatrix::operator()"); 00142 #endif 00143 return *(m.begin()+(row-1)*ncol+col-1); 00144 } 00145 00146 00147 HepMatrix::HepMatrix(const HepSymMatrix &m1) 00148 : m(m1.nrow*m1.nrow), nrow(m1.nrow), ncol(m1.nrow) 00149 { 00150 size_ = nrow * ncol; 00151 00152 mcIter sjk = m1.m.begin(); 00153 // j >= k 00154 for(int j=0; j!=nrow; ++j) { 00155 for(int k=0; k<=j; ++k) { 00156 m[j*ncol+k] = *sjk; 00157 // we could copy the diagonal element twice or check 00158 // doing the check may be a tiny bit faster, 00159 // so we choose that option for now 00160 if(k!=j) m[k*nrow+j] = *sjk; 00161 ++sjk; 00162 } 00163 } 00164 } 00165 00166 HepMatrix::HepMatrix(const HepDiagMatrix &m1) 00167 : m(m1.nrow*m1.nrow), nrow(m1.nrow), ncol(m1.nrow) 00168 { 00169 size_ = nrow * ncol; 00170 00171 int n = num_row(); 00172 mIter mrr; 00173 mcIter mr = m1.m.begin(); 00174 for(int r=0;r<n;r++) { 00175 mrr = m.begin()+(n+1)*r; 00176 *mrr = *(mr++); 00177 } 00178 } 00179 00180 HepMatrix::HepMatrix(const HepVector &m1) 00181 : m(m1.nrow), nrow(m1.nrow), ncol(1) 00182 { 00183 00184 size_ = nrow; 00185 m = m1.m; 00186 } 00187 00188 00189 // 00190 // 00191 // Sub matrix 00192 // 00193 // 00194 00195 HepMatrix HepMatrix::sub(int min_row, int max_row, 00196 int min_col,int max_col) const 00197 #ifdef HEP_GNU_OPTIMIZED_RETURN 00198 return mret(max_row-min_row+1,max_col-min_col+1); 00199 { 00200 #else 00201 { 00202 HepMatrix mret(max_row-min_row+1,max_col-min_col+1); 00203 #endif 00204 if(max_row > num_row() || max_col >num_col()) 00205 error("HepMatrix::sub: Index out of range"); 00206 mIter a = mret.m.begin(); 00207 int nc = num_col(); 00208 mcIter b1 = m.begin() + (min_row - 1) * nc + min_col - 1; 00209 int rowsize = mret.num_row(); 00210 for(int irow=1; irow<=rowsize; ++irow) { 00211 mcIter brc = b1; 00212 for(int icol=0; icol<mret.num_col(); ++icol) { 00213 *(a++) = *(brc++); 00214 } 00215 if(irow<rowsize) b1 += nc; 00216 } 00217 return mret; 00218 } 00219 00220 void HepMatrix::sub(int row,int col,const HepMatrix &m1) 00221 { 00222 if(row <1 || row+m1.num_row()-1 > num_row() || 00223 col <1 || col+m1.num_col()-1 > num_col() ) 00224 error("HepMatrix::sub: Index out of range"); 00225 mcIter a = m1.m.begin(); 00226 int nc = num_col(); 00227 mIter b1 = m.begin() + (row - 1) * nc + col - 1; 00228 int rowsize = m1.num_row(); 00229 for(int irow=1; irow<=rowsize; ++irow) { 00230 mIter brc = b1; 00231 for(int icol=0; icol<m1.num_col(); ++icol) { 00232 *(brc++) = *(a++); 00233 } 00234 if(irow<rowsize) b1 += nc; 00235 } 00236 } 00237 00238 // 00239 // Direct sum of two matricies 00240 // 00241 00242 HepMatrix dsum(const HepMatrix &m1, const HepMatrix &m2) 00243 #ifdef HEP_GNU_OPTIMIZED_RETURN 00244 return mret(m1.num_row() + m2.num_row(), m1.num_col() + m2.num_col(), 00245 0); 00246 { 00247 #else 00248 { 00249 HepMatrix mret(m1.num_row() + m2.num_row(), m1.num_col() + m2.num_col(), 00250 0); 00251 #endif 00252 mret.sub(1,1,m1); 00253 mret.sub(m1.num_row()+1,m1.num_col()+1,m2); 00254 return mret; 00255 } 00256 00257 /* ----------------------------------------------------------------------- 00258 This section contains support routines for matrix.h. This section contains 00259 The two argument functions +,-. They call the copy constructor and +=,-=. 00260 ----------------------------------------------------------------------- */ 00261 HepMatrix HepMatrix::operator- () const 00262 #ifdef HEP_GNU_OPTIMIZED_RETURN 00263 return m2(nrow, ncol); 00264 { 00265 #else 00266 { 00267 HepMatrix m2(nrow, ncol); 00268 #endif 00269 mcIter a=m.begin(); 00270 mIter b=m2.m.begin(); 00271 mcIter e=m.end(); 00272 for(;a<e; a++, b++) (*b) = -(*a); 00273 return m2; 00274 } 00275 00276 00277 00278 HepMatrix operator+(const HepMatrix &m1,const HepMatrix &m2) 00279 #ifdef HEP_GNU_OPTIMIZED_RETURN 00280 return mret(m1.nrow, m1.ncol); 00281 { 00282 #else 00283 { 00284 HepMatrix mret(m1.nrow, m1.ncol); 00285 #endif 00286 CHK_DIM_2(m1.num_row(),m2.num_row(), m1.num_col(),m2.num_col(),+); 00287 SIMPLE_TOP(+) 00288 return mret; 00289 } 00290 00291 // 00292 // operator - 00293 // 00294 00295 HepMatrix operator-(const HepMatrix &m1,const HepMatrix &m2) 00296 #ifdef HEP_GNU_OPTIMIZED_RETURN 00297 return mret(m1.num_row(), m1.num_col()); 00298 { 00299 #else 00300 { 00301 HepMatrix mret(m1.num_row(), m1.num_col()); 00302 #endif 00303 CHK_DIM_2(m1.num_row(),m2.num_row(), 00304 m1.num_col(),m2.num_col(),-); 00305 SIMPLE_TOP(-) 00306 return mret; 00307 } 00308 00309 /* ----------------------------------------------------------------------- 00310 This section contains support routines for matrix.h. This file contains 00311 The two argument functions *,/. They call copy constructor and then /=,*=. 00312 ----------------------------------------------------------------------- */ 00313 00314 HepMatrix operator/( 00315 const HepMatrix &m1,double t) 00316 #ifdef HEP_GNU_OPTIMIZED_RETURN 00317 return mret(m1); 00318 { 00319 #else 00320 { 00321 HepMatrix mret(m1); 00322 #endif 00323 mret /= t; 00324 return mret; 00325 } 00326 00327 HepMatrix operator*(const HepMatrix &m1,double t) 00328 #ifdef HEP_GNU_OPTIMIZED_RETURN 00329 return mret(m1); 00330 { 00331 #else 00332 { 00333 HepMatrix mret(m1); 00334 #endif 00335 mret *= t; 00336 return mret; 00337 } 00338 00339 HepMatrix operator*(double t,const HepMatrix &m1) 00340 #ifdef HEP_GNU_OPTIMIZED_RETURN 00341 return mret(m1); 00342 { 00343 #else 00344 { 00345 HepMatrix mret(m1); 00346 #endif 00347 mret *= t; 00348 return mret; 00349 } 00350 00351 HepMatrix operator*(const HepMatrix &m1,const HepMatrix &m2) 00352 #ifdef HEP_GNU_OPTIMIZED_RETURN 00353 return mret(m1.nrow,m2.ncol,0); 00354 { 00355 #else 00356 { 00357 // initialize matrix to 0.0 00358 HepMatrix mret(m1.nrow,m2.ncol,0); 00359 #endif 00360 CHK_DIM_1(m1.ncol,m2.nrow,*); 00361 00362 int m1cols = m1.ncol; 00363 int m2cols = m2.ncol; 00364 00365 for (int i=0; i<m1.nrow; i++) 00366 { 00367 for (int j=0; j<m1cols; j++) 00368 { 00369 register double temp = m1.m[i*m1cols+j]; 00370 HepMatrix::mIter pt = mret.m.begin() + i*m2cols; 00371 00372 // Loop over k (the column index in matrix m2) 00373 HepMatrix::mcIter pb = m2.m.begin() + m2cols*j; 00374 const HepMatrix::mcIter pblast = pb + m2cols; 00375 while (pb < pblast) 00376 { 00377 (*pt) += temp * (*pb); 00378 pb++; 00379 pt++; 00380 } 00381 } 00382 } 00383 00384 return mret; 00385 } 00386 00387 /* ----------------------------------------------------------------------- 00388 This section contains the assignment and inplace operators =,+=,-=,*=,/=. 00389 ----------------------------------------------------------------------- */ 00390 00391 HepMatrix & HepMatrix::operator+=(const HepMatrix &m2) 00392 { 00393 CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),+=); 00394 SIMPLE_BOP(+=) 00395 return (*this); 00396 } 00397 00398 HepMatrix & HepMatrix::operator-=(const HepMatrix &m2) 00399 { 00400 CHK_DIM_2(num_row(),m2.num_row(),num_col(),m2.num_col(),-=); 00401 SIMPLE_BOP(-=) 00402 return (*this); 00403 } 00404 00405 HepMatrix & HepMatrix::operator/=(double t) 00406 { 00407 SIMPLE_UOP(/=) 00408 return (*this); 00409 } 00410 00411 HepMatrix & HepMatrix::operator*=(double t) 00412 { 00413 SIMPLE_UOP(*=) 00414 return (*this); 00415 } 00416 00417 HepMatrix & HepMatrix::operator=(const HepMatrix &m1) 00418 { 00419 if(m1.nrow*m1.ncol != size_) //??fixme?? m1.size != size 00420 { 00421 size_ = m1.nrow * m1.ncol; 00422 m.resize(size_); //??fixme?? if (size < m1.size) m.resize(m1.size); 00423 } 00424 nrow = m1.nrow; 00425 ncol = m1.ncol; 00426 m = m1.m; 00427 return (*this); 00428 } 00429 00430 // HepMatrix & HepMatrix::operator=(const HepRotation &m2) 00431 // is now in Matrix=Rotation.cc 00432 00433 // Print the Matrix. 00434 00435 std::ostream& operator<<(std::ostream &s, const HepMatrix &q) 00436 { 00437 s << "\n"; 00438 /* Fixed format needs 3 extra characters for field, while scientific needs 7 */ 00439 int width; 00440 if(s.flags() & std::ios::fixed) 00441 width = s.precision()+3; 00442 else 00443 width = s.precision()+7; 00444 for(int irow = 1; irow<= q.num_row(); irow++) 00445 { 00446 for(int icol = 1; icol <= q.num_col(); icol++) 00447 { 00448 s.width(width); 00449 s << q(irow,icol) << " "; 00450 } 00451 s << std::endl; 00452 } 00453 return s; 00454 } 00455 00456 HepMatrix HepMatrix::T() const 00457 #ifdef HEP_GNU_OPTIMIZED_RETURN 00458 return mret(ncol,nrow); 00459 { 00460 #else 00461 { 00462 HepMatrix mret(ncol,nrow); 00463 #endif 00464 mcIter pme = m.begin(); 00465 mIter pt = mret.m.begin(); 00466 for( int nr=0; nr<nrow; ++nr) { 00467 for( int nc=0; nc<ncol; ++nc) { 00468 pt = mret.m.begin() + nr + nrow*nc; 00469 (*pt) = (*pme); 00470 ++pme; 00471 } 00472 } 00473 return mret; 00474 } 00475 00476 HepMatrix HepMatrix::apply(double (*f)(double, int, int)) const 00477 #ifdef HEP_GNU_OPTIMIZED_RETURN 00478 return mret(num_row(),num_col()); 00479 { 00480 #else 00481 { 00482 HepMatrix mret(num_row(),num_col()); 00483 #endif 00484 mcIter a = m.begin(); 00485 mIter b = mret.m.begin(); 00486 for(int ir=1;ir<=num_row();ir++) { 00487 for(int ic=1;ic<=num_col();ic++) { 00488 *(b++) = (*f)(*(a++), ir, ic); 00489 } 00490 } 00491 return mret; 00492 } 00493 00494 int HepMatrix::dfinv_matrix(int *ir) { 00495 if (num_col()!=num_row()) 00496 error("dfinv_matrix: Matrix is not NxN"); 00497 int n = num_col(); 00498 if (n==1) return 0; 00499 00500 double s31, s32; 00501 register double s33, s34; 00502 00503 mIter m11 = m.begin(); 00504 mIter m12 = m11 + 1; 00505 mIter m21 = m11 + n; 00506 mIter m22 = m12 + n; 00507 *m21 = -(*m22) * (*m11) * (*m21); 00508 *m12 = -(*m12); 00509 if (n>2) { 00510 mIter mimim = m11 + n + 1; 00511 for (int i=3;i<=n;i++) { 00512 // calculate these to avoid pointing off the end of the storage array 00513 mIter mi = m11 + (i-1) * n; 00514 mIter mii= m11 + (i-1) * n + i - 1; 00515 int im2 = i - 2; 00516 mIter mj = m11; 00517 mIter mji = mj + i - 1; 00518 mIter mij = mi; 00519 for (int j=1;j<=im2;j++) { 00520 s31 = 0.0; 00521 s32 = *mji; 00522 mIter mkj = mj + j - 1; 00523 mIter mik = mi + j - 1; 00524 mIter mjkp = mj + j; 00525 mIter mkpi = mj + n + i - 1; 00526 for (int k=j;k<=im2;k++) { 00527 s31 += (*mkj) * (*(mik++)); 00528 s32 += (*(mjkp++)) * (*mkpi); 00529 mkj += n; 00530 mkpi += n; 00531 } // for k 00532 *mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31)); 00533 *mji = -s32; 00534 mj += n; 00535 mji += n; 00536 mij++; 00537 } // for j 00538 *(mii-1) = -(*mii) * (*mimim) * (*(mii-1)); 00539 *(mimim+1) = -(*(mimim+1)); 00540 mimim += (n+1); 00541 } // for i 00542 } // n>2 00543 mIter mi = m11; 00544 mIter mii = m11; 00545 for (int i=1;i<n;i++) { 00546 int ni = n - i; 00547 mIter mij = mi; 00548 int j; 00549 for (j=1; j<=i;j++) { 00550 s33 = *mij; 00551 // change initial definition of mikj to avoid pointing off the end of the storage array 00552 mIter mikj = mi + j - 1; 00553 mIter miik = mii + 1; 00554 mIter min_end = mi + n; 00555 for (;miik<min_end;) { 00556 // iterate by n as we enter the loop to avoid pointing off the end of the storage array 00557 mikj += n; 00558 s33 += (*mikj) * (*(miik++)); 00559 } 00560 *(mij++) = s33; 00561 } 00562 for (j=1;j<=ni;j++) { 00563 s34 = 0.0; 00564 mIter miik = mii + j; 00565 for (int k=j;k<=ni;k++) { 00566 // calculate mikij here to avoid pointing off the end of the storage array 00567 mIter mikij = mii + k * n + j; 00568 s34 += *mikij * (*(miik++)); 00569 } 00570 *(mii+j) = s34; 00571 } 00572 mi += n; 00573 mii += (n+1); 00574 } // for i 00575 int nxch = ir[n]; 00576 if (nxch==0) return 0; 00577 for (int mm=1;mm<=nxch;mm++) { 00578 int k = nxch - mm + 1; 00579 int ij = ir[k]; 00580 int i = ij >> 12; 00581 int j = ij%4096; 00582 for (k=1; k<=n;k++) { 00583 // avoid setting the iterator beyond the end of the storage vector 00584 mIter mki = m11 + (k-1)*n + i - 1; 00585 mIter mkj = m11 + (k-1)*n + j - 1; 00586 // 2/24/05 David Sachs fix of improper swap bug that was present 00587 // for many years: 00588 double ti = *mki; // 2/24/05 00589 *mki = *mkj; 00590 *mkj = ti; // 2/24/05 00591 } 00592 } // for mm 00593 return 0; 00594 } 00595 00596 int HepMatrix::dfact_matrix(double &det, int *ir) { 00597 if (ncol!=nrow) 00598 error("dfact_matrix: Matrix is not NxN"); 00599 00600 int ifail, jfail; 00601 int n = ncol; 00602 00603 double tf; 00604 double g1 = 1.0e-19, g2 = 1.0e19; 00605 00606 double p, q, t; 00607 double s11, s12; 00608 00609 double epsilon = 8*DBL_EPSILON; 00610 // could be set to zero (like it was before) 00611 // but then the algorithm often doesn't detect 00612 // that a matrix is singular 00613 00614 int normal = 0, imposs = -1; 00615 int jrange = 0, jover = 1, junder = -1; 00616 ifail = normal; 00617 jfail = jrange; 00618 int nxch = 0; 00619 det = 1.0; 00620 mIter mj = m.begin(); 00621 mIter mjj = mj; 00622 for (int j=1;j<=n;j++) { 00623 int k = j; 00624 p = (fabs(*mjj)); 00625 if (j!=n) { 00626 // replace mij with calculation of position 00627 for (int i=j+1;i<n;i++) { 00628 q = (fabs(*(mj + n*(i-j) + j - 1))); 00629 if (q > p) { 00630 k = i; 00631 p = q; 00632 } 00633 } // for i 00634 if (k==j) { 00635 if (p <= epsilon) { 00636 det = 0; 00637 ifail = imposs; 00638 jfail = jrange; 00639 return ifail; 00640 } 00641 det = -det; // in this case the sign of the determinant 00642 // must not change. So I change it twice. 00643 } // k==j 00644 mIter mjl = mj; 00645 mIter mkl = m.begin() + (k-1)*n; 00646 for (int l=1;l<=n;l++) { 00647 tf = *mjl; 00648 *(mjl++) = *mkl; 00649 *(mkl++) = tf; 00650 } 00651 nxch = nxch + 1; // this makes the determinant change its sign 00652 ir[nxch] = (((j)<<12)+(k)); 00653 } else { // j!=n 00654 if (p <= epsilon) { 00655 det = 0.0; 00656 ifail = imposs; 00657 jfail = jrange; 00658 return ifail; 00659 } 00660 } // j!=n 00661 det *= *mjj; 00662 *mjj = 1.0 / *mjj; 00663 t = (fabs(det)); 00664 if (t < g1) { 00665 det = 0.0; 00666 if (jfail == jrange) jfail = junder; 00667 } else if (t > g2) { 00668 det = 1.0; 00669 if (jfail==jrange) jfail = jover; 00670 } 00671 // calculate mk and mkjp so we don't point off the end of the vector 00672 if (j!=n) { 00673 mIter mjk = mj + j; 00674 for (k=j+1;k<=n;k++) { 00675 mIter mk = mj + n*(k-j); 00676 mIter mkjp = mk + j; 00677 s11 = - (*mjk); 00678 s12 = - (*mkjp); 00679 if (j!=1) { 00680 mIter mik = m.begin() + k - 1; 00681 mIter mijp = m.begin() + j; 00682 mIter mki = mk; 00683 mIter mji = mj; 00684 for (int i=1;i<j;i++) { 00685 s11 += (*mik) * (*(mji++)); 00686 s12 += (*mijp) * (*(mki++)); 00687 mik += n; 00688 mijp += n; 00689 } // for i 00690 } // j!=1 00691 *(mjk++) = -s11 * (*mjj); 00692 *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12)); 00693 } // for k 00694 } // j!=n 00695 // avoid setting the iterator beyond the end of the vector 00696 if(j!=n) { 00697 mj += n; 00698 mjj += (n+1); 00699 } 00700 } // for j 00701 if (nxch%2==1) det = -det; 00702 if (jfail !=jrange) det = 0.0; 00703 ir[n] = nxch; 00704 return 0; 00705 } 00706 00707 void HepMatrix::invert(int &ierr) { 00708 if(ncol != nrow) 00709 error("HepMatrix::invert: Matrix is not NxN"); 00710 00711 static int max_array = 20; 00712 static int *ir = new int [max_array+1]; 00713 00714 if (ncol > max_array) { 00715 delete [] ir; 00716 max_array = nrow; 00717 ir = new int [max_array+1]; 00718 } 00719 double t1, t2, t3; 00720 double det, temp, s; 00721 int ifail; 00722 switch(nrow) { 00723 case 3: 00724 double c11,c12,c13,c21,c22,c23,c31,c32,c33; 00725 ifail = 0; 00726 c11 = (*(m.begin()+4)) * (*(m.begin()+8)) - (*(m.begin()+5)) * (*(m.begin()+7)); 00727 c12 = (*(m.begin()+5)) * (*(m.begin()+6)) - (*(m.begin()+3)) * (*(m.begin()+8)); 00728 c13 = (*(m.begin()+3)) * (*(m.begin()+7)) - (*(m.begin()+4)) * (*(m.begin()+6)); 00729 c21 = (*(m.begin()+7)) * (*(m.begin()+2)) - (*(m.begin()+8)) * (*(m.begin()+1)); 00730 c22 = (*(m.begin()+8)) * (*m.begin()) - (*(m.begin()+6)) * (*(m.begin()+2)); 00731 c23 = (*(m.begin()+6)) * (*(m.begin()+1)) - (*(m.begin()+7)) * (*m.begin()); 00732 c31 = (*(m.begin()+1)) * (*(m.begin()+5)) - (*(m.begin()+2)) * (*(m.begin()+4)); 00733 c32 = (*(m.begin()+2)) * (*(m.begin()+3)) - (*m.begin()) * (*(m.begin()+5)); 00734 c33 = (*m.begin()) * (*(m.begin()+4)) - (*(m.begin()+1)) * (*(m.begin()+3)); 00735 t1 = fabs(*m.begin()); 00736 t2 = fabs(*(m.begin()+3)); 00737 t3 = fabs(*(m.begin()+6)); 00738 if (t1 >= t2) { 00739 if (t3 >= t1) { 00740 temp = *(m.begin()+6); 00741 det = c23*c12-c22*c13; 00742 } else { 00743 temp = *(m.begin()); 00744 det = c22*c33-c23*c32; 00745 } 00746 } else if (t3 >= t2) { 00747 temp = *(m.begin()+6); 00748 det = c23*c12-c22*c13; 00749 } else { 00750 temp = *(m.begin()+3); 00751 det = c13*c32-c12*c33; 00752 } 00753 if (det==0) { 00754 ierr = 1; 00755 return; 00756 } 00757 { 00758 double s1 = temp/det; 00759 mIter mm = m.begin(); 00760 *(mm++) = s1*c11; 00761 *(mm++) = s1*c21; 00762 *(mm++) = s1*c31; 00763 *(mm++) = s1*c12; 00764 *(mm++) = s1*c22; 00765 *(mm++) = s1*c32; 00766 *(mm++) = s1*c13; 00767 *(mm++) = s1*c23; 00768 *(mm) = s1*c33; 00769 } 00770 break; 00771 case 2: 00772 ifail = 0; 00773 det = (*m.begin())*(*(m.begin()+3)) - (*(m.begin()+1))*(*(m.begin()+2)); 00774 if (det==0) { 00775 ierr = 1; 00776 return; 00777 } 00778 s = 1.0/det; 00779 temp = s*(*(m.begin()+3)); 00780 *(m.begin()+1) *= -s; 00781 *(m.begin()+2) *= -s; 00782 *(m.begin()+3) = s*(*m.begin()); 00783 *(m.begin()) = temp; 00784 break; 00785 case 1: 00786 ifail = 0; 00787 if ((*(m.begin()))==0) { 00788 ierr = 1; 00789 return; 00790 } 00791 *(m.begin()) = 1.0/(*(m.begin())); 00792 break; 00793 case 4: 00794 invertHaywood4(ierr); 00795 return; 00796 case 5: 00797 invertHaywood5(ierr); 00798 return; 00799 case 6: 00800 invertHaywood6(ierr); 00801 return; 00802 default: 00803 ifail = dfact_matrix(det, ir); 00804 if(ifail) { 00805 ierr = 1; 00806 return; 00807 } 00808 dfinv_matrix(ir); 00809 break; 00810 } 00811 ierr = 0; 00812 return; 00813 } 00814 00815 double HepMatrix::determinant() const { 00816 static int max_array = 20; 00817 static int *ir = new int [max_array+1]; 00818 if(ncol != nrow) 00819 error("HepMatrix::determinant: Matrix is not NxN"); 00820 if (ncol > max_array) { 00821 delete [] ir; 00822 max_array = nrow; 00823 ir = new int [max_array+1]; 00824 } 00825 double det; 00826 HepMatrix mt(*this); 00827 int i = mt.dfact_matrix(det, ir); 00828 if(i==0) return det; 00829 return 0; 00830 } 00831 00832 double HepMatrix::trace() const { 00833 double t = 0.0; 00834 for (mcIter d = m.begin(); d < m.end(); d += (ncol+1) ) 00835 t += *d; 00836 return t; 00837 } 00838 00839 } // namespace CLHEP