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