16#include "CLHEP/Matrix/defs.h"
17#include "CLHEP/Random/Random.h"
18#include "CLHEP/Matrix/Matrix.h"
19#include "CLHEP/Matrix/SymMatrix.h"
20#include "CLHEP/Matrix/DiagMatrix.h"
21#include "CLHEP/Matrix/Vector.h"
23#ifdef HEP_DEBUG_INLINE
24#include "CLHEP/Matrix/Matrix.icc"
31#define SIMPLE_UOP(OPER) \
34 for(;a!=e; a++) (*a) OPER t;
36#define SIMPLE_BOP(OPER) \
37 HepMatrix::mIter a=m.begin(); \
38 HepMatrix::mcIter b=hm2.m.begin(); \
39 HepMatrix::mIter e=m.end(); \
40 for(;a!=e; a++, b++) (*a) OPER (*b);
42#define SIMPLE_TOP(OPER) \
43 HepMatrix::mcIter a=hm1.m.begin(); \
44 HepMatrix::mcIter b=hm2.m.begin(); \
45 HepMatrix::mIter t=mret.m.begin(); \
46 HepMatrix::mcIter e=hm1.m.end(); \
47 for(;a!=e; a++, b++, t++) (*t) = (*a) OPER (*b);
51#define CHK_DIM_2(r1,r2,c1,c2,fun) \
52 if (r1!=r2 || c1!=c2) { \
53 HepGenMatrix::error("Range error in Matrix function " #fun "(1)."); \
56#define CHK_DIM_1(c1,r2,fun) \
58 HepGenMatrix::error("Range error in Matrix function " #fun "(2)."); \
64 : m(p*q), nrow(p), ncol(q)
70 : m(p*q), nrow(p), ncol(q)
84 for(
int step=0; step < size_; step+=(ncol+1) ) *(
a+step) = 1.0;
86 error(
"Invalid dimension in HepMatrix(int,int,1).");
91 error(
"Matrix: initialization must be either 0 or 1.");
97 : m(p*q), nrow(p), ncol(q)
103 for(;
a<
b;
a++) *
a = r();
112 :
HepGenMatrix(hm1), m(hm1.size_), nrow(hm1.nrow), ncol(hm1.ncol), size_(hm1.size_)
130#ifdef MATRIX_BOUND_CHECK
132 error(
"Range error in HepMatrix::operator()");
134 return *(m.begin()+(row-1)*ncol+col-1);
139#ifdef MATRIX_BOUND_CHECK
141 error(
"Range error in HepMatrix::operator()");
143 return *(m.begin()+(row-1)*ncol+col-1);
148 : m(hm1.nrow*hm1.nrow), nrow(hm1.nrow), ncol(hm1.nrow)
152 mcIter sjk = hm1.m.begin();
154 for(
int j=0; j!=nrow; ++j) {
155 for(
int k=0; k<=j; ++k) {
160 if(k!=j) m[k*nrow+j] = *sjk;
167 : m(hm1.nrow*hm1.nrow), nrow(hm1.nrow), ncol(hm1.nrow)
173 mcIter mr = hm1.m.begin();
174 for(
int r=0;r<n;r++) {
175 mrr = m.begin()+(n+1)*r;
181 : m(hm1.nrow), nrow(hm1.nrow), ncol(1)
196 int min_col,
int max_col)
const
197#ifdef HEP_GNU_OPTIMIZED_RETURN
198return mret(max_row-min_row+1,max_col-min_col+1);
202 HepMatrix mret(max_row-min_row+1,max_col-min_col+1);
204 if(max_row > num_row() || max_col >num_col())
205 error(
"HepMatrix::sub: Index out of range");
208 mcIter b1 = m.begin() + (min_row - 1) * nc + min_col - 1;
209 int rowsize = mret.num_row();
210 for(
int irow=1; irow<=rowsize; ++irow) {
212 for(
int icol=0; icol<mret.num_col(); ++icol) {
215 if(irow<rowsize) b1 += nc;
224 error(
"HepMatrix::sub: Index out of range");
227 mIter b1 = m.begin() + (row - 1) * nc + col - 1;
229 for(
int irow=1; irow<=rowsize; ++irow) {
231 for(
int icol=0; icol<hm1.
num_col(); ++icol) {
234 if(irow<rowsize) b1 += nc;
243#ifdef HEP_GNU_OPTIMIZED_RETURN
262#ifdef HEP_GNU_OPTIMIZED_RETURN
263 return hm2(nrow, ncol);
272 for(;
a<e;
a++,
b++) (*
b) = -(*a);
279#ifdef HEP_GNU_OPTIMIZED_RETURN
280 return mret(hm1.nrow, hm1.ncol);
296#ifdef HEP_GNU_OPTIMIZED_RETURN
297 return mret(hm1.num_row(), hm1.num_col());
301 HepMatrix mret(hm1.num_row(), hm1.num_col());
304 hm1.num_col(),hm2.
num_col(),-);
316#ifdef HEP_GNU_OPTIMIZED_RETURN
328#ifdef HEP_GNU_OPTIMIZED_RETURN
340#ifdef HEP_GNU_OPTIMIZED_RETURN
352#ifdef HEP_GNU_OPTIMIZED_RETURN
353 return mret(hm1.nrow,hm2.ncol,0);
362 int m1cols = hm1.ncol;
363 int m2cols = hm2.ncol;
365 for (
int i=0; i<hm1.nrow; i++)
367 for (
int j=0; j<m1cols; j++)
369 register double temp = hm1.m[i*m1cols+j];
377 (*pt) += temp * (*pb);
419 if(hm1.nrow*hm1.ncol != size_)
421 size_ = hm1.nrow * hm1.ncol;
440 if(os.flags() & std::ios::fixed)
441 width = os.precision()+3;
443 width = os.precision()+7;
444 for(
int irow = 1; irow<= q.
num_row(); irow++)
446 for(
int icol = 1; icol <= q.
num_col(); icol++)
449 os << q(irow,icol) <<
" ";
457#ifdef HEP_GNU_OPTIMIZED_RETURN
458return mret(ncol,nrow);
465 mIter pt = mret.m.begin();
466 for(
int nr=0; nr<nrow; ++nr) {
467 for(
int nc=0; nc<ncol; ++nc) {
468 pt = mret.m.begin() + nr + nrow*nc;
477#ifdef HEP_GNU_OPTIMIZED_RETURN
478return mret(num_row(),num_col());
486 for(
int ir=1;ir<=num_row();ir++) {
487 for(
int ic=1;ic<=num_col();ic++) {
488 *(
b++) = (*
f)(*(
a++), ir, ic);
494int HepMatrix::dfinv_matrix(
int *ir) {
496 error(
"dfinv_matrix: Matrix is not NxN");
501 register double s33, s34;
503 mIter hm11 = m.begin();
504 mIter hm12 = hm11 + 1;
507 *hm21 = -(*hm22) * (*hm11) * (*hm21);
510 mIter mimim = hm11 +
n + 1;
511 for (
int i=3;i<=
n;i++) {
513 mIter mi = hm11 + (i-1) * n;
514 mIter mii= hm11 + (i-1) * n + i - 1;
517 mIter mji = mj + i - 1;
519 for (
int j=1;j<=ihm2;j++) {
522 mIter mkj = mj + j - 1;
523 mIter mik = mi + j - 1;
525 mIter mkpi = mj +
n + i - 1;
526 for (
int k=j;k<=ihm2;k++) {
527 s31 += (*mkj) * (*(mik++));
528 s32 += (*(mjkp++)) * (*mkpi);
532 *mij = -(*mii) * (((*(mij-
n)))*( (*(mii-1)))+(s31));
538 *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
539 *(mimim+1) = -(*(mimim+1));
545 for (
int i=1;i<
n;i++) {
549 for (j=1; j<=i;j++) {
552 mIter mikj = mi + j - 1;
553 mIter miik = mii + 1;
555 for (;miik<min_end;) {
558 s33 += (*mikj) * (*(miik++));
562 for (j=1;j<=ni;j++) {
564 mIter miik = mii + j;
565 for (
int k=j;k<=ni;k++) {
567 mIter mikij = mii + k *
n + j;
568 s34 += *mikij * (*(miik++));
576 if (nxch==0)
return 0;
577 for (
int hmm=1;hmm<=nxch;hmm++) {
578 int k = nxch - hmm + 1;
582 for (k=1; k<=
n;k++) {
584 mIter mki = hm11 + (k-1)*n + i - 1;
585 mIter mkj = hm11 + (k-1)*n + j - 1;
596int HepMatrix::dfact_matrix(
double &det,
int *ir) {
598 error(
"dfact_matrix: Matrix is not NxN");
604 double g1 = 1.0e-19, g2 = 1.0e19;
609 double epsilon = 8*DBL_EPSILON;
614 int normal = 0, imposs = -1;
615 int jrange = 0, jover = 1, junder = -1;
620 mIter mj = m.begin();
622 for (
int j=1;j<=
n;j++) {
627 for (
int i=j+1;i<
n;i++) {
628 q = (fabs(*(mj + n*(i-j) + j - 1)));
645 mIter mkl = m.begin() + (k-1)*n;
646 for (
int l=1;l<=
n;l++) {
652 ir[nxch] = (((j)<<12)+(k));
666 if (jfail == jrange) jfail = junder;
669 if (jfail==jrange) jfail = jover;
674 for (k=j+1;k<=
n;k++) {
680 mIter mik = m.begin() + k - 1;
681 mIter mijp = m.begin() + j;
684 for (
int i=1;i<j;i++) {
685 s11 += (*mik) * (*(mji++));
686 s12 += (*mijp) * (*(mki++));
691 *(mjk++) = -s11 * (*mjj);
692 *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
701 if (nxch%2==1) det = -det;
702 if (jfail !=jrange) det = 0.0;
709 error(
"HepMatrix::invert: Matrix is not NxN");
711 static int max_array = 20;
712 static int *ir =
new int [max_array+1];
714 if (ncol > max_array) {
717 ir =
new int [max_array+1];
720 double det, temp, sd;
724 double c11,c12,c13,c21,c22,c23,c31,c32,c33;
726 c11 = (*(m.begin()+4)) * (*(m.begin()+8)) - (*(m.begin()+5)) * (*(m.begin()+7));
727 c12 = (*(m.begin()+5)) * (*(m.begin()+6)) - (*(m.begin()+3)) * (*(m.begin()+8));
728 c13 = (*(m.begin()+3)) * (*(m.begin()+7)) - (*(m.begin()+4)) * (*(m.begin()+6));
729 c21 = (*(m.begin()+7)) * (*(m.begin()+2)) - (*(m.begin()+8)) * (*(m.begin()+1));
730 c22 = (*(m.begin()+8)) * (*m.begin()) - (*(m.begin()+6)) * (*(m.begin()+2));
731 c23 = (*(m.begin()+6)) * (*(m.begin()+1)) - (*(m.begin()+7)) * (*m.begin());
732 c31 = (*(m.begin()+1)) * (*(m.begin()+5)) - (*(m.begin()+2)) * (*(m.begin()+4));
733 c32 = (*(m.begin()+2)) * (*(m.begin()+3)) - (*m.begin()) * (*(m.begin()+5));
734 c33 = (*m.begin()) * (*(m.begin()+4)) - (*(m.begin()+1)) * (*(m.begin()+3));
735 t1 = fabs(*m.begin());
736 t2 = fabs(*(m.begin()+3));
737 t3 = fabs(*(m.begin()+6));
740 temp = *(m.begin()+6);
741 det = c23*c12-c22*c13;
744 det = c22*c33-c23*c32;
746 }
else if (t3 >= t2) {
747 temp = *(m.begin()+6);
748 det = c23*c12-c22*c13;
750 temp = *(m.begin()+3);
751 det = c13*c32-c12*c33;
758 double s1 = temp/det;
759 mIter hmm = m.begin();
773 det = (*m.begin())*(*(m.begin()+3)) - (*(m.begin()+1))*(*(m.begin()+2));
779 temp = sd*(*(m.begin()+3));
780 *(m.begin()+1) *= -sd;
781 *(m.begin()+2) *= -sd;
782 *(m.begin()+3) = sd*(*m.begin());
787 if ((*(m.begin()))==0) {
791 *(m.begin()) = 1.0/(*(m.begin()));
803 ifail = dfact_matrix(det, ir);
816 static int max_array = 20;
817 static int *ir =
new int [max_array+1];
819 error(
"HepMatrix::determinant: Matrix is not NxN");
820 if (ncol > max_array) {
823 ir =
new int [max_array+1];
827 int i = mt.dfact_matrix(det, ir);
834 for (
mcIter d = m.begin(); d < m.end(); d += (ncol+1) )
#define CHK_DIM_2(r1, r2, c1, c2, fun)
#define CHK_DIM_1(c1, r2, fun)
std::vector< double, Alloc< double, 25 > >::const_iterator mcIter
std::vector< double, Alloc< double, 25 > >::iterator mIter
static void error(const char *s)
virtual void invertHaywood4(int &ierr)
HepMatrix sub(int min_row, int max_row, int min_col, int max_col) const
HepMatrix & operator=(const HepMatrix &)
HepMatrix apply(double(*f)(double, int, int)) const
virtual int num_col() const
virtual const double & operator()(int row, int col) const
HepMatrix & operator/=(double t)
virtual int num_row() const
HepMatrix & operator+=(const HepMatrix &)
HepMatrix operator-() const
HepMatrix & operator*=(double t)
HepMatrix & operator-=(const HepMatrix &)
double determinant() const
virtual void invertHaywood6(int &ierr)
virtual int num_size() const
virtual void invertHaywood5(int &ierr)
Hep3Vector operator+(const Hep3Vector &, const Hep3Vector &)
HepLorentzRotation operator*(const HepRotation &r, const HepLorentzRotation <)
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2)