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