14 #include "CLHEP/Matrix/defs.h"
15 #include "CLHEP/Random/Random.h"
16 #include "CLHEP/Matrix/DiagMatrix.h"
17 #include "CLHEP/Matrix/Matrix.h"
18 #include "CLHEP/Matrix/SymMatrix.h"
19 #include "CLHEP/Matrix/Vector.h"
21 #ifdef HEP_DEBUG_INLINE
22 #include "CLHEP/Matrix/DiagMatrix.icc"
29 #define SIMPLE_UOP(OPER) \
30 HepMatrix::mIter a=m.begin(); \
31 HepMatrix::mIter e=m.begin()+num_size(); \
32 for(;a<e; a++) (*a) OPER t;
34 #define SIMPLE_BOP(OPER) \
35 HepMatrix::mIter a=m.begin(); \
36 HepMatrix::mcIter b=m2.m.begin(); \
37 HepMatrix::mIter e=m.begin()+num_size(); \
38 for(;a<e; a++, b++) (*a) OPER (*b);
40 #define SIMPLE_TOP(OPER) \
41 HepMatrix::mcIter a=m1.m.begin(); \
42 HepMatrix::mcIter b=m2.m.begin(); \
43 HepMatrix::mIter t=mret.m.begin(); \
44 HepMatrix::mcIter e=m1.m.begin()+m1.nrow; \
45 for( ;a<e; a++, b++, t++) (*t) = (*a) OPER (*b);
47 #define CHK_DIM_2(r1,r2,c1,c2,fun) \
48 if (r1!=r2 || c1!=c2) { \
49 HepGenMatrix::error("Range error in DiagMatrix function " #fun "(1)."); \
52 #define CHK_DIM_1(c1,r2,fun) \
54 HepGenMatrix::error("Range error in DiagMatrix function " #fun "(2)."); \
59 #if defined(__sun) || !defined(__GNUG__)
63 double HepDiagMatrix::zero = 0;
65 const double HepDiagMatrix::zero = 0;
88 for( ; a<
b; a++) *a = 1.0;
92 error(
"DiagMatrix: initialization must be either 0 or 1.");
101 for(;a<
b;a++) *a = r();
122 #ifdef HEP_GNU_OPTIMIZED_RETURN
123 return mret(max_row-min_row+1);
129 if(max_row > num_row())
130 error(
"HepDiagMatrix::sub: Index out of range");
134 for(;a<e;) *(a++) = *(b++);
142 error(
"HepDiagMatrix::sub: Index out of range");
146 for(;a<e;) *(a++) = *(b++);
153 error(
"HepDiagMatrix::sub: Index out of range");
157 for(;a<e;) *(b++) = *(a++);
166 #ifdef HEP_GNU_OPTIMIZED_RETURN
180 #ifdef HEP_GNU_OPTIMIZED_RETURN
190 for(;a<e; a++, b++) (*b) = -(*a);
197 #ifdef HEP_GNU_OPTIMIZED_RETURN
211 #ifdef HEP_GNU_OPTIMIZED_RETURN
225 #ifdef HEP_GNU_OPTIMIZED_RETURN
226 return mret(m1.nrow);
238 #ifdef HEP_GNU_OPTIMIZED_RETURN
251 #ifdef HEP_GNU_OPTIMIZED_RETURN
268 #ifdef HEP_GNU_OPTIMIZED_RETURN
281 #ifdef HEP_GNU_OPTIMIZED_RETURN
295 #ifdef HEP_GNU_OPTIMIZED_RETURN
296 return mret(m1.nrow);
307 #ifdef HEP_GNU_OPTIMIZED_RETURN
320 #ifdef HEP_GNU_OPTIMIZED_RETURN
340 #ifdef HEP_GNU_OPTIMIZED_RETURN
352 #ifdef HEP_GNU_OPTIMIZED_RETURN
364 #ifdef HEP_GNU_OPTIMIZED_RETURN
376 #ifdef HEP_GNU_OPTIMIZED_RETURN
386 for(
int irow=1;irow<=m1.
num_row();irow++) {
388 for(
int icol=1;icol<=m1.
num_col();icol++) {
389 *(mir++) = *(mit1++) * (*(mcc++));
396 #ifdef HEP_GNU_OPTIMIZED_RETURN
407 for(
int irow=1;irow<=m2.
num_row();irow++) {
408 for(
int icol=1;icol<=m2.
num_col();icol++) {
409 *(mir++) = *(mit1++) * (*mrr);
417 #ifdef HEP_GNU_OPTIMIZED_RETURN
429 for(;a<e;) *(a++) = *(b++) * (*(c++));
434 #ifdef HEP_GNU_OPTIMIZED_RETURN
444 for(
int icol=1;icol<=m1.
num_col();icol++) {
445 *(mir++) = *(mi1++) * *(mi2++);
458 mIter mrr = m.begin();
460 for(
int r=1;r<=
n;r++) {
462 if(r<n) mrr += (n+1);
472 for(
int i=1;i<=
num_row();i++) {
490 mIter mrr = m.begin();
492 for(
int r=1;r<=
n;r++) {
494 if(r<n) mrr += (n+1);
504 for(
int i=1;i<=
num_row();i++) {
532 if(m1.nrow*m1.nrow != size_)
534 size_ = m1.nrow * m1.nrow;
541 mIter mrr = m.begin();
543 for(
int r=1;r<=
n;r++) {
545 if(r<n) mrr += (n+1);
568 if(s.flags() & std::ios::fixed)
569 width = s.precision()+3;
571 width = s.precision()+7;
572 for(
int irow = 1; irow<= q.
num_row(); irow++)
574 for(
int icol = 1; icol <= q.
num_col(); icol++)
577 s << q(irow,icol) <<
" ";
585 apply(
double (*
f)(
double,
int,
int)) const
586 #ifdef HEP_GNU_OPTIMIZED_RETURN
587 return mret(num_row());
595 for(
int ir=1;ir<=num_row();ir++) {
596 *(b++) = (*
f)(*(a++), ir, ir);
610 for(
int r=1;r<=nrow;r++) {
612 if(r<nrow) a += (nrow+1);
625 for(
int r=1;r<=nrow;r++) {
627 if(r<nrow) a += (r+1);
632 #ifdef HEP_GNU_OPTIMIZED_RETURN
644 for(
int r=1;r<=mret.num_row();r++) {
647 for(
int c=1;c<=r;c++) {
649 register double tmp = 0;
651 for(
int i=0;i<m1.
num_col();i++)
652 tmp+=*(mr++) * *(mc++) * *(mi++);
661 register double mret;
665 mret = *(mv)* *(mv)* *(mi++);
667 for(
int i=2;i<=m1.
num_row();i++) {
668 mret+=*(mv)* *(mv)* *(mi++);
675 #ifdef HEP_GNU_OPTIMIZED_RETURN
686 for(
int r=1;r<=mret.num_row();r++)
687 for(
int c=1;c<=r;c++)
690 register double tmp = m1(1,r)*m1(1,c)* *(mi++);
691 for(
int i=2;i<=m1.
num_row();i++)
692 tmp+=m1(i,r)*m1(i,c)* *(mi++);
693 mret.
fast(r,c) = tmp;
704 if(*(mm++)==0)
return;