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 #ifdef GNUPRAGMA 00038 #pragma implementation 00039 #endif 00040 00041 #include <string.h> 00042 00043 #include "CLHEP/Matrix/defs.h" 00044 #include "CLHEP/Random/Random.h" 00045 #include "CLHEP/Vector/ThreeVector.h" 00046 #include "CLHEP/Matrix/Vector.h" 00047 #include "CLHEP/Matrix/Matrix.h" 00048 00049 #ifdef HEP_DEBUG_INLINE 00050 #include "CLHEP/Matrix/Vector.icc" 00051 #endif 00052 00053 namespace CLHEP { 00054 00055 // Simple operation for all elements 00056 00057 #define SIMPLE_UOP(OPER) \ 00058 register HepGenMatrix::mIter a=m.begin(); \ 00059 register HepGenMatrix::mIter e=m.begin()+num_size(); \ 00060 for(;a<e; a++) (*a) OPER t; 00061 00062 #define SIMPLE_BOP(OPER) \ 00063 register mIter a=m.begin(); \ 00064 register mcIter b=m2.m.begin(); \ 00065 register mcIter e=m.begin()+num_size(); \ 00066 for(;a<e; a++, b++) (*a) OPER (*b); 00067 00068 #define SIMPLE_TOP(OPER) \ 00069 register HepGenMatrix::mcIter a=m1.m.begin(); \ 00070 register HepGenMatrix::mcIter b=m2.m.begin(); \ 00071 register HepGenMatrix::mIter t=mret.m.begin(); \ 00072 register HepGenMatrix::mcIter e=m1.m.begin()+m1.num_size(); \ 00073 for( ;a<e; a++, b++, t++) (*t) = (*a) OPER (*b); 00074 00075 #define CHK_DIM_2(r1,r2,c1,c2,fun) \ 00076 if (r1!=r2 || c1!=c2) { \ 00077 HepGenMatrix::error("Range error in Vector function " #fun "(1)."); \ 00078 } 00079 00080 #define CHK_DIM_1(c1,r2,fun) \ 00081 if (c1!=r2) { \ 00082 HepGenMatrix::error("Range error in Vector function " #fun "(2)."); \ 00083 } 00084 00085 // Constructors. (Default constructors are inlined and in .icc file) 00086 00087 HepVector::HepVector(int p) 00088 : m(p), nrow(p) 00089 { 00090 } 00091 00092 HepVector::HepVector(int p, int init) 00093 : m(p), nrow(p) 00094 { 00095 switch (init) 00096 { 00097 case 0: 00098 m.assign(p,0); 00099 break; 00100 00101 case 1: 00102 { 00103 mIter e = m.begin() + nrow; 00104 for (mIter i=m.begin(); i<e; i++) *i = 1.0; 00105 break; 00106 } 00107 00108 default: 00109 error("Vector: initialization must be either 0 or 1."); 00110 } 00111 } 00112 00113 HepVector::HepVector(int p, HepRandom &r) 00114 : m(p), nrow(p) 00115 { 00116 HepGenMatrix::mIter a = m.begin(); 00117 HepGenMatrix::mIter b = m.begin() + nrow; 00118 for(;a<b;a++) *a = r(); 00119 } 00120 00121 00122 // 00123 // Destructor 00124 // 00125 HepVector::~HepVector() { 00126 } 00127 00128 HepVector::HepVector(const HepVector &m1) 00129 : m(m1.nrow), nrow(m1.nrow) 00130 { 00131 m = m1.m; 00132 } 00133 00134 // 00135 // Copy constructor from the class of other precision 00136 // 00137 00138 00139 HepVector::HepVector(const HepMatrix &m1) 00140 : m(m1.nrow), nrow(m1.nrow) 00141 { 00142 if (m1.num_col() != 1) 00143 error("Vector::Vector(Matrix) : Matrix is not Nx1"); 00144 00145 m = m1.m; 00146 } 00147 00148 // trivial methods 00149 00150 inline int HepVector::num_row() const {return nrow;} 00151 inline int HepVector::num_size() const {return nrow;} 00152 inline int HepVector::num_col() const { return 1; } 00153 00154 // operator() 00155 00156 #ifdef MATRIX_BOUND_CHECK 00157 inline double & HepVector::operator()(int row, int col) 00158 { 00159 if( col!=1 || row<1 || row>nrow) 00160 error("Range error in HepVector::operator(i,j)"); 00161 #else 00162 inline double & HepVector::operator()(int row, int) 00163 { 00164 #endif 00165 00166 return *(m.begin()+(row-1)); 00167 } 00168 00169 #ifdef MATRIX_BOUND_CHECK 00170 inline const double & HepVector::operator()(int row, int col) const 00171 { 00172 if( col!=1 || row<1 || row>nrow) 00173 error("Range error in HepVector::operator(i,j)"); 00174 #else 00175 inline const double & HepVector::operator()(int row, int) const 00176 { 00177 #endif 00178 00179 return *(m.begin()+(row-1)); 00180 } 00181 00182 // Sub matrix 00183 00184 HepVector HepVector::sub(int min_row, int max_row) const 00185 #ifdef HEP_GNU_OPTIMIZED_RETURN 00186 return vret(max_row-min_row+1); 00187 { 00188 #else 00189 { 00190 HepVector vret(max_row-min_row+1); 00191 #endif 00192 if(max_row > num_row()) 00193 error("HepVector::sub: Index out of range"); 00194 HepGenMatrix::mIter a = vret.m.begin(); 00195 HepGenMatrix::mcIter b = m.begin() + min_row - 1; 00196 HepGenMatrix::mIter e = vret.m.begin() + vret.num_row(); 00197 for(;a<e;) *(a++) = *(b++); 00198 return vret; 00199 } 00200 00201 HepVector HepVector::sub(int min_row, int max_row) 00202 { 00203 HepVector vret(max_row-min_row+1); 00204 if(max_row > num_row()) 00205 error("HepVector::sub: Index out of range"); 00206 HepGenMatrix::mIter a = vret.m.begin(); 00207 HepGenMatrix::mIter b = m.begin() + min_row - 1; 00208 HepGenMatrix::mIter e = vret.m.begin() + vret.num_row(); 00209 for(;a<e;) *(a++) = *(b++); 00210 return vret; 00211 } 00212 00213 void HepVector::sub(int row,const HepVector &v1) 00214 { 00215 if(row <1 || row+v1.num_row()-1 > num_row()) 00216 error("HepVector::sub: Index out of range"); 00217 HepGenMatrix::mcIter a = v1.m.begin(); 00218 HepGenMatrix::mIter b = m.begin() + row - 1; 00219 HepGenMatrix::mcIter e = v1.m.begin() + v1.num_row(); 00220 for(;a<e;) *(b++) = *(a++); 00221 } 00222 00223 // 00224 // Direct sum of two matricies 00225 // 00226 00227 HepVector dsum(const HepVector &m1, 00228 const HepVector &m2) 00229 #ifdef HEP_GNU_OPTIMIZED_RETURN 00230 return mret(m1.num_row() + m2.num_row(), 0); 00231 { 00232 #else 00233 { 00234 HepVector mret(m1.num_row() + m2.num_row(), 00235 0); 00236 #endif 00237 mret.sub(1,m1); 00238 mret.sub(m1.num_row()+1,m2); 00239 return mret; 00240 } 00241 00242 /* ----------------------------------------------------------------------- 00243 This section contains support routines for matrix.h. This section contains 00244 The two argument functions +,-. They call the copy constructor and +=,-=. 00245 ----------------------------------------------------------------------- */ 00246 HepVector HepVector::operator- () const 00247 #ifdef HEP_GNU_OPTIMIZED_RETURN 00248 return m2(nrow); 00249 { 00250 #else 00251 { 00252 HepVector m2(nrow); 00253 #endif 00254 register HepGenMatrix::mcIter a=m.begin(); 00255 register HepGenMatrix::mIter b=m2.m.begin(); 00256 register HepGenMatrix::mcIter e=m.begin()+num_size(); 00257 for(;a<e; a++, b++) (*b) = -(*a); 00258 return m2; 00259 } 00260 00261 00262 00263 HepVector operator+(const HepMatrix &m1,const HepVector &m2) 00264 #ifdef HEP_GNU_OPTIMIZED_RETURN 00265 return mret(m2); 00266 { 00267 #else 00268 { 00269 HepVector mret(m2); 00270 #endif 00271 CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),1,+); 00272 mret += m1; 00273 return mret; 00274 } 00275 00276 HepVector operator+(const HepVector &m1,const HepMatrix &m2) 00277 #ifdef HEP_GNU_OPTIMIZED_RETURN 00278 return mret(m1); 00279 { 00280 #else 00281 { 00282 HepVector mret(m1); 00283 #endif 00284 CHK_DIM_2(m1.num_row(),m2.num_row(),1,m2.num_col(),+); 00285 mret += m2; 00286 return mret; 00287 } 00288 00289 HepVector operator+(const HepVector &m1,const HepVector &m2) 00290 #ifdef HEP_GNU_OPTIMIZED_RETURN 00291 return mret(m1.num_row()); 00292 { 00293 #else 00294 { 00295 HepVector mret(m1.num_row()); 00296 #endif 00297 CHK_DIM_1(m1.num_row(),m2.num_row(),+); 00298 SIMPLE_TOP(+) 00299 return mret; 00300 } 00301 00302 // 00303 // operator - 00304 // 00305 00306 HepVector operator-(const HepMatrix &m1,const HepVector &m2) 00307 #ifdef HEP_GNU_OPTIMIZED_RETURN 00308 return mret; 00309 { 00310 #else 00311 { 00312 HepVector mret; 00313 #endif 00314 CHK_DIM_2(m1.num_row(),m2.num_row(),m1.num_col(),1,-); 00315 mret = m1; 00316 mret -= m2; 00317 return mret; 00318 } 00319 00320 HepVector operator-(const HepVector &m1,const HepMatrix &m2) 00321 #ifdef HEP_GNU_OPTIMIZED_RETURN 00322 return mret(m1); 00323 { 00324 #else 00325 { 00326 HepVector mret(m1); 00327 #endif 00328 CHK_DIM_2(m1.num_row(),m2.num_row(),1,m2.num_col(),-); 00329 mret -= m2; 00330 return mret; 00331 } 00332 00333 HepVector operator-(const HepVector &m1,const HepVector &m2) 00334 #ifdef HEP_GNU_OPTIMIZED_RETURN 00335 return mret(m1.num_row()); 00336 { 00337 #else 00338 { 00339 HepVector mret(m1.num_row()); 00340 #endif 00341 CHK_DIM_1(m1.num_row(),m2.num_row(),-); 00342 SIMPLE_TOP(-) 00343 return mret; 00344 } 00345 00346 /* ----------------------------------------------------------------------- 00347 This section contains support routines for matrix.h. This file contains 00348 The two argument functions *,/. They call copy constructor and then /=,*=. 00349 ----------------------------------------------------------------------- */ 00350 00351 HepVector operator/( 00352 const HepVector &m1,double t) 00353 #ifdef HEP_GNU_OPTIMIZED_RETURN 00354 return mret(m1); 00355 { 00356 #else 00357 { 00358 HepVector mret(m1); 00359 #endif 00360 mret /= t; 00361 return mret; 00362 } 00363 00364 HepVector operator*(const HepVector &m1,double t) 00365 #ifdef HEP_GNU_OPTIMIZED_RETURN 00366 return mret(m1); 00367 { 00368 #else 00369 { 00370 HepVector mret(m1); 00371 #endif 00372 mret *= t; 00373 return mret; 00374 } 00375 00376 HepVector operator*(double t,const HepVector &m1) 00377 #ifdef HEP_GNU_OPTIMIZED_RETURN 00378 return mret(m1); 00379 { 00380 #else 00381 { 00382 HepVector mret(m1); 00383 #endif 00384 mret *= t; 00385 return mret; 00386 } 00387 00388 HepVector operator*(const HepMatrix &m1,const HepVector &m2) 00389 #ifdef HEP_GNU_OPTIMIZED_RETURN 00390 return mret(m1.num_row()); 00391 { 00392 #else 00393 { 00394 HepVector mret(m1.num_row()); 00395 #endif 00396 CHK_DIM_1(m1.num_col(),m2.num_row(),*); 00397 HepGenMatrix::mcIter m1p,m2p,vp; 00398 HepGenMatrix::mIter m3p; 00399 double temp; 00400 m3p=mret.m.begin(); 00401 for(m1p=m1.m.begin();m1p<m1.m.begin()+m1.num_row()*m1.num_col();m1p=m2p) 00402 { 00403 temp=0; 00404 vp=m2.m.begin(); 00405 m2p=m1p; 00406 while(m2p<m1p+m1.num_col()) 00407 temp+=(*(m2p++))*(*(vp++)); 00408 *(m3p++)=temp; 00409 } 00410 return mret; 00411 } 00412 00413 HepMatrix operator*(const HepVector &m1,const HepMatrix &m2) 00414 #ifdef HEP_GNU_OPTIMIZED_RETURN 00415 return mret(m1.num_row(),m2.num_col()); 00416 { 00417 #else 00418 { 00419 HepMatrix mret(m1.num_row(),m2.num_col()); 00420 #endif 00421 CHK_DIM_1(1,m2.num_row(),*); 00422 HepGenMatrix::mcIter m1p; 00423 HepMatrix::mcIter m2p; 00424 HepMatrix::mIter mrp=mret.m.begin(); 00425 for(m1p=m1.m.begin();m1p<m1.m.begin()+m1.num_row();m1p++) 00426 for(m2p=m2.m.begin();m2p<m2.m.begin()+m2.num_col();m2p++) 00427 *(mrp++)=*m1p*(*m2p); 00428 return mret; 00429 } 00430 00431 /* ----------------------------------------------------------------------- 00432 This section contains the assignment and inplace operators =,+=,-=,*=,/=. 00433 ----------------------------------------------------------------------- */ 00434 00435 HepMatrix & HepMatrix::operator+=(const HepVector &m2) 00436 { 00437 CHK_DIM_2(num_row(),m2.num_row(),num_col(),1,+=); 00438 SIMPLE_BOP(+=) 00439 return (*this); 00440 } 00441 00442 HepVector & HepVector::operator+=(const HepMatrix &m2) 00443 { 00444 CHK_DIM_2(num_row(),m2.num_row(),1,m2.num_col(),+=); 00445 SIMPLE_BOP(+=) 00446 return (*this); 00447 } 00448 00449 HepVector & HepVector::operator+=(const HepVector &m2) 00450 { 00451 CHK_DIM_1(num_row(),m2.num_row(),+=); 00452 SIMPLE_BOP(+=) 00453 return (*this); 00454 } 00455 00456 HepMatrix & HepMatrix::operator-=(const HepVector &m2) 00457 { 00458 CHK_DIM_2(num_row(),m2.num_row(),num_col(),1,-=); 00459 SIMPLE_BOP(-=) 00460 return (*this); 00461 } 00462 00463 HepVector & HepVector::operator-=(const HepMatrix &m2) 00464 { 00465 CHK_DIM_2(num_row(),m2.num_row(),1,m2.num_col(),-=); 00466 SIMPLE_BOP(-=) 00467 return (*this); 00468 } 00469 00470 HepVector & HepVector::operator-=(const HepVector &m2) 00471 { 00472 CHK_DIM_1(num_row(),m2.num_row(),-=); 00473 SIMPLE_BOP(-=) 00474 return (*this); 00475 } 00476 00477 HepVector & HepVector::operator/=(double t) 00478 { 00479 SIMPLE_UOP(/=) 00480 return (*this); 00481 } 00482 00483 HepVector & HepVector::operator*=(double t) 00484 { 00485 SIMPLE_UOP(*=) 00486 return (*this); 00487 } 00488 00489 HepMatrix & HepMatrix::operator=(const HepVector &m1) 00490 { 00491 if(m1.nrow != size) 00492 { 00493 size = m1.nrow; 00494 m.resize(size); 00495 } 00496 nrow = m1.nrow; 00497 ncol = 1; 00498 m = m1.m; 00499 return (*this); 00500 } 00501 00502 HepVector & HepVector::operator=(const HepVector &m1) 00503 { 00504 if(m1.nrow != nrow) 00505 { 00506 nrow = m1.nrow; 00507 m.resize(nrow); 00508 } 00509 m = m1.m; 00510 return (*this); 00511 } 00512 00513 HepVector & HepVector::operator=(const HepMatrix &m1) 00514 { 00515 if (m1.num_col() != 1) 00516 error("Vector::operator=(Matrix) : Matrix is not Nx1"); 00517 00518 if(m1.nrow != nrow) 00519 { 00520 nrow = m1.nrow; 00521 m.resize(nrow); 00522 } 00523 m = m1.m; 00524 return (*this); 00525 } 00526 00527 HepVector & HepVector::operator=(const Hep3Vector &v) 00528 { 00529 if(nrow != 3) 00530 { 00531 nrow = 3; 00532 m.resize(nrow); 00533 } 00534 m[0] = v.x(); 00535 m[1] = v.y(); 00536 m[2] = v.z(); 00537 return (*this); 00538 } 00539 00540 // 00541 // Copy constructor from the class of other precision 00542 // 00543 00544 00545 // Print the Matrix. 00546 00547 std::ostream& operator<<(std::ostream &s, const HepVector &q) 00548 { 00549 s << std::endl; 00550 /* Fixed format needs 3 extra characters for field, while scientific needs 7 */ 00551 int width; 00552 if(s.flags() & std::ios::fixed) 00553 width = s.precision()+3; 00554 else 00555 width = s.precision()+7; 00556 for(int irow = 1; irow<= q.num_row(); irow++) 00557 { 00558 s.width(width); 00559 s << q(irow) << std::endl; 00560 } 00561 return s; 00562 } 00563 00564 HepMatrix HepVector::T() const 00565 #ifdef HEP_GNU_OPTIMIZED_RETURN 00566 return mret(1,num_row()); 00567 { 00568 #else 00569 { 00570 HepMatrix mret(1,num_row()); 00571 #endif 00572 mret.m = m; 00573 return mret; 00574 } 00575 00576 double dot(const HepVector &v1,const HepVector &v2) 00577 { 00578 if(v1.num_row()!=v2.num_row()) 00579 HepGenMatrix::error("v1 and v2 need to be the same size in dot(HepVector, HepVector)"); 00580 double d= 0; 00581 HepGenMatrix::mcIter a = v1.m.begin(); 00582 HepGenMatrix::mcIter b = v2.m.begin(); 00583 HepGenMatrix::mcIter e = a + v1.num_size(); 00584 for(;a<e;) d += (*(a++)) * (*(b++)); 00585 return d; 00586 } 00587 00588 HepVector HepVector:: 00589 apply(double (*f)(double, int)) const 00590 #ifdef HEP_GNU_OPTIMIZED_RETURN 00591 return mret(num_row()); 00592 { 00593 #else 00594 { 00595 HepVector mret(num_row()); 00596 #endif 00597 HepGenMatrix::mcIter a = m.begin(); 00598 HepGenMatrix::mIter b = mret.m.begin(); 00599 for(int ir=1;ir<=num_row();ir++) { 00600 *(b++) = (*f)(*(a++), ir); 00601 } 00602 return mret; 00603 } 00604 00605 void HepVector::invert(int &) { 00606 error("HepVector::invert: You cannot invert a Vector"); 00607 } 00608 00609 HepVector solve(const HepMatrix &a, const HepVector &v) 00610 #ifdef HEP_GNU_OPTIMIZED_RETURN 00611 return vret(v); 00612 { 00613 #else 00614 { 00615 HepVector vret(v); 00616 #endif 00617 static int max_array = 20; 00618 static int *ir = new int [max_array+1]; 00619 00620 if(a.ncol != a.nrow) 00621 HepGenMatrix::error("Matrix::solve Matrix is not NxN"); 00622 if(a.ncol != v.nrow) 00623 HepGenMatrix::error("Matrix::solve Vector has wrong number of rows"); 00624 00625 int n = a.ncol; 00626 if (n > max_array) { 00627 delete [] ir; 00628 max_array = n; 00629 ir = new int [max_array+1]; 00630 } 00631 double det; 00632 HepMatrix mt(a); 00633 int i = mt.dfact_matrix(det, ir); 00634 if (i!=0) { 00635 for (i=1;i<=n;i++) vret(i) = 0; 00636 return vret; 00637 } 00638 double s21, s22; 00639 int nxch = ir[n]; 00640 if (nxch!=0) { 00641 for (int mm=1;mm<=nxch;mm++) { 00642 int ij = ir[mm]; 00643 i = ij >> 12; 00644 int j = ij%4096; 00645 double te = vret(i); 00646 vret(i) = vret(j); 00647 vret(j) = te; 00648 } 00649 } 00650 vret(1) = mt(1,1) * vret(1); 00651 if (n!=1) { 00652 for (i=2;i<=n;i++) { 00653 s21 = -vret(i); 00654 for (int j=1;j<i;j++) { 00655 s21 += mt(i,j) * vret(j); 00656 } 00657 vret(i) = -mt(i,i)*s21; 00658 } 00659 for (i=1;i<n;i++) { 00660 int nmi = n-i; 00661 s22 = -vret(nmi); 00662 for (int j=1;j<=i;j++) { 00663 s22 += mt(nmi,n-j+1) * vret(n-j+1); 00664 } 00665 vret(nmi) = -s22; 00666 } 00667 } 00668 return vret; 00669 } 00670 00671 } // namespace CLHEP