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
27namespace 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
69HepMatrix::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
120int HepMatrix::num_row() const { return nrow;}
121
122int HepMatrix::num_col() const { return ncol;}
123
124int HepMatrix::num_size() const { return size_;}
125
126// operator()
127
128double & 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
137const 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
195HepMatrix HepMatrix::sub(int min_row, int max_row,
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);
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
220void 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
242HepMatrix 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
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
295HepMatrix 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
314HepMatrix operator/(
315const 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
327HepMatrix 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
339HepMatrix 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
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
435std::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
458return 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
476HepMatrix HepMatrix::apply(double (*f)(double, int, int)) const
477#ifdef HEP_GNU_OPTIMIZED_RETURN
478return 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
494int 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
596int 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
707void 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
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
832double 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
#define CHK_DIM_2(r1, r2, c1, c2, fun)
Definition: Matrix.cc:51
#define SIMPLE_BOP(OPER)
Definition: Matrix.cc:36
#define SIMPLE_UOP(OPER)
Definition: Matrix.cc:31
#define SIMPLE_TOP(OPER)
Definition: Matrix.cc:42
#define CHK_DIM_1(c1, r2, fun)
Definition: Matrix.cc:56
std::vector< double, Alloc< double, 25 > >::const_iterator mcIter
std::vector< double, Alloc< double, 25 > >::iterator mIter
static void error(const char *s)
Definition: GenMatrix.cc:73
virtual ~HepMatrix()
Definition: Matrix.cc:108
virtual void invertHaywood4(int &ierr)
HepMatrix sub(int min_row, int max_row, int min_col, int max_col) const
Definition: Matrix.cc:195
HepMatrix & operator=(const HepMatrix &)
Definition: Matrix.cc:417
HepMatrix apply(double(*f)(double, int, int)) const
Definition: Matrix.cc:476
virtual int num_col() const
Definition: Matrix.cc:122
virtual const double & operator()(int row, int col) const
Definition: Matrix.cc:137
HepMatrix & operator/=(double t)
Definition: Matrix.cc:405
virtual int num_row() const
Definition: Matrix.cc:120
HepMatrix & operator+=(const HepMatrix &)
Definition: Matrix.cc:391
HepMatrix operator-() const
Definition: Matrix.cc:261
HepMatrix & operator*=(double t)
Definition: Matrix.cc:411
HepMatrix & operator-=(const HepMatrix &)
Definition: Matrix.cc:398
HepMatrix T() const
Definition: Matrix.cc:456
double determinant() const
Definition: Matrix.cc:815
virtual void invertHaywood6(int &ierr)
virtual int num_size() const
Definition: Matrix.cc:124
virtual void invertHaywood5(int &ierr)
double trace() const
Definition: Matrix.cc:832
void f(void g())
Definition: excDblThrow.cc:38
Hep3Vector operator+(const Hep3Vector &, const Hep3Vector &)
HepLorentzRotation operator*(const HepRotation &r, const HepLorentzRotation &lt)
std::ostream & operator<<(std::ostream &os, const HepAxisAngle &aa)
Definition: AxisAngle.cc:86
HepDiagMatrix dsum(const HepDiagMatrix &s1, const HepDiagMatrix &s2)
Definition: DiagMatrix.cc:164
@ b
@ a
void init()
Definition: testRandom.cc:27