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