CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

Matrix.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 // ---------------------------------------------------------------------------
3 //
4 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
5 //
6 
7 #ifdef GNUPRAGMA
8 #pragma implementation
9 #endif
10 
11 #include <string.h>
12 #include <float.h> // for DBL_EPSILON
13 #include <cmath>
14 #include <stdlib.h>
15 
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"
22 
23 #ifdef HEP_DEBUG_INLINE
24 #include "CLHEP/Matrix/Matrix.icc"
25 #endif
26 
27 namespace CLHEP {
28 
29 // Simple operation for all elements
30 
31 #define SIMPLE_UOP(OPER) \
32  mIter a=m.begin(); \
33  mIter e=m.end(); \
34  for(;a!=e; a++) (*a) OPER t;
35 
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);
41 
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);
48 
49 // Static functions.
50 
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)."); \
54  }
55 
56 #define CHK_DIM_1(c1,r2,fun) \
57  if (c1!=r2) { \
58  HepGenMatrix::error("Range error in Matrix function " #fun "(2)."); \
59  }
60 
61 // Constructors. (Default constructors are inlined and in .icc file)
62 
64  : m(p*q), nrow(p), ncol(q)
65 {
66  size_ = nrow * ncol;
67 }
68 
69 HepMatrix::HepMatrix(int p,int q,int init)
70  : m(p*q), nrow(p), ncol(q)
71 {
72  size_ = nrow * ncol;
73 
74  if (size_ > 0) {
75  switch(init)
76  {
77  case 0:
78  break;
79 
80  case 1:
81  {
82  if ( ncol == nrow ) {
83  mIter a = m.begin();
84  for( int step=0; step < size_; step+=(ncol+1) ) *(a+step) = 1.0;
85  } else {
86  error("Invalid dimension in HepMatrix(int,int,1).");
87  }
88  break;
89  }
90  default:
91  error("Matrix: initialization must be either 0 or 1.");
92  }
93  }
94 }
95 
97  : m(p*q), nrow(p), ncol(q)
98 {
99  size_ = nrow * ncol;
100 
101  mIter a = m.begin();
102  mIter b = m.end();
103  for(; a<b; a++) *a = r();
104 }
105 //
106 // Destructor
107 //
109 }
110 
112  : HepGenMatrix(hm1), m(hm1.size_), nrow(hm1.nrow), ncol(hm1.ncol), size_(hm1.size_)
113 {
114  m = hm1.m;
115 
116 }
117 
118 // trivial
119 
120 int HepMatrix::num_row() const { return nrow;}
121 
122 int HepMatrix::num_col() const { return ncol;}
123 
124 int HepMatrix::num_size() const { return size_;}
125 
126 // operator()
127 
128 double & HepMatrix::operator()(int row, int col)
129 {
130 #ifdef MATRIX_BOUND_CHECK
131  if(row<1 || row>num_row() || col<1 || col>num_col())
132  error("Range error in HepMatrix::operator()");
133 #endif
134  return *(m.begin()+(row-1)*ncol+col-1);
135 }
136 
137 const double & HepMatrix::operator()(int row, int col) const
138 {
139 #ifdef MATRIX_BOUND_CHECK
140  if(row<1 || row>num_row() || col<1 || col>num_col())
141  error("Range error in HepMatrix::operator()");
142 #endif
143  return *(m.begin()+(row-1)*ncol+col-1);
144 }
145 
146 
148  : m(hm1.nrow*hm1.nrow), nrow(hm1.nrow), ncol(hm1.nrow)
149 {
150  size_ = nrow * ncol;
151 
152  mcIter sjk = hm1.m.begin();
153  // j >= k
154  for(int j=0; j!=nrow; ++j) {
155  for(int k=0; k<=j; ++k) {
156  m[j*ncol+k] = *sjk;
157  // we could copy the diagonal element twice or check
158  // doing the check may be a tiny bit faster,
159  // so we choose that option for now
160  if(k!=j) m[k*nrow+j] = *sjk;
161  ++sjk;
162  }
163  }
164 }
165 
167  : m(hm1.nrow*hm1.nrow), nrow(hm1.nrow), ncol(hm1.nrow)
168 {
169  size_ = nrow * ncol;
170 
171  int n = num_row();
172  mIter mrr;
173  mcIter mr = hm1.m.begin();
174  for(int r=0;r<n;r++) {
175  mrr = m.begin()+(n+1)*r;
176  *mrr = *(mr++);
177  }
178 }
179 
181  : m(hm1.nrow), nrow(hm1.nrow), ncol(1)
182 {
183 
184  size_ = nrow;
185  m = hm1.m;
186 }
187 
188 
189 //
190 //
191 // Sub matrix
192 //
193 //
194 
195 HepMatrix HepMatrix::sub(int min_row, int max_row,
196  int min_col,int max_col) const
197 #ifdef HEP_GNU_OPTIMIZED_RETURN
198 return mret(max_row-min_row+1,max_col-min_col+1);
199 {
200 #else
201 {
202  HepMatrix mret(max_row-min_row+1,max_col-min_col+1);
203 #endif
204  if(max_row > num_row() || max_col >num_col())
205  error("HepMatrix::sub: Index out of range");
206  mIter a = mret.m.begin();
207  int nc = num_col();
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) {
211  mcIter brc = b1;
212  for(int icol=0; icol<mret.num_col(); ++icol) {
213  *(a++) = *(brc++);
214  }
215  if(irow<rowsize) b1 += nc;
216  }
217  return mret;
218 }
219 
220 void HepMatrix::sub(int row,int col,const HepMatrix &hm1)
221 {
222  if(row <1 || row+hm1.num_row()-1 > num_row() ||
223  col <1 || col+hm1.num_col()-1 > num_col() )
224  error("HepMatrix::sub: Index out of range");
225  mcIter a = hm1.m.begin();
226  int nc = num_col();
227  mIter b1 = m.begin() + (row - 1) * nc + col - 1;
228  int rowsize = hm1.num_row();
229  for(int irow=1; irow<=rowsize; ++irow) {
230  mIter brc = b1;
231  for(int icol=0; icol<hm1.num_col(); ++icol) {
232  *(brc++) = *(a++);
233  }
234  if(irow<rowsize) b1 += nc;
235  }
236 }
237 
238 //
239 // Direct sum of two matricies
240 //
241 
242 HepMatrix dsum(const HepMatrix &hm1, const HepMatrix &hm2)
243 #ifdef HEP_GNU_OPTIMIZED_RETURN
244  return mret(hm1.num_row() + hm2.num_row(), hm1.num_col() + hm2.num_col(),
245  0);
246 {
247 #else
248 {
249  HepMatrix mret(hm1.num_row() + hm2.num_row(), hm1.num_col() + hm2.num_col(),
250  0);
251 #endif
252  mret.sub(1,1,hm1);
253  mret.sub(hm1.num_row()+1,hm1.num_col()+1,hm2);
254  return mret;
255 }
256 
257 /* -----------------------------------------------------------------------
258  This section contains support routines for matrix.h. This section contains
259  The two argument functions +,-. They call the copy constructor and +=,-=.
260  ----------------------------------------------------------------------- */
262 #ifdef HEP_GNU_OPTIMIZED_RETURN
263  return hm2(nrow, ncol);
264 {
265 #else
266 {
267  HepMatrix hm2(nrow, ncol);
268 #endif
269  mcIter a=m.begin();
270  mIter b=hm2.m.begin();
271  mcIter e=m.end();
272  for(;a<e; a++, b++) (*b) = -(*a);
273  return hm2;
274 }
275 
276 
277 
278 HepMatrix operator+(const HepMatrix &hm1,const HepMatrix &hm2)
279 #ifdef HEP_GNU_OPTIMIZED_RETURN
280  return mret(hm1.nrow, hm1.ncol);
281 {
282 #else
283 {
284  HepMatrix mret(hm1.nrow, hm1.ncol);
285 #endif
286  CHK_DIM_2(hm1.num_row(),hm2.num_row(), hm1.num_col(),hm2.num_col(),+);
287  SIMPLE_TOP(+)
288  return mret;
289 }
290 
291 //
292 // operator -
293 //
294 
295 HepMatrix operator-(const HepMatrix &hm1,const HepMatrix &hm2)
296 #ifdef HEP_GNU_OPTIMIZED_RETURN
297  return mret(hm1.num_row(), hm1.num_col());
298 {
299 #else
300 {
301  HepMatrix mret(hm1.num_row(), hm1.num_col());
302 #endif
303  CHK_DIM_2(hm1.num_row(),hm2.num_row(),
304  hm1.num_col(),hm2.num_col(),-);
305  SIMPLE_TOP(-)
306  return mret;
307 }
308 
309 /* -----------------------------------------------------------------------
310  This section contains support routines for matrix.h. This file contains
311  The two argument functions *,/. They call copy constructor and then /=,*=.
312  ----------------------------------------------------------------------- */
313 
315 const HepMatrix &hm1,double t)
316 #ifdef HEP_GNU_OPTIMIZED_RETURN
317  return mret(hm1);
318 {
319 #else
320 {
321  HepMatrix mret(hm1);
322 #endif
323  mret /= t;
324  return mret;
325 }
326 
327 HepMatrix operator*(const HepMatrix &hm1,double t)
328 #ifdef HEP_GNU_OPTIMIZED_RETURN
329  return mret(hm1);
330 {
331 #else
332 {
333  HepMatrix mret(hm1);
334 #endif
335  mret *= t;
336  return mret;
337 }
338 
339 HepMatrix operator*(double t,const HepMatrix &hm1)
340 #ifdef HEP_GNU_OPTIMIZED_RETURN
341  return mret(hm1);
342 {
343 #else
344 {
345  HepMatrix mret(hm1);
346 #endif
347  mret *= t;
348  return mret;
349 }
350 
351 HepMatrix operator*(const HepMatrix &hm1,const HepMatrix &hm2)
352 #ifdef HEP_GNU_OPTIMIZED_RETURN
353  return mret(hm1.nrow,hm2.ncol,0);
354 {
355 #else
356 {
357  // initialize matrix to 0.0
358  HepMatrix mret(hm1.nrow,hm2.ncol,0);
359 #endif
360  CHK_DIM_1(hm1.ncol,hm2.nrow,*);
361 
362  int m1cols = hm1.ncol;
363  int m2cols = hm2.ncol;
364 
365  for (int i=0; i<hm1.nrow; i++)
366  {
367  for (int j=0; j<m1cols; j++)
368  {
369  register double temp = hm1.m[i*m1cols+j];
370  HepMatrix::mIter pt = mret.m.begin() + i*m2cols;
371 
372  // Loop over k (the column index in matrix hm2)
373  HepMatrix::mcIter pb = hm2.m.begin() + m2cols*j;
374  const HepMatrix::mcIter pblast = pb + m2cols;
375  while (pb < pblast)
376  {
377  (*pt) += temp * (*pb);
378  pb++;
379  pt++;
380  }
381  }
382  }
383 
384  return mret;
385 }
386 
387 /* -----------------------------------------------------------------------
388  This section contains the assignment and inplace operators =,+=,-=,*=,/=.
389  ----------------------------------------------------------------------- */
390 
392 {
393  CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),+=);
394  SIMPLE_BOP(+=)
395  return (*this);
396 }
397 
399 {
400  CHK_DIM_2(num_row(),hm2.num_row(),num_col(),hm2.num_col(),-=);
401  SIMPLE_BOP(-=)
402  return (*this);
403 }
404 
406 {
407  SIMPLE_UOP(/=)
408  return (*this);
409 }
410 
412 {
413  SIMPLE_UOP(*=)
414  return (*this);
415 }
416 
418 {
419  if(hm1.nrow*hm1.ncol != size_) //??fixme?? hm1.size != size
420  {
421  size_ = hm1.nrow * hm1.ncol;
422  m.resize(size_); //??fixme?? if (size < hm1.size) m.resize(hm1.size);
423  }
424  nrow = hm1.nrow;
425  ncol = hm1.ncol;
426  m = hm1.m;
427  return (*this);
428 }
429 
430 // HepMatrix & HepMatrix::operator=(const HepRotation &hm2)
431 // is now in Matrix=Rotation.cc
432 
433 // Print the Matrix.
434 
435 std::ostream& operator<<(std::ostream &os, const HepMatrix &q)
436 {
437  os << "\n";
438 /* Fixed format needs 3 extra characters for field, while scientific needs 7 */
439  int width;
440  if(os.flags() & std::ios::fixed)
441  width = os.precision()+3;
442  else
443  width = os.precision()+7;
444  for(int irow = 1; irow<= q.num_row(); irow++)
445  {
446  for(int icol = 1; icol <= q.num_col(); icol++)
447  {
448  os.width(width);
449  os << q(irow,icol) << " ";
450  }
451  os << std::endl;
452  }
453  return os;
454 }
455 
457 #ifdef HEP_GNU_OPTIMIZED_RETURN
458 return mret(ncol,nrow);
459 {
460 #else
461 {
462  HepMatrix mret(ncol,nrow);
463 #endif
464  mcIter pme = m.begin();
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;
469  (*pt) = (*pme);
470  ++pme;
471  }
472  }
473  return mret;
474 }
475 
476 HepMatrix HepMatrix::apply(double (*f)(double, int, int)) const
477 #ifdef HEP_GNU_OPTIMIZED_RETURN
478 return mret(num_row(),num_col());
479 {
480 #else
481 {
482  HepMatrix mret(num_row(),num_col());
483 #endif
484  mcIter a = m.begin();
485  mIter b = mret.m.begin();
486  for(int ir=1;ir<=num_row();ir++) {
487  for(int ic=1;ic<=num_col();ic++) {
488  *(b++) = (*f)(*(a++), ir, ic);
489  }
490  }
491  return mret;
492 }
493 
494 int HepMatrix::dfinv_matrix(int *ir) {
495  if (num_col()!=num_row())
496  error("dfinv_matrix: Matrix is not NxN");
497  int n = num_col();
498  if (n==1) return 0;
499 
500  double s31, s32;
501  register double s33, s34;
502 
503  mIter hm11 = m.begin();
504  mIter hm12 = hm11 + 1;
505  mIter hm21 = hm11 + n;
506  mIter hm22 = hm12 + n;
507  *hm21 = -(*hm22) * (*hm11) * (*hm21);
508  *hm12 = -(*hm12);
509  if (n>2) {
510  mIter mimim = hm11 + n + 1;
511  for (int i=3;i<=n;i++) {
512  // calculate these to avoid pointing off the end of the storage array
513  mIter mi = hm11 + (i-1) * n;
514  mIter mii= hm11 + (i-1) * n + i - 1;
515  int ihm2 = i - 2;
516  mIter mj = hm11;
517  mIter mji = mj + i - 1;
518  mIter mij = mi;
519  for (int j=1;j<=ihm2;j++) {
520  s31 = 0.0;
521  s32 = *mji;
522  mIter mkj = mj + j - 1;
523  mIter mik = mi + j - 1;
524  mIter mjkp = mj + j;
525  mIter mkpi = mj + n + i - 1;
526  for (int k=j;k<=ihm2;k++) {
527  s31 += (*mkj) * (*(mik++));
528  s32 += (*(mjkp++)) * (*mkpi);
529  mkj += n;
530  mkpi += n;
531  } // for k
532  *mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31));
533  *mji = -s32;
534  mj += n;
535  mji += n;
536  mij++;
537  } // for j
538  *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
539  *(mimim+1) = -(*(mimim+1));
540  mimim += (n+1);
541  } // for i
542  } // n>2
543  mIter mi = hm11;
544  mIter mii = hm11;
545  for (int i=1;i<n;i++) {
546  int ni = n - i;
547  mIter mij = mi;
548  int j;
549  for (j=1; j<=i;j++) {
550  s33 = *mij;
551  // change initial definition of mikj to avoid pointing off the end of the storage array
552  mIter mikj = mi + j - 1;
553  mIter miik = mii + 1;
554  mIter min_end = mi + n;
555  for (;miik<min_end;) {
556  // iterate by n as we enter the loop to avoid pointing off the end of the storage array
557  mikj += n;
558  s33 += (*mikj) * (*(miik++));
559  }
560  *(mij++) = s33;
561  }
562  for (j=1;j<=ni;j++) {
563  s34 = 0.0;
564  mIter miik = mii + j;
565  for (int k=j;k<=ni;k++) {
566  // calculate mikij here to avoid pointing off the end of the storage array
567  mIter mikij = mii + k * n + j;
568  s34 += *mikij * (*(miik++));
569  }
570  *(mii+j) = s34;
571  }
572  mi += n;
573  mii += (n+1);
574  } // for i
575  int nxch = ir[n];
576  if (nxch==0) return 0;
577  for (int hmm=1;hmm<=nxch;hmm++) {
578  int k = nxch - hmm + 1;
579  int ij = ir[k];
580  int i = ij >> 12;
581  int j = ij%4096;
582  for (k=1; k<=n;k++) {
583  // avoid setting the iterator beyond the end of the storage vector
584  mIter mki = hm11 + (k-1)*n + i - 1;
585  mIter mkj = hm11 + (k-1)*n + j - 1;
586  // 2/24/05 David Sachs fix of improper swap bug that was present
587  // for many years:
588  double ti = *mki; // 2/24/05
589  *mki = *mkj;
590  *mkj = ti; // 2/24/05
591  }
592  } // for hmm
593  return 0;
594 }
595 
596 int HepMatrix::dfact_matrix(double &det, int *ir) {
597  if (ncol!=nrow)
598  error("dfact_matrix: Matrix is not NxN");
599 
600  int ifail, jfail;
601  int n = ncol;
602 
603  double tf;
604  double g1 = 1.0e-19, g2 = 1.0e19;
605 
606  double p, q, t;
607  double s11, s12;
608 
609  double epsilon = 8*DBL_EPSILON;
610  // could be set to zero (like it was before)
611  // but then the algorithm often doesn't detect
612  // that a matrix is singular
613 
614  int normal = 0, imposs = -1;
615  int jrange = 0, jover = 1, junder = -1;
616  ifail = normal;
617  jfail = jrange;
618  int nxch = 0;
619  det = 1.0;
620  mIter mj = m.begin();
621  mIter mjj = mj;
622  for (int j=1;j<=n;j++) {
623  int k = j;
624  p = (fabs(*mjj));
625  if (j!=n) {
626  // replace mij with calculation of position
627  for (int i=j+1;i<n;i++) {
628  q = (fabs(*(mj + n*(i-j) + j - 1)));
629  if (q > p) {
630  k = i;
631  p = q;
632  }
633  } // for i
634  if (k==j) {
635  if (p <= epsilon) {
636  det = 0;
637  ifail = imposs;
638  jfail = jrange;
639  return ifail;
640  }
641  det = -det; // in this case the sign of the determinant
642  // must not change. So I change it twice.
643  } // k==j
644  mIter mjl = mj;
645  mIter mkl = m.begin() + (k-1)*n;
646  for (int l=1;l<=n;l++) {
647  tf = *mjl;
648  *(mjl++) = *mkl;
649  *(mkl++) = tf;
650  }
651  nxch = nxch + 1; // this makes the determinant change its sign
652  ir[nxch] = (((j)<<12)+(k));
653  } else { // j!=n
654  if (p <= epsilon) {
655  det = 0.0;
656  ifail = imposs;
657  jfail = jrange;
658  return ifail;
659  }
660  } // j!=n
661  det *= *mjj;
662  *mjj = 1.0 / *mjj;
663  t = (fabs(det));
664  if (t < g1) {
665  det = 0.0;
666  if (jfail == jrange) jfail = junder;
667  } else if (t > g2) {
668  det = 1.0;
669  if (jfail==jrange) jfail = jover;
670  }
671  // calculate mk and mkjp so we don't point off the end of the vector
672  if (j!=n) {
673  mIter mjk = mj + j;
674  for (k=j+1;k<=n;k++) {
675  mIter mk = mj + n*(k-j);
676  mIter mkjp = mk + j;
677  s11 = - (*mjk);
678  s12 = - (*mkjp);
679  if (j!=1) {
680  mIter mik = m.begin() + k - 1;
681  mIter mijp = m.begin() + j;
682  mIter mki = mk;
683  mIter mji = mj;
684  for (int i=1;i<j;i++) {
685  s11 += (*mik) * (*(mji++));
686  s12 += (*mijp) * (*(mki++));
687  mik += n;
688  mijp += n;
689  } // for i
690  } // j!=1
691  *(mjk++) = -s11 * (*mjj);
692  *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
693  } // for k
694  } // j!=n
695  // avoid setting the iterator beyond the end of the vector
696  if(j!=n) {
697  mj += n;
698  mjj += (n+1);
699  }
700  } // for j
701  if (nxch%2==1) det = -det;
702  if (jfail !=jrange) det = 0.0;
703  ir[n] = nxch;
704  return 0;
705 }
706 
707 void HepMatrix::invert(int &ierr) {
708  if(ncol != nrow)
709  error("HepMatrix::invert: Matrix is not NxN");
710 
711  static int max_array = 20;
712  static int *ir = new int [max_array+1];
713 
714  if (ncol > max_array) {
715  delete [] ir;
716  max_array = nrow;
717  ir = new int [max_array+1];
718  }
719  double t1, t2, t3;
720  double det, temp, sd;
721  int ifail;
722  switch(nrow) {
723  case 3:
724  double c11,c12,c13,c21,c22,c23,c31,c32,c33;
725  ifail = 0;
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));
738  if (t1 >= t2) {
739  if (t3 >= t1) {
740  temp = *(m.begin()+6);
741  det = c23*c12-c22*c13;
742  } else {
743  temp = *(m.begin());
744  det = c22*c33-c23*c32;
745  }
746  } else if (t3 >= t2) {
747  temp = *(m.begin()+6);
748  det = c23*c12-c22*c13;
749  } else {
750  temp = *(m.begin()+3);
751  det = c13*c32-c12*c33;
752  }
753  if (det==0) {
754  ierr = 1;
755  return;
756  }
757  {
758  double s1 = temp/det;
759  mIter hmm = m.begin();
760  *(hmm++) = s1*c11;
761  *(hmm++) = s1*c21;
762  *(hmm++) = s1*c31;
763  *(hmm++) = s1*c12;
764  *(hmm++) = s1*c22;
765  *(hmm++) = s1*c32;
766  *(hmm++) = s1*c13;
767  *(hmm++) = s1*c23;
768  *(hmm) = s1*c33;
769  }
770  break;
771  case 2:
772  ifail = 0;
773  det = (*m.begin())*(*(m.begin()+3)) - (*(m.begin()+1))*(*(m.begin()+2));
774  if (det==0) {
775  ierr = 1;
776  return;
777  }
778  sd = 1.0/det;
779  temp = sd*(*(m.begin()+3));
780  *(m.begin()+1) *= -sd;
781  *(m.begin()+2) *= -sd;
782  *(m.begin()+3) = sd*(*m.begin());
783  *(m.begin()) = temp;
784  break;
785  case 1:
786  ifail = 0;
787  if ((*(m.begin()))==0) {
788  ierr = 1;
789  return;
790  }
791  *(m.begin()) = 1.0/(*(m.begin()));
792  break;
793  case 4:
794  invertHaywood4(ierr);
795  return;
796  case 5:
797  invertHaywood5(ierr);
798  return;
799  case 6:
800  invertHaywood6(ierr);
801  return;
802  default:
803  ifail = dfact_matrix(det, ir);
804  if(ifail) {
805  ierr = 1;
806  return;
807  }
808  dfinv_matrix(ir);
809  break;
810  }
811  ierr = 0;
812  return;
813 }
814 
815 double HepMatrix::determinant() const {
816  static int max_array = 20;
817  static int *ir = new int [max_array+1];
818  if(ncol != nrow)
819  error("HepMatrix::determinant: Matrix is not NxN");
820  if (ncol > max_array) {
821  delete [] ir;
822  max_array = nrow;
823  ir = new int [max_array+1];
824  }
825  double det;
826  HepMatrix mt(*this);
827  int i = mt.dfact_matrix(det, ir);
828  if(i==0) return det;
829  return 0;
830 }
831 
832 double HepMatrix::trace() const {
833  double t = 0.0;
834  for (mcIter d = m.begin(); d < m.end(); d += (ncol+1) )
835  t += *d;
836  return t;
837 }
838 
839 } // namespace CLHEP
virtual void invertHaywood5(int &ierr)
HepMatrix & operator=(const HepMatrix &)
Definition: Matrix.cc:417
double trace() const
Definition: Matrix.cc:832
HepLorentzVector operator/(const HepLorentzVector &, double a)
#define CHK_DIM_1(c1, r2, fun)
Definition: Matrix.cc:56
HepMatrix apply(double(*f)(double, int, int)) const
Definition: Matrix.cc:476
HepMatrix operator-() const
Definition: Matrix.cc:261
friend HepMatrix operator+(const HepMatrix &hm1, const HepMatrix &hm2)
Definition: Matrix.cc:278
HepMatrix sub(int min_row, int max_row, int min_col, int max_col) const
Definition: Matrix.cc:195
std::vector< double, Alloc< double, 25 > >::const_iterator mcIter
friend HepMatrix operator*(const HepMatrix &hm1, const HepMatrix &hm2)
Definition: Matrix.cc:351
virtual void invertHaywood6(int &ierr)
virtual ~HepMatrix()
Definition: Matrix.cc:108
#define SIMPLE_UOP(OPER)
Definition: Matrix.cc:31
HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2)
Definition: DiagMatrix.cc:164
virtual int num_col() const
Definition: Matrix.cc:122
#define SIMPLE_BOP(OPER)
Definition: Matrix.cc:36
HepMatrix & operator/=(double t)
Definition: Matrix.cc:405
void init()
Definition: testRandom.cc:27
std::vector< double, Alloc< double, 25 > >::iterator mIter
HepMatrix T() const
Definition: Matrix.cc:456
HepMatrix & operator*=(double t)
Definition: Matrix.cc:411
virtual int num_size() const
Definition: Matrix.cc:124
void f(void g())
Definition: excDblThrow.cc:38
virtual const double & operator()(int row, int col) const
Definition: Matrix.cc:137
static void error(const char *s)
Definition: GenMatrix.cc:73
virtual void invertHaywood4(int &ierr)
virtual int num_row() const
Definition: Matrix.cc:120
HepMatrix & operator-=(const HepMatrix &)
Definition: Matrix.cc:398
double determinant() const
Definition: Matrix.cc:815
HepMatrix & operator+=(const HepMatrix &)
Definition: Matrix.cc:391
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition: AxisAngle.cc:86
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition: Matrix.cc:51
#define SIMPLE_TOP(OPER)
Definition: Matrix.cc:42