StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StMatrix.hh
1 /***************************************************************************
2  *
3  * $Id: StMatrix.hh,v 1.21 2016/01/22 17:10:50 smirnovd Exp $
4  *
5  * Author: Original code from CLHEP by Mike Smyth
6  * Modified April 17, 1998 Brian Lasiuk (templated version)
7  *
8  ***************************************************************************
9  *
10  * Description:
11  * StMatrix class, does basic matrix operations
12  *
13  * Remarks: Since not all compilers support member templates
14  * we have to specialize the templated member on these
15  * platforms. If member templates are not supported the
16  * ST_NO_MEMBER_TEMPLATES flag has to be set. tu.
17  *
18  ***************************************************************************
19  *
20  * $Log: StMatrix.hh,v $
21  * Revision 1.21 2016/01/22 17:10:50 smirnovd
22  * StarClassLibrary: Removed deprecated storage class specifier 'register'
23  *
24  * This keyword is deprecated since C++11 and serves no purpose
25  *
26  * "
27  * The register specifier is only allowed for objects declared at block scope and
28  * in function parameter lists. It indicates automatic storage duration, which is
29  * the default for these kinds of declarations. Additionally, the presence of this
30  * keyword may be used as a hint for the optimizer to store the value of this
31  * variable in a CPU register.
32  * "
33  *
34  * Revision 1.20 2012/06/11 15:29:26 fisyak
35  * std namespace
36  *
37  * Revision 1.19 2006/04/07 22:02:34 ullrich
38  * Fixed bug in dfinv and dfact that were affecting
39  * determinant(), invert(), and inverse().
40  *
41  * Revision 1.18 2006/01/09 23:47:27 fisyak
42  * Add missing methods (found by Zhangbu) to Cint dictionary
43  *
44  * Revision 1.17 2005/09/22 20:09:20 fisyak
45  * Make StLorentzVector persistent
46  *
47  * Revision 1.16 2005/07/06 18:49:56 fisyak
48  * Replace StHelixD, StLorentzVectorD,StLorentzVectorF,StMatrixD,StMatrixF,StPhysicalHelixD,StThreeVectorD,StThreeVectorF by templated version
49  *
50  * Revision 1.15 2004/01/14 22:37:27 fisyak
51  * unsigned int => size_t to make alpha happy
52  *
53  * Revision 1.14 2003/09/02 17:59:35 perev
54  * gcc 3.2 updates + WarnOff
55  *
56  * Revision 1.13 2002/07/02 22:20:53 ullrich
57  * Add dummy return statements to get rid of warnings under Linux.
58  *
59  * Revision 1.12 2001/10/31 15:11:36 ullrich
60  * Rewrote swap() to work as non-friend to StMatrix.
61  *
62  * Revision 1.11 2001/10/31 00:33:34 ullrich
63  * Remove macro ifdef for GCC which is not needed anymore.
64  *
65  * Revision 1.10 2000/02/02 18:31:05 lasiuk
66  * restore files
67  *
68  * Revision 1.9 2000/02/01 16:02:59 lasiuk
69  * namespace std is different on SUN CC5 and KCC. Redefine macros!
70  *
71  * Revision 1.8 2000/01/31 20:53:45 lasiuk
72  * using std::swap
73  *
74  * Revision 1.7 1999/12/21 15:14:13 ullrich
75  * Modified to cope with new compiler version on Sun (CC5.0).
76  *
77  * Revision 1.6 1999/12/07 23:43:04 ullrich
78  * Modified to get rid of warnings on Linux.
79  *
80  * Revision 1.5 1999/04/14 23:12:07 fisyak
81  * Add __CINT__ to handle references
82  *
83  * Revision 1.4 1999/03/04 18:12:24 ullrich
84  * Added namespace 'std'.
85  *
86  * Revision 1.3 1999/02/17 11:38:54 ullrich
87  * Removed specialization for 'long double'.
88  *
89  * Revision 1.2 1999/02/14 23:11:43 fisyak
90  * Fixes for Rootcint
91  *
92  * Revision 1.1 1999/01/30 03:59:03 fisyak
93  * Root Version of StarClassLibrary
94  *
95  * Revision 1.1 1999/01/23 00:27:56 ullrich
96  * Initial Revision
97  *
98  **************************************************************************/
99 // This code was modified from the CLHEP files:
100 // GenMatrixD.h/.icc/.cc MatrixD.h/.icc/.cc
101 // Comments that were kept from the original code are preceding
102 // by ///.
103 
123 // The Matrix class provides the tools to do matrix arithmetic.
124 // The Matrices are templates with the default type being "double"
125 // A Matrix is declared by:
126 //
127 // StMatrix<double> m1(n,m); // n x m matrix
128 // StMatrix<> m2(n,m); // default type is double
129 //
130 // A Matrix may be initialized by giving it a third integer
131 // argument, either 0 or 1:
132 // - Zero (0) means initialize all elements to 0.
133 // - One (1) means initialize to the identity matrix if the matrix is
134 // square (NxN); otherwise an exception is thrown.
135 // - The default behavior is to initialize all elements to zero (0)
136 //
137 // StMatrix<double> m2(n, m, 0); // all zeros
138 // StMatrix<double> m1(n, n, 1); // diagnol elements are 1
139 // StMatrix<double> m1(n, m, 1); // error if m!=n
140 //
141 // Copying and Assignment:
142 //
143 // StMatrix<double> m1(m2); // m2 copied to m1
144 // StMatrix<double> m1 = m2; // m1 assigned value of m1
145 //
146 // Element Access
147 // Single elements of the matrix can be accessed in two differernt ways:
148 // - C array style operators []
149 //
150 // For a pXq matrix:
151 // --> m1[i][j]:
152 // 0 <= i <= (p-1)
153 // 0 <= j <= (q-1)
154 //
155 // - Fortran array style operators ()
156 // --> m1(i,j):
157 // 1 <= i <= p
158 // 1 <= j <= q
159 //
166 //
168 // #define MATRIX_BOUND_CHECK
169 // before including StMatrix.hh
171 // Access Functions:
172 // Properties of the matrices can be extracted via the following
173 // access functions:
174 //
175 // m1.num_row(); // number of rows
176 // m1.numRow(); // "
177 // m1.num_col(); // number of columns
178 // m1.numCol(); // "
179 //
181 // - Scaler operations:
182 // A matrix may be multiplied or divided by a scaler
183 // via the /= and *= operators:
184 // m1 /= 5; // divides all elements in matrix m1 by 5
185 //
186 // - Matrix operations:
187 // Two matrices can be added, subtracted, and multiplied with the
188 // operators +=, -=, *=
189 //
190 // m1 += m2; // same as m1 = m1+m2;
191 //
196 // mt = m.transpose() ; // "
197 //
198 // unsigned int ierr; // unsigned int == unsigned int
200 // // ierr will be non-zero if the matrix is singular.
202 // // faster (avoids copy operator)
205 // // m2(row_min:row_max, col_min:col_max)
207 // // is replaced with m1.
210 // // applies f to every element of m1 and places it in m.
211 //
214 //
215 // There also exists functionality to multiply a
216 // 1) StThreeVector
217 // 2) StLorentzVector
218 // 3) Std C++ Library vector<>
219 // by a matrix.
220 //
221 // StThreeVector<> a(1,1,1);
222 // StMatrix<> M1(3,3,1);
223 // StThreeVector<> b = a*M1; // multiplication by identity matrix
224 //
225 // This can be used for simple rotations. The same syntax is also used
226 // for the LorentzVector and Std C++ vector. In order to use the
227 // StThreeVector and StLorentzVector capablilties simply define:
228 //
229 // #define WITH_ST_THREEVECTOR // for use with StThreeVector
230 // #define WITH_ST_LORENTZVECTOR // for use with StLorentzVector
239 
240 #ifndef ST_MATRIX_HH
241 #define ST_MATRIX_HH
242 
243 #include <Stiostream.h>
244 #include <math.h>
245 #include <vector>
246 #include <float.h>
247 #if !defined(ST_NO_NAMESPACES)
248 using std::vector;
249 #endif
250 
251 #ifndef ST_NO_EXCEPTIONS
252 # include <stdexcept>
253 # if !defined(ST_NO_NAMESPACES)
254  using std::out_of_range;
255  using std::domain_error;
256 # endif
257 #endif
258 
259 #include "StThreeVector.hh"
260 #include "StLorentzVector.hh"
261 
262 // in types.h assumed:
263 // typedef unsigned int size_t;
264 
265 template<class DataType> class StMatrix {
266 public:
267  StMatrix();
268  StMatrix(size_t p, size_t q, size_t init=0);
269 
270  // Copy constructor.
271 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
272  template<class X>
273  StMatrix(const StMatrix<X>&);
274  StMatrix(const StMatrix<DataType>&);
275 #else
276  StMatrix(const StMatrix<float>&);
277  StMatrix(const StMatrix<double>&);
278 #endif
279 
280  // Assignment operators.
281 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
282  template<class X>
283  StMatrix<DataType>& operator=(const StMatrix<X>&);
284  StMatrix<DataType>& operator=(const StMatrix<DataType>&);
285 #else
286  StMatrix<DataType>& operator=(const StMatrix<float>&);
287  StMatrix<DataType>& operator=(const StMatrix<double>&);
288 #endif
289 
290  // Destructor. --
291  //Problem with LINUX (virtual table error when exception thrown)
292  virtual ~StMatrix();
293 
294  //
295  // access functions
296  //
297  unsigned int numRow() const;
298  unsigned int numCol() const;
299  unsigned int numSize() const;
300  unsigned int num_row() const;
301  unsigned int num_col() const;
302  unsigned int num_size() const;
303  //
304  // Element access:
305  // ** Note that the indexing starts from:
306  // For (,) --> (1,1). **
307  // For [][] --> [0][0]
308  const DataType& operator()(size_t row, size_t col) const;
309  DataType& operator()(size_t row, size_t col);
310 
311  // classes for implementing m[i][j]
312  class StMatrixRow {
313  public:
314  StMatrixRow(StMatrix<DataType>&, size_t);
315  DataType& operator[](size_t);
316  private:
317 #ifndef __CINT__
318  StMatrix<DataType>& _a;
319 #else
320  StMatrix<DataType>* _a;
321 #endif
322  size_t _r;
323  };
324  class StMatrixRowConst {
325  public:
326  StMatrixRowConst (const StMatrix<DataType>&, size_t);
327  const DataType & operator[](size_t) const;
328  private:
329 #ifndef __CINT__
330  const StMatrix<DataType>& _a;
331 #else
332  const StMatrix<DataType>* _a;
333 #endif
334  size_t _r;
335  };
336 
337  StMatrixRow operator[] (size_t r);
338  const StMatrixRowConst operator[] (size_t) const;
339 
340  // Scaler Operations
341  StMatrix<DataType>& operator*=(double t);
342  StMatrix<DataType>& operator/=(double t);
343 
344  // Matrix Operations
345 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
346  template<class X> StMatrix<DataType>& operator+=(const StMatrix<X>&);
347  template<class X> StMatrix<DataType>& operator-=(const StMatrix<X>&);
348  template<class X> StMatrix<DataType> dot(const StMatrix<X>&);
349 #else
350  StMatrix<DataType>& operator+=(const StMatrix<float>&);
351  StMatrix<DataType>& operator+=(const StMatrix<double>&);
352 
353  StMatrix<DataType>& operator-=(const StMatrix<float>&);
354  StMatrix<DataType>& operator-=(const StMatrix<double>&);
355 
356  StMatrix<DataType> dot(const StMatrix<float>&);
357  StMatrix<DataType> dot(const StMatrix<double>&);
358 #endif
359 
360  // unary Operations
361  StMatrix<DataType> operator+ () const;
362  StMatrix<DataType> operator- () const;
363 
364  // Equality
365  bool operator==(const StMatrix<DataType>&) const;
366  bool operator!=(const StMatrix<DataType>&) const;
367 
368  // Apply a function to all elements of the matrix.
369  StMatrix<DataType> apply(DataType (*f)(DataType, size_t, size_t)) const;
370 
371  StMatrix<DataType> T() const;
372  StMatrix<DataType> transpose() const;
373 
374  StMatrix<DataType> sub(size_t min_row, size_t max_row, size_t min_col, size_t max_col) const;
375 
376  void sub(size_t row, size_t col, const StMatrix<DataType> &m1);
377 
378  StMatrix<DataType> inverse(size_t& ierr) const;
379  void invert(size_t& ierr);
380 
381  DataType determinant() const;
382 
383  // Must be of same type to swap
384  static void swap(unsigned int&, unsigned int&);
385  static void swap(DataType *&, DataType *&);
386 
387 private:
388  // Friend classes.
389  friend class StMatrixRow;
390  friend class StMatrixRowConst;
391 
392  // If successful, the return code is 0.
393  // On return, det is the determinant and ir[] is row-unsigned interchange
394  // matrix. See CERNLIB's DFACT routine.
395  int dfact(DataType &det, int *ir); // factorize the matrix.
396 
397  // invert the matrix. See CERNLIB DFINV.
398  int dfinv(int *ir);
399 
400 protected:
401  DataType *mElement;
402  unsigned int mRow, mCol;
403  unsigned int mSize;
404 #ifdef __ROOT__
405  ClassDef(StMatrix,3)
406 #endif /* __ROOT__ */
407 };
408 
409 #ifndef __CINT__
410 // Constructors.
411 
412 template<class DataType>
413 inline StMatrix<DataType>::StMatrix()
414  : mElement(0), mRow(0), mCol(0), mSize(0) {}
415 
416 template<class DataType>
417 StMatrix<DataType>::StMatrix(size_t p,size_t q, size_t init)
418  : mRow(p), mCol(q)
419 {
420  mSize = mRow*mCol;
421  mElement = new DataType[mSize];
422 
423  DataType *a = mElement;
424  DataType *b = mElement + mSize;
425  for( ; a<b; a++)
426  *a = 0;
427 
428  if (mSize > 0) {
429  switch(init)
430  {
431  case 0:
432  break;
433  case 1:
434  if ( mCol == mRow ) {
435  a = mElement;
436  b = mElement + mSize;
437  for( ; a<b; a+=(mCol+1)) *a = 1.0;
438  }
439  else {
440 #ifndef ST_NO_EXCEPTIONS
441  throw domain_error("StMatrix<T>::StMatrix(): Matrix must be NxN");
442 #else
443  cerr << "StMatrix<T>::StMatrix(): Matrix must be NxN" << endl;
444 #endif
445  }
446  break;
447  default:
448 #ifndef ST_NO_EXCEPTIONS
449  throw domain_error("StMatrix<T>::StMatrix(p,q,init): init must be 0 or 1");
450 #else
451  cerr << "StMatrix<T>::StMatrix(p,q,init): init must be 0 or 1" << endl;
452 #endif
453  break;
454  }
455  }
456 }
457 
458 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
459 template<class DataType>
460 template<class X>
461 StMatrix<DataType>::StMatrix(const StMatrix<X>& m1)
462  : mRow(m1.numRow()), mCol(m1.numCol()), mSize(m1.numSize())
463 {
464  mElement = new DataType[mSize];
465 
466  for(unsigned int ii=0; ii<mRow; ii++)
467  for(unsigned int jj=0; jj<mCol; jj++)
468  *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
469 }
470 
471 template<class DataType>
472 StMatrix<DataType>::StMatrix(const StMatrix<DataType>& m1)
473  : mRow(m1.numRow()), mCol(m1.numCol()), mSize(m1.numSize())
474 {
475  mElement = new DataType[mSize];
476 
477  for(unsigned int ii=0; ii<mRow; ii++)
478  for(unsigned int jj=0; jj<mCol; jj++)
479  *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
480 }
481 #else
482 template<class DataType>
483 StMatrix<DataType>::StMatrix(const StMatrix<float>& m1)
484  : mRow(m1.numRow()), mCol(m1.numCol()), mSize(m1.numSize())
485 {
486  mElement = new DataType[mSize];
487 
488  for(unsigned int ii=0; ii<mRow; ii++)
489  for(unsigned int jj=0; jj<mCol; jj++)
490  *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
491 }
492 
493 template<class DataType>
494 StMatrix<DataType>::StMatrix(const StMatrix<double>& m1)
495  : mRow(m1.numRow()), mCol(m1.numCol()), mSize(m1.numSize())
496 {
497  mElement = new DataType[mSize];
498 
499  for(unsigned int ii=0; ii<mRow; ii++)
500  for(unsigned int jj=0; jj<mCol; jj++)
501  *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
502 }
503 
504 #endif
505 
506 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
507 template<class DataType>
508 template<class X>
509 StMatrix<DataType>& StMatrix<DataType>::operator=(const StMatrix<X>& m1)
510 {
511  if ((void *)&m1 == (void *)this) {
512  return *this;
513  }
514  else {
515  delete [] mElement;
516  mSize = m1.numRow()*m1.numCol();
517  mElement = new DataType[mSize];
518 
519  mRow = m1.numRow();
520  mCol = m1.numCol();
521 
522  for(unsigned int ii=0; ii<mRow; ii++)
523  for(unsigned int jj=0; jj<mCol; jj++) {
524  *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
525  }
526  return (*this);
527  }
528 }
529 
530 template<class DataType>
531 StMatrix<DataType>& StMatrix<DataType>::operator=(const StMatrix<DataType>& m1)
532 {
533  if ((void *)&m1 == (void *)this) {
534  return *this;
535  }
536  else {
537  delete [] mElement;
538  mSize = m1.numRow()*m1.numCol();
539  mElement = new DataType[mSize];
540 
541  mRow = m1.numRow();
542  mCol = m1.numCol();
543 
544  for(unsigned int ii=0; ii<mRow; ii++)
545  for(unsigned int jj=0; jj<mCol; jj++) {
546  *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
547  }
548  return (*this);
549  }
550 }
551 #else // ST_NO_MEMBER_TEMPLATES
552 template<class DataType>
553 StMatrix<DataType>& StMatrix<DataType>::operator=(const StMatrix<float>& m1)
554 {
555  if ((void *)&m1 == (void *)this) {
556  return *this;
557  }
558  else {
559  delete [] mElement;
560  mSize = m1.numRow()*m1.numCol();
561  mElement = new DataType[mSize];
562 
563  mRow = m1.numRow();
564  mCol = m1.numCol();
565 
566  for(unsigned int ii=0; ii<mRow; ii++)
567  for(unsigned int jj=0; jj<mCol; jj++) {
568  *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
569  }
570  return (*this);
571  }
572 }
573 template<class DataType>
574 StMatrix<DataType>& StMatrix<DataType>::operator=(const StMatrix<double>& m1)
575 {
576  if ((void *)&m1 == (void *)this) {
577  return *this;
578  }
579  else {
580  delete [] mElement;
581  mSize = m1.numRow()*m1.numCol();
582  mElement = new DataType[mSize];
583 
584  mRow = m1.numRow();
585  mCol = m1.numCol();
586 
587  for(unsigned int ii=0; ii<mRow; ii++)
588  for(unsigned int jj=0; jj<mCol; jj++) {
589  *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
590  }
591  return (*this);
592  }
593 }
594 
595 #endif
596 
597 // Destructor
598 template<class DataType>
599 StMatrix<DataType>::~StMatrix() {
600  delete [] mElement;
601 }
602 
603 //
604 // access functions
605 //
606 template<class DataType>
607 unsigned int StMatrix<DataType>::numRow() const { return mRow;}
608 
609 template<class DataType>
610 unsigned int StMatrix<DataType>::numCol() const { return mCol;}
611 
612 template<class DataType>
613 unsigned int StMatrix<DataType>::numSize() const { return mSize;}
614 
615 // backward compatibility
616 template<class DataType>
617 unsigned int StMatrix<DataType>::num_row() const { return mRow;}
618 
619 template<class DataType>
620 unsigned int StMatrix<DataType>::num_col() const { return mCol;}
621 
622 template<class DataType>
623 unsigned int StMatrix<DataType>::num_size() const { return mSize;}
624 
625 //
626 // Element Access
627 //
628 template<class DataType>
629 const DataType& StMatrix<DataType>::operator()(size_t row, size_t col) const
630 {
631 #ifdef MATRIX_BOUND_CHECK
632  if(row<1 || row>numRow() || col<1 || col>numCol()) {
633 #ifndef ST_NO_EXCEPTIONS
634  throw out_of_range("StMatrix<DataType>::operator(): const Bad Index");
635 #else
636  cerr << "StMatrix<DataType>::operator(): const Bad Index" << endl;
637 #endif
638  }
639 #endif
640  return *(mElement+(row-1)*mCol+col-1);
641 }
642 
643 template<class DataType>
644 DataType& StMatrix<DataType>::operator()(size_t row, size_t col)
645 {
646 #ifdef MATRIX_BOUND_CHECK
647  if(row<1 || row>mRow || col<1 || col>mCol) {
648 #ifndef ST_NO_EXCEPTIONS
649  throw out_of_range("StMatrix<DataType>::operator(): Bad Index");
650 #else
651  cerr << "StMatrix<DataType>::operator(): Bad Index" << endl;
652 #endif
653  }
654 #endif
655  return *(mElement+(row-1)*mCol+col-1);
656 }
657 
658 /* -----------------------------------------------------------------------
659  This section contains the assignment and inplace operators =,+=,-=,*=,/=.
660  ----------------------------------------------------------------------- */
661 
662 template<class DataType>
663 StMatrix<DataType>& StMatrix<DataType>::operator*=(double fact)
664 {
665  for(unsigned int ii=0; ii<mRow; ii++)
666  for(unsigned int jj=0; jj<mCol; jj++)
667  *(mElement+(ii)*mCol+jj) *= fact;
668 
669  return (*this);
670 }
671 
672 template<class DataType>
673 StMatrix<DataType> & StMatrix<DataType>::operator/=(double fact)
674 {
675  if(fact == 0) {
676 #ifndef ST_NO_EXCEPTIONS
677  throw out_of_range("StMatrix<T>::operator/=(): Cannot divide by zero!");
678 #else
679  cerr << "StMatrix<T>::operator/=(): Cannot divide by zero!" << endl;
680 #endif
681  }
682  for(unsigned int ii=0; ii<mCol; ii++)
683  for(unsigned int jj=0; jj<mRow; jj++)
684  *(mElement+(ii)*mCol+jj) /= fact;
685 
686  return (*this);
687 }
688 
689 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
690 template<class DataType>
691 template<class X>
692 StMatrix<DataType>& StMatrix<DataType>::operator+=(const StMatrix<X>& m2)
693 {
694  if(mRow == m2.numRow() && mCol == m2.numCol()) {
695  for(unsigned int ii=0; ii<mRow; ii++)
696  for(unsigned int jj=0; jj<mCol; jj++)
697  *(mElement+(ii)*mCol+jj) += m2(ii+1,jj+1);
698  }
699  else {
700 #ifndef ST_NO_EXCEPTIONS
701  throw out_of_range("StMatrix<T>::operator+=(): Matrices are not same size!");
702 #else
703  cerr << "StMatrix<T>::operator+=(): Matrices are not same size!" << endl;
704 #endif
705  }
706  return (*this);
707 }
708 
709 template<class DataType>
710 template<class X>
711 StMatrix<DataType>& StMatrix<DataType>::operator-=(const StMatrix<X>& m2)
712 {
713  if(mRow == m2.numRow() && mCol == m2.numCol()) {
714  for(unsigned int ii=0; ii<mRow; ii++)
715  for(unsigned int jj=0; jj<mCol; jj++)
716  *(mElement+(ii)*mCol+jj) -= m2(ii+1,jj+1);
717  }
718  else {
719 #ifndef ST_NO_EXCEPTIONS
720  throw out_of_range("StMatrix<T>::operator-=(): Matrices are not same size!");
721 #else
722  cerr << "StMatrix<T>::operator-=(): Matrices are not same size!" << endl;
723 #endif
724  }
725  return (*this);
726 }
727 
728 template<class DataType>
729 template<class X>
730 StMatrix<DataType> StMatrix<DataType>::dot(const StMatrix<X>& m2)
731 {
732  if(mCol == m2.numRow() ) {
733  StMatrix<DataType> mret(mRow, m2.numCol(), 0);
734  for(unsigned int i=0; i<mRow; i++)
735  for(unsigned int j=0; j<m2.numCol(); j++) {
736  for(unsigned int kk=0; kk<mCol; kk++)
737  mret(i+1, j+1) += (*(mElement+(i)*mCol+kk))*m2(kk+1,j+1);
738  }
739  return mret;
740  }
741  else {
742 #ifndef ST_NO_EXCEPTIONS
743  throw out_of_range("StMatrix<T>::dot(): incompatible matrix sizes");
744 #else
745  cerr << "StMatrix<T>::dot(): incompatible matrix sizes" << endl;
746 #endif
747  return StMatrix();
748  }
749 }
750 #else // ST_NO_MEMBER_TEMPLATES
751 // operator+=
752 template<class DataType>
753 StMatrix<DataType>& StMatrix<DataType>::operator+=(const StMatrix<float>& m2)
754 {
755  if(mRow == m2.numRow() && mCol == m2.numCol()) {
756  for(unsigned int ii=0; ii<mRow; ii++)
757  for(unsigned int jj=0; jj<mCol; jj++)
758  *(mElement+(ii)*mCol+jj) += m2(ii+1,jj+1);
759  }
760  else {
761 #ifndef ST_NO_EXCEPTIONS
762  throw out_of_range("StMatrix<float>::operator+=(): Matrices are not same size!");
763 #else
764  cerr << "StMatrix<float>::operator+=(): Matrices are not same size!" << endl;
765 #endif
766  }
767  return (*this);
768 }
769 
770 template<class DataType>
771 StMatrix<DataType>& StMatrix<DataType>::operator+=(const StMatrix<double>& m2)
772 {
773  if(mRow == m2.numRow() && mCol == m2.numCol()) {
774  for(unsigned int ii=0; ii<mRow; ii++)
775  for(unsigned int jj=0; jj<mCol; jj++)
776  *(mElement+(ii)*mCol+jj) += m2(ii+1,jj+1);
777  }
778  else {
779 #ifndef ST_NO_EXCEPTIONS
780  throw out_of_range("StMatrix<double>::operator+=(): Matrices are not same size!");
781 #else
782  cerr << "StMatrix<double>::operator+= Matrices are not same size!" << endl;
783 #endif
784  }
785  return (*this);
786 }
787 
788 // operator -=
789 template<class DataType>
790 StMatrix<DataType>& StMatrix<DataType>::operator-=(const StMatrix<float>& m2)
791 {
792  if(mRow == m2.numRow() && mCol == m2.numCol()) {
793  for(unsigned int ii=0; ii<mRow; ii++)
794  for(unsigned int jj=0; jj<mCol; jj++)
795  *(mElement+(ii)*mCol+jj) -= m2(ii+1,jj+1);
796  }
797  else {
798 #ifndef ST_NO_EXCEPTIONS
799  throw out_of_range("StMatrix<float>::operator-=(): Matrices are not same size!");
800 #else
801  cerr << "StMatrix<float>::operator-=(): Matrices are not same size!" << endl;
802 #endif
803  }
804  return (*this);
805 }
806 
807 template<class DataType>
808 StMatrix<DataType>& StMatrix<DataType>::operator-=(const StMatrix<double>& m2)
809 {
810  if(mRow == m2.numRow() && mCol == m2.numCol()) {
811  for(unsigned int ii=0; ii<mRow; ii++)
812  for(unsigned int jj=0; jj<mCol; jj++)
813  *(mElement+(ii)*mCol+jj) -= m2(ii+1,jj+1);
814  }
815  else {
816 #ifndef ST_NO_EXCEPTIONS
817  throw out_of_range("StMatrix<double>::operator-=(): Matrices are not same size!");
818 #else
819  cerr << "StMatrix<double>::operator-=(): Matrices are not same size!" << endl;
820 #endif
821  }
822  return (*this);
823 }
824 
825 // operator::dot
826 template<class DataType>
827 StMatrix<DataType> StMatrix<DataType>::dot(const StMatrix<float>& m2)
828 {
829  if(mCol == m2.numRow() ) {
830  StMatrix<DataType> mret(mRow, m2.numCol(), 0);
831 
832  for(unsigned int i=0; i<mRow; i++)
833  for(unsigned int j=0; j<m2.numCol(); j++) {
834  for(unsigned int kk=0; kk<mCol; kk++)
835  mret(i+1, j+1) += (*(mElement+(i)*mCol+kk))*m2(kk+1,j+1);
836  }
837  return mret;
838  }
839  else {
840 #ifndef ST_NO_EXCEPTIONS
841  throw out_of_range("StMatrix<float>::dot(): Incompatible matrix sizes");
842 #else
843  cerr << "StMatrix<float>::dot(): Incompatible matrix sizes" << endl;
844 #endif
845  return StMatrix();
846  }
847 }
848 
849 template<class DataType>
850 StMatrix<DataType> StMatrix<DataType>::dot(const StMatrix<double>& m2)
851 {
852  if(mCol == m2.numRow() ) {
853  StMatrix<DataType> mret(mRow, m2.numCol(), 0);
854 
855  for(unsigned int i=0; i<mRow; i++)
856  for(unsigned int j=0; j<m2.numCol(); j++) {
857  for(unsigned int kk=0; kk<mCol; kk++)
858  mret(i+1, j+1) += (*(mElement+(i)*mCol+kk))*m2(kk+1,j+1);
859  }
860  return mret;
861  }
862  else {
863 #ifndef ST_NO_EXCEPTIONS
864  throw out_of_range("StMatrix<double>::dot(): Incompatible matrix sizes");
865 #else
866  cerr << "StMatrix<double>::dot(): Incompatible matrix sizes" << endl;
867 #endif
868  return StMatrix();
869  }
870 }
871 
873 #endif // ST_NO_MEMBER_TEMPLATES
874 
875 template<class DataType>
876 StMatrix<DataType> StMatrix<DataType>::operator+ () const
877 {
878  return *this;
879 }
880 
881 template<class DataType>
882 StMatrix<DataType> StMatrix<DataType>::operator- () const
883 {
884  // this can be streamlined
885  StMatrix<DataType> mret(*this);
886  return mret*=-1;
887 }
888 
889 template<class DataType>
890 bool StMatrix<DataType>::operator== (const StMatrix<DataType>& m1) const
891 {
892  if (mCol == m1.numCol() && mRow == m1.numRow()) {
893  for(unsigned int ii=0; ii<mRow; ii++)
894  for(unsigned int jj=0; jj<mCol; jj++) {
895  if(*(mElement+(ii)*mCol+jj) != m1(ii+1,jj+1))
896  return false;
897  }
898  return true;
899  }
900  else
901  return false; // not even right size
902 }
903 
904 template<class DataType>
905 bool StMatrix<DataType>::operator!= (const StMatrix<DataType>& m1) const
906 {
907  return !(*this == m1);
908 }
909 
910 //Apply a function to all elements
911 template<class DataType>
912 StMatrix<DataType> StMatrix<DataType>::apply(DataType (*f)(DataType, size_t, size_t)) const
913 {
914  StMatrix<DataType> mret(mRow, mCol);
915  DataType *a = mElement;
916  for(unsigned int ir=1; ir<=mRow; ir++) {
917  for(unsigned int ic=1; ic<=mCol; ic++) {
918  mret(ir,ic) = (*f)(*(a++), ir, ic);
919  }
920  }
921  return mret;
922 }
923 
924 template<class DataType>
925 StMatrix<DataType> StMatrix<DataType>::T() const
926 {
927  StMatrix<DataType> mret(mCol,mRow);
928  DataType *pl = mElement + mSize;
929  DataType *pme = mElement;
930  DataType *pt = mret.mElement;
931  DataType *ptl = mret.mElement + mSize;
932  for (; pme < pl; pme++, pt+=mRow)
933  {
934  if (pt >= ptl)
935  pt -= (mSize-1);
936  (*pt) = (*pme);
937  }
938  return mret;
939 }
940 template<class DataType>
941 StMatrix<DataType> StMatrix<DataType>::transpose() const
942 {
943  return this->T();
944 }
945 
946 //
947 // Sub matrix
948 //
949 template<class DataType>
950 StMatrix<DataType> StMatrix<DataType>::sub(size_t min_row, size_t max_row,
951  size_t min_col,size_t max_col) const
952 {
953  StMatrix<DataType> mret(max_row-min_row+1,max_col-min_col+1);
954  if(max_row > mRow || max_col > mCol) {
955 #ifndef ST_NO_EXCEPTIONS
956  throw out_of_range("StMatrix<DataType>::sub(): Index out of range");
957 #else
958  cerr << "StMatrix<DataType>::sub(): Index out of range" << endl;
959 #endif
960  }
961  int nrows = mret.numRow();
962  int ncols = mret.numCol();
963  for(int ii=0; ii<nrows; ii++)
964  for(int jj=0; jj<ncols; jj++) {
965  mret(ii+1, jj+1) =
966  *(mElement+(min_row+ii-1)*mCol+(min_col+jj-1));
967  }
968  return mret;
969 }
970 
971 template<class DataType>
972 void StMatrix<DataType>::sub(size_t row,size_t col,const StMatrix<DataType> &m1)
973 {
974  if((row <1) ||
975  (row+m1.numRow()-1 > mRow) ||
976  (col <1) ||
977  (col+m1.numCol()-1 > mCol)) {
978 #ifndef ST_NO_EXCEPTIONS
979  throw out_of_range("StMatrix<DataType>::sub(): Index out of range");
980 #else
981  cerr << "StMatrix<DataType>::sub(): Index out of range" << endl;
982 #endif
983  }
984  DataType *a = m1.mElement;
985  unsigned int nc = numCol();
986  DataType *b1 = mElement + (row - 1) * nc + col - 1;
987 
988  for(unsigned int irow=1; irow<=m1.numRow(); irow++) {
989  DataType *brc = b1;
990  for(unsigned int icol=1; icol<=m1.numCol(); icol++) {
991  *(brc++) = *(a++);
992  }
993  b1 += nc;
994  }
995 }
996 
997 // Contents of the matrix are replaced by the inverse--less overhead then inverse().
998 template<class DataType>
999 StMatrix<DataType> StMatrix<DataType>::inverse(size_t &ierr) const
1000 {
1001  StMatrix<DataType> tmp(*this);
1002  tmp.invert(ierr);
1003  return tmp;
1004 }
1005 
1006 template<class DataType>
1007 void StMatrix<DataType>::invert(size_t &ierr) {
1008  if(mCol != mRow) {
1009 #ifndef ST_NO_EXCEPTIONS
1010  throw domain_error("StMatrix<DataType>::invert(): not a NxN matrix");
1011 #else
1012  cerr << "StMatrix<DataType>::invert(): not a NxN matrix" << endl;
1013 #endif
1014  }
1015  static unsigned int max_array = 20;
1016  static int *ir = new int [max_array+1];
1017 
1018  if (mCol > max_array) {
1019  delete [] ir;
1020  max_array = mRow;
1021  ir = new int [max_array+1];
1022  }
1023  DataType t1, t2, t3;
1024  DataType det, temp, s;
1025  int ifail;
1026  switch(mRow) {
1027  case 3:
1028  {
1029  DataType c11,c12,c13,c21,c22,c23,c31,c32,c33;
1030  ifail = 0;
1031  c11 = (*(mElement+4)) * (*(mElement+8)) - (*(mElement+5)) * (*(mElement+7));
1032  c12 = (*(mElement+5)) * (*(mElement+6)) - (*(mElement+3)) * (*(mElement+8));
1033  c13 = (*(mElement+3)) * (*(mElement+7)) - (*(mElement+4)) * (*(mElement+6));
1034  c21 = (*(mElement+7)) * (*(mElement+2)) - (*(mElement+8)) * (*(mElement+1));
1035  c22 = (*(mElement+8)) * (*mElement) - (*(mElement+6)) * (*(mElement+2));
1036  c23 = (*(mElement+6)) * (*(mElement+1)) - (*(mElement+7)) * (*mElement);
1037  c31 = (*(mElement+1)) * (*(mElement+5)) - (*(mElement+2)) * (*(mElement+4));
1038  c32 = (*(mElement+2)) * (*(mElement+3)) - (*mElement) * (*(mElement+5));
1039  c33 = (*mElement) * (*(mElement+4)) - (*(mElement+1)) * (*(mElement+3));
1040  t1 = fabs(*mElement);
1041  t2 = fabs(*(mElement+3));
1042  t3 = fabs(*(mElement+6));
1043  if (t1 >= t2) {
1044  if (t3 >= t1) {
1045  temp = *(mElement+6);
1046  det = c23*c12-c22*c13;
1047  } else {
1048  temp = *mElement;
1049  det = c22*c33-c23*c32;
1050  }
1051  } else if (t3 >= t2) {
1052  temp = *(mElement+6);
1053  det = c23*c12-c22*c13;
1054  } else {
1055  temp = *(mElement+3);
1056  det = c13*c32-c12*c33;
1057  }
1058  if (det==0) {
1059  ierr = 1;
1060  return;
1061  }
1062  {
1063  s = temp/det;
1064  DataType *tmp = mElement;
1065  *(tmp++) = s*c11;
1066  *(tmp++) = s*c21;
1067  *(tmp++) = s*c31;
1068  *(tmp++) = s*c12;
1069  *(tmp++) = s*c22;
1070  *(tmp++) = s*c32;
1071  *(tmp++) = s*c13;
1072  *(tmp++) = s*c23;
1073  *(tmp) = s*c33;
1074  }
1075  }
1076  break;
1077  case 2:
1078  ifail = 0;
1079  det = (*mElement)*(*(mElement+3)) - (*(mElement+1))*(*(mElement+2));
1080  if (det==0) {
1081  ierr = 1;
1082  return;
1083  }
1084  s = 1.0/det;
1085  temp = s*(*(mElement+3));
1086  *(mElement+1) *= -s;
1087  *(mElement+2) *= -s;
1088  *(mElement+3) = s*(*mElement);
1089  *mElement = temp;
1090  break;
1091  case 1:
1092  ifail = 0;
1093  if ((*mElement)==0) {
1094  ierr = 0;
1095  return;
1096  }
1097  *mElement = 1.0/(*mElement);
1098  break;
1099  default:
1100  ifail = dfact(det, ir);
1101  if(ifail) {
1102  ierr = 1;
1103  return;
1104  }
1105  dfinv(ir);
1106  break;
1107  }
1108  ierr = 0;
1109  return;
1110 }
1111 
1112 template<class DataType>
1113 DataType StMatrix<DataType>::determinant() const {
1114  static unsigned int max_array = 20;
1115  static int *ir = new int [max_array+1];
1116  if(mCol != mRow) {
1117 #ifndef ST_NO_EXCEPTIONS
1118  throw out_of_range("StMatrix<DataType>::determinant(): not a NxN matrix");
1119 #else
1120  cerr << "StMatrix<DataType>::determinant(): not a NxN matrix" << endl;
1121 #endif
1122  }
1123  if (mCol > max_array) {
1124  delete [] ir;
1125  max_array = mRow;
1126  ir = new int [max_array+1];
1127  }
1128  DataType det;
1129  StMatrix<DataType> mt(*this);
1130  int i = mt.dfact(det, ir);
1131  if(i==0) return det;
1132  return 0;
1133 }
1134 
1135 //
1136 // Access operators
1137 //
1138 template<class DataType>
1139 StMatrix<DataType>::StMatrixRow::StMatrixRow(StMatrix<DataType>& a, size_t r)
1140  : _a(a) {
1141  _r = r;
1142 }
1143 
1144 template<class DataType>
1145 DataType& StMatrix<DataType>::StMatrixRow::operator[](size_t c) {
1146 #ifdef MATRIX_BOUND_CHECK
1147  if (_r<0 || _r>=_a.numRow() || c<0 || c>=_a.numCol()) {
1148 #ifndef ST_NO_EXCEPTIONS
1149  throw out_of_range("StMatrix<DataType>::operator[]: index out of range");
1150 #else
1151  cerr << "StMatrix<DataType>::operator[]: index out of range" << endl;
1152 #endif
1153  }
1154 #endif
1155  return *(_a.mElement+_r*_a.mCol+c);
1156 }
1157 
1158 template<class DataType>
1159 StMatrix<DataType>::StMatrixRowConst::StMatrixRowConst(const StMatrix<DataType>& a, size_t r)
1160  : _a(a) {
1161  _r = r;
1162 }
1163 
1164 template<class DataType>
1165 const DataType& StMatrix<DataType>::StMatrixRowConst::operator[](size_t c) const
1166 {
1167 #ifdef MATRIX_BOUND_CHECK
1168  if (_r<0 || _r>=_a.numRow() || c<0 || c>=_a.numCol()) {
1169 #ifndef ST_NO_EXCEPTIONS
1170  throw out_of_range("StMatrix<DataType>::operator[]: const index out of range");
1171 #else
1172  cerr << "StMatrix<DataType>::operator[]: const index out of range" << endl;
1173 #endif
1174  }
1175 #endif
1176  return *(_a.mElement+_r*_a.mCol+c);
1177 }
1178 
1179 template<class DataType>
1180 typename StMatrix<DataType>::StMatrixRow StMatrix<DataType>::operator[] (size_t r)
1181 {
1182  StMatrixRow b(*this,r);
1183  return b;
1184 }
1185 
1186 template<class DataType>
1187 const typename StMatrix<DataType>::StMatrixRowConst StMatrix<DataType>::operator[] (size_t r) const
1188 {
1189  StMatrixRowConst b(*this,r);
1190  return b;
1191 }
1192 
1193 template<class DataType>
1194 void StMatrix<DataType>::swap(unsigned int &i, unsigned int &j)
1195 {
1196  unsigned int tmp=i;
1197  i=j;
1198  j=tmp;
1199 }
1200 
1201 template<class DataType>
1202 void StMatrix<DataType>::swap(DataType *&i, DataType *&j)
1203 {
1204  DataType *tmp=i;
1205  i=j;
1206  j=tmp;
1207 }
1208 
1209 // This function swaps two Matrices doing a full copy
1210 template<class DataType>
1211 void swap(StMatrix<DataType>& m1, StMatrix<DataType>& m2) {
1212  StMatrix<DataType> tmp(m1);
1213  m1 = m2;
1214  m2 = tmp;
1215 }
1216 
1217 // This function swaps two Matrices without doing a full copy.
1218 // The version below had some problems with the new GCC compiler
1219 // (friendship with StMatrix caused trouble) and was replaced
1220 // by the version above. tu
1221 // template<class DataType>
1222 // void swap(StMatrix<DataType>& m1, StMatrix<DataType>& m2) {
1223 // StMatrix<DataType>::swap(m1.mElement, m2.mElement);
1224 // StMatrix<DataType>::swap(m1.mRow, m2.mRow);
1225 // StMatrix<DataType>::swap(m1.mCol, m2.mCol);
1226 // StMatrix<DataType>::swap(m1.mSize, m2.mSize);
1227 // }
1228 
1229 
1230 template<class DataType>
1231 int StMatrix<DataType>::dfact(DataType &det, int *ir)
1232 {
1233  if (mCol!=mRow) {
1234 #ifndef ST_NO_EXCEPTIONS
1235  throw domain_error("StMatrix<DataType>::dfact(): Matrix not NxN");
1236 #else
1237  cerr << "StMatrix<DataType>::dfact(): Matrix not NxN" << endl;
1238 #endif
1239  }
1240 
1241  int ifail, jfail;
1242  int n = mCol;
1243 
1244  DataType tf;
1245  DataType g1 = 1.0e-19;
1246  DataType g2 = 1.0e19;
1247 
1248  DataType p, q, t;
1249  DataType s11, s12;
1250 
1251 
1252  DataType epsilon;
1253  if (sizeof(DataType) == sizeof(double))
1254  epsilon = 8*DBL_EPSILON;
1255  else
1256  epsilon = 8*FLT_EPSILON;
1257  // could be set to zero (like it was before)
1258  // but then the algorithm often doesn't detect
1259  // that a matrix is singular
1260 
1261  int normal = 0, imposs = -1;
1262  int jrange = 0, jover = 1, junder = -1;
1263  ifail = normal;
1264  jfail = jrange;
1265  int nxch = 0;
1266  det = 1.0;
1267  DataType *mj = mElement;
1268  DataType *mjj = mj;
1269  for (int j=1;j<=n;j++) {
1270  int k = j;
1271  p = (fabs(*mjj));
1272  if (j!=n) {
1273  DataType *mij = mj + n + j - 1;
1274  for (int i=j+1;i<=n;i++) {
1275  q = (fabs(*(mij)));
1276  if (q > p) {
1277  k = i;
1278  p = q;
1279  }
1280  mij += n;
1281  }
1282  if (k==j) {
1283  if (p <= epsilon) {
1284  det = 0;
1285  ifail = imposs;
1286  jfail = jrange;
1287  return ifail;
1288  }
1289  det = -det; // in this case the sign of the determinant
1290  // must not change. So I change it twice.
1291  }
1292  DataType *mjl = mj;
1293  DataType *mkl = mElement + (k-1)*n;
1294  for (int l=1;l<=n;l++) {
1295  tf = *mjl;
1296  *(mjl++) = *mkl;
1297  *(mkl++) = tf;
1298  }
1299  nxch = nxch + 1; // this makes the determinant change its sign
1300  ir[nxch] = (((j)<<12)+(k));
1301  } else {
1302  if (p <= epsilon) {
1303  det = 0.0;
1304  ifail = imposs;
1305  jfail = jrange;
1306  return ifail;
1307  }
1308  }
1309  det *= *mjj;
1310  *mjj = 1.0 / *mjj;
1311  t = (fabs(det));
1312  if (t < g1) {
1313  det = 0.0;
1314  if (jfail == jrange) jfail = junder;
1315  }
1316  else if (t > g2) {
1317  det = 1.0;
1318  if (jfail==jrange) jfail = jover;
1319  }
1320  if (j!=n) {
1321  DataType *mk = mj + n;
1322  DataType *mkjp = mk + j;
1323  DataType *mjk = mj + j;
1324  for (k=j+1;k<=n;k++) {
1325  s11 = - (*mjk);
1326  s12 = - (*mkjp);
1327  if (j!=1) {
1328  DataType *mik = mElement + k - 1;
1329  DataType *mijp = mElement + j;
1330  DataType *mki = mk;
1331  DataType *mji = mj;
1332  for (int i=1;i<j;i++) {
1333  s11 += (*mik) * (*(mji++));
1334  s12 += (*mijp) * (*(mki++));
1335  mik += n;
1336  mijp += n;
1337  }
1338  }
1339  *(mjk++) = -s11 * (*mjj);
1340  *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
1341  mk += n;
1342  mkjp += n;
1343  }
1344  }
1345  mj += n;
1346  mjj += (n+1);
1347  }
1348  if (nxch%2==1) det = -det;
1349  if (jfail !=jrange) det = 0.0;
1350  ir[n] = nxch;
1351  return 0;
1352 }
1353 
1354 template<class DataType>
1355 int StMatrix<DataType>::dfinv(int *ir)
1356 {
1357  if (mCol != mRow) {
1358 #ifndef ST_NO_EXCEPTIONS
1359  throw domain_error("StMatrix<DataType>::dfinv(): Matrix not NxN");
1360 #else
1361  cerr << "StMatrix<DataType>::dfinv(): Matrix not NxN" << endl;
1362 #endif
1363  }
1364 
1365  int n = mCol;
1366  if (n==1) return 0;
1367 
1368  DataType s31, s32;
1369  DataType s33, s34;
1370 
1371  DataType *m11 = mElement;
1372  DataType *m12 = m11 + 1;
1373  DataType *m21 = m11 + n;
1374  DataType *m22 = m12 + n;
1375  *m21 = -(*m22) * (*m11) * (*m21);
1376  *m12 = -(*m12);
1377  if (n>2) {
1378  DataType *mi = mElement + 2 * n;
1379  DataType *mii= mElement + 2 * n + 2;
1380  DataType *mimim = mElement + n + 1;
1381  for (int i=3;i<=n;i++) {
1382  int im2 = i - 2;
1383  DataType *mj = mElement;
1384  DataType *mji = mj + i - 1;
1385  DataType *mij = mi;
1386  for (int j=1;j<=im2;j++) {
1387  s31 = 0.0;
1388  s32 = *mji;
1389  DataType *mkj = mj + j - 1;
1390  DataType *mik = mi + j - 1;
1391  DataType *mjkp = mj + j;
1392  DataType *mkpi = mj + n + i - 1;
1393  for (int k=j;k<=im2;k++) {
1394  s31 += (*mkj) * (*(mik++));
1395  s32 += (*(mjkp++)) * (*mkpi);
1396  mkj += n;
1397  mkpi += n;
1398  }
1399  *mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31));
1400  *mji = -s32;
1401  mj += n;
1402  mji += n;
1403  mij++;
1404  }
1405  *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
1406  *(mimim+1) = -(*(mimim+1));
1407  mi += n;
1408  mimim += (n+1);
1409  mii += (n+1);
1410  }
1411  }
1412  DataType *mi = mElement;
1413  DataType *mii = mElement;
1414  for (int i=1;i<n;i++) {
1415  int ni = n - i;
1416  DataType *mij = mi;
1417  int j;
1418  for (j=1; j<=i;j++) {
1419  s33 = *mij;
1420  DataType *mikj = mi + n + j - 1;
1421  DataType *miik = mii + 1;
1422  DataType *min_end = mi + n;
1423  for (;miik<min_end;) {
1424  s33 += (*mikj) * (*(miik++));
1425  mikj += n;
1426  }
1427  *(mij++) = s33;
1428  }
1429  for (j=1;j<=ni;j++) {
1430  s34 = 0.0;
1431  DataType *miik = mii + j;
1432  DataType *mikij = mii + j * n + j;
1433  for (int k=j;k<=ni;k++) {
1434  s34 += *mikij * (*(miik++));
1435  mikij += n;
1436  }
1437  *(mii+j) = s34;
1438  }
1439  mi += n;
1440  mii += (n+1);
1441  }
1442  int nxch = ir[n];
1443  if (nxch==0) return 0;
1444  for (int mm=1;mm<=nxch;mm++) {
1445  int k = nxch - mm + 1;
1446  int ij = ir[k];
1447  int i = ij >> 12;
1448  int j = ij%4096;
1449  DataType *mki = mElement + i - 1;
1450  DataType *mkj =mElement + j - 1;
1451  for (k=1; k<=n;k++) {
1452  DataType ti = *mki;
1453  *mki = *mkj;
1454  *mkj = ti;
1455  mki += n;
1456  mkj += n;
1457  }
1458  }
1459  return 0;
1460 }
1461 
1462 
1463 
1464 template<class DataType>
1465 StMatrix<DataType> operator/(const StMatrix<DataType>& m1, double fact)
1466 {
1467  StMatrix<DataType> mret(m1);
1468  mret /= fact;
1469  return mret;
1470 }
1471 
1472 template<class DataType>
1473 StMatrix<DataType> operator*(const StMatrix<DataType>& m1, double fact)
1474 {
1475  StMatrix<DataType> mret(m1);
1476  mret *= fact;
1477  return mret;
1478 }
1479 
1480 template<class DataType>
1481 StMatrix<DataType> operator*(double fact,const StMatrix<DataType>& m1)
1482 {
1483  StMatrix<DataType> mret(m1);
1484  mret *= fact;
1485  return mret;
1486 }
1487 
1488 // Private Members
1489 // Direct sum of two matricies
1490 template<class DataType>
1491 StMatrix<DataType> dsum(const StMatrix<DataType> &m1, const StMatrix<DataType> &m2)
1492 {
1493  StMatrix<DataType> mret(m1.numRow() + m2.numRow(), m1.numCol() + m2.numCol());
1494  mret.sub(1,1,m1);
1495  mret.sub(m1.numRow()+1, m1.numCol()+1, m2);
1496  return mret;
1497 }
1498 
1499 #endif /* ! __CINT__ */
1500 #ifdef __CINT__
1501 template<> StMatrix<double> operator*(const StMatrix<double>& m1,const StMatrix<double>& m2);
1502 template<> StMatrix<double> operator*(const StMatrix<double>& m1,const StMatrix<float>& m2);
1503 template<> StMatrix<double> operator*(const StMatrix<float>& m1,const StMatrix<double>& m2);
1504 template<> StMatrix<float> operator*(const StMatrix<float>& m1,const StMatrix<float>& m2);
1505 
1506 template<> vector<double> operator*(const StMatrix<double>& m1, const vector<double>& v);
1507 template<> vector<double> operator*(const vector<double>& v, const StMatrix<double>& m1);
1508 template<> vector<double> operator*(const StMatrix<double>& m1, const vector<float>& v);
1509 template<> vector<double> operator*(const vector<float>& v, const StMatrix<double>& m1);
1510 template<> vector<double> operator*(const StMatrix<float>& m1, const vector<double>& v);
1511 template<> vector<double> operator*(const vector<double>& v, const StMatrix<float>& m1);
1512 template<> vector<float> operator*(const StMatrix<float>& m1, const vector<float>& v);
1513 template<> vector<float> operator*(const vector<float>& v, const StMatrix<float>& m1);
1514 
1515 template<> StThreeVector<double> operator*(const StMatrix<double>& m1, const StThreeVector<double>& v);
1516 template<> StThreeVector<double> operator*(const StThreeVector<double>& v, const StMatrix<double>& m1);
1517 template<> StThreeVector<double> operator*(const StMatrix<double>& m1, const StThreeVector<float>& v);
1518 template<> StThreeVector<double> operator*(const StThreeVector<float>& v, const StMatrix<double>& m1);
1519 template<> StThreeVector<double> operator*(const StMatrix<float>& m1, const StThreeVector<double>& v);
1520 template<> StThreeVector<double> operator*(const StThreeVector<double>& v, const StMatrix<float>& m1);
1521 template<> StThreeVector<float> operator*(const StMatrix<float>& m1, const StThreeVector<float>& v);
1522 template<> StThreeVector<float> operator*(const StThreeVector<float>& v, const StMatrix<float>& m1);
1523 
1524 template<> StLorentzVector<double> operator*(const StMatrix<double>& m1, const StLorentzVector<double>& v);
1525 template<> StLorentzVector<double> operator*(const StLorentzVector<double>& v, const StMatrix<double>& m1);
1526 template<> StLorentzVector<double> operator*(const StMatrix<double>& m1, const StLorentzVector<float>& v);
1527 template<> StLorentzVector<double> operator*(const StLorentzVector<float>& v, const StMatrix<double>& m1);
1528 template<> StLorentzVector<double> operator*(const StMatrix<float>& m1, const StLorentzVector<double>& v);
1529 template<> StLorentzVector<double> operator*(const StLorentzVector<double>& v, const StMatrix<float>& m1);
1530 template<> StLorentzVector<float> operator*(const StMatrix<float>& m1, const StLorentzVector<float>& v);
1531 template<> StLorentzVector<float> operator*(const StLorentzVector<float>& v, const StMatrix<float>& m1);
1532 
1533 template<> StMatrix<double> operator+(const StMatrix<double>& m1,const StMatrix<double>& m2);
1534 template<> StMatrix<double> operator+(const StMatrix<float>& m1,const StMatrix<double>& m2);
1535 template<> StMatrix<double> operator+(const StMatrix<double>& m1,const StMatrix<float>& m2);
1536 template<> StMatrix<float> operator+(const StMatrix<float>& m1,const StMatrix<float>& m2);
1537 template<> StMatrix<double> operator-(const StMatrix<double>& m1,const StMatrix<double>& m2);
1538 template<> StMatrix<double> operator-(const StMatrix<float>& m1,const StMatrix<double>& m2);
1539 template<> StMatrix<double> operator-(const StMatrix<double>& m1,const StMatrix<float>& m2);
1540 template<> StMatrix<float> operator-(const StMatrix<float>& m1,const StMatrix<float>& m2);
1541 
1542 template<> ostream& operator<<(ostream& s, const StMatrix<double>& q);
1543 template<> ostream& operator<<(ostream& s, const StMatrix<float>& q);
1544 
1545 template<> double norm_infinity(const StMatrix<double>& m1);
1546 template<> double normInfinity(const StMatrix<double>& m1);
1547 template<> double norm1(const StMatrix<double>& m1);
1548 template<> float norm_infinity(const StMatrix<float>& m1);
1549 template<> float normInfinity(const StMatrix<float>& m1);
1550 template<> float norm1(const StMatrix<float>& m1);
1551 #else /* ! __CINT__ */
1552 // Non-Member
1553 template<class DataType, class X>
1554 StMatrix<DataType> operator*(const StMatrix<DataType>& m1,const StMatrix<X>& m2)
1555 {
1556  return StMatrix<DataType>(m1).dot(m2);
1557 }
1558 
1559 // STL vector
1560 template <class DataType, class X>
1561 #ifdef ST_NO_TEMPLATE_DEF_ARGS
1562 vector<DataType, allocator<DataType> >
1563 operator*(const StMatrix<X>& m1, const vector<DataType, allocator<DataType> >& v)
1564 #else
1565 vector<DataType> operator*(const StMatrix<X>& m1, const vector<DataType>& v)
1566 #endif
1567 {
1568  vector<DataType> mret(m1.numRow());
1569  if(m1.numCol() == v.size() ) {
1570  for (unsigned int i=0; i<v.size(); i++) {
1571  DataType temp = 0;
1572  for(unsigned int j=0; j<m1.numCol(); j++)
1573  temp += m1(i+1,j+1)*v[j];
1574  mret[i]=temp;
1575  }
1576  }
1577  else {
1578 #ifndef ST_NO_EXCEPTIONS
1579  throw out_of_range("operator*(): StMatrix * STL vector Sizes are incompatible");
1580 #else
1581  cerr << "operator*(): StMatrix * STL vector Sizes are incompatible" << endl;
1582 #endif
1583  }
1584  return mret;
1585 }
1586 
1587 template <class DataType, class X>
1588 #ifdef ST_NO_TEMPLATE_DEF_ARGS
1589 vector<DataType, allocator<DataType> >
1590 operator*(const vector<DataType, allocator<DataType> >& v, const StMatrix<X>& m1)
1591 #else
1592 vector<DataType> operator*(const vector<DataType>& v, const StMatrix<X>& m1)
1593 #endif
1594 {
1595 
1596  if(v.size() == m1.numRow()) {
1597  vector<DataType> mret(m1.numCol());
1598  for (unsigned int i=0; i<m1.numCol(); i++) {
1599  DataType temp = 0;
1600  for(unsigned int j=0; j<m1.numRow(); j++)
1601  temp+=v[j]*m1(j+1,i+1);
1602  mret[i]=temp;
1603  }
1604  return mret;
1605  }
1606  else {
1607 #ifndef ST_NO_EXCEPTIONS
1608  throw out_of_range("operator*(): STL vector * StMatrix : Matrix Sizes are incompatible");
1609 #else
1610  cerr << "operator*(): STL vector * StMatrix : Matrix Sizes are incompatible" << endl;
1611 #endif
1612  return vector<DataType>();
1613  }
1614 }
1615 
1616 template <class DataType, class X>
1617 StThreeVector<DataType> operator*(const StMatrix<X>& m1, const StThreeVector<DataType>& v3)
1618 {
1619  if(m1.numRow() == 3 && m1.numCol() == 3) {
1620  return StThreeVector<DataType>(m1[0][0]*v3.x()+m1[0][1]*v3.y()+m1[0][2]*v3.z(),
1621  m1[1][0]*v3.x()+m1[1][1]*v3.y()+m1[1][2]*v3.z(),
1622  m1[2][0]*v3.x()+m1[2][1]*v3.y()+m1[2][2]*v3.z());
1623  }
1624  else {
1625 #ifndef ST_NO_EXCEPTIONS
1626  throw out_of_range("operator*(): StMatrix<> * StThreeVector<> : Matrix Must be 3x3.");
1627 
1628 #else
1629  cerr << "StMatrix<> * StThreeVector<>: Matrix Must be 3x3" << endl;
1630 #endif
1631  return StThreeVector<DataType>();
1632  }
1633 }
1634 
1635 
1636 template <class DataType, class X>
1637 StThreeVector<DataType> operator*(const StThreeVector<DataType>& v3, const StMatrix<X>& m1)
1638 {
1639  if(m1.numRow() == 3 && m1.numCol() == 3) {
1640  return StThreeVector<DataType>(m1[0][0]*v3.x()+m1[1][0]*v3.y()+m1[2][0]*v3.z(),
1641  m1[0][1]*v3.x()+m1[1][1]*v3.y()+m1[2][1]*v3.z(),
1642  m1[0][2]*v3.x()+m1[1][2]*v3.y()+m1[2][2]*v3.z());
1643  }
1644  else {
1645 #ifndef ST_NO_EXCEPTIONS
1646  throw out_of_range("operator*(): StThreeVector<> * StMatrix<>: Matrix Must be 3x3.");
1647 
1648 #else
1649  cerr << "operator*(): StThreeVector<> * StMatrix<>: Matrix Must be 3x3" << endl;
1650 #endif
1651  return StThreeVector<DataType>();
1652  }
1653 }
1654 
1655 template <class DataType, class X>
1656 StLorentzVector<DataType> operator*(const StMatrix<X>& m1, const StLorentzVector<DataType>& v4)
1657 {
1658  if (m1.numRow() == 4 && m1.numCol() == 4) {
1659  return StLorentzVector<DataType>(m1[0][0]*v4.x()+m1[0][1]*v4.y()+m1[0][2]*v4.z()+m1[0][3]*v4.t(),
1660  m1[1][0]*v4.x()+m1[1][1]*v4.y()+m1[1][2]*v4.z()+m1[1][3]*v4.t(),
1661  m1[2][0]*v4.x()+m1[2][1]*v4.y()+m1[2][2]*v4.z()+m1[2][3]*v4.t(),
1662  m1[3][0]*v4.x()+m1[3][1]*v4.y()+m1[3][2]*v4.z()+m1[3][3]*v4.t());
1663  }
1664  else {
1665 #ifndef ST_NO_EXCEPTIONS
1666  throw out_of_range("operator*(): StMatrix<> * StLorentzVector<> : Matrix Must be 4x4.");
1667 
1668 #else
1669  cerr << "operator*(): StMatrix<> * StLorentzVector<>: Matrix Must be 4x4" << endl;
1670 #endif
1671  return StLorentzVector<DataType>();
1672  }
1673 }
1674 
1675 template <class DataType, class X>
1676 StLorentzVector<DataType> operator*(const StLorentzVector<DataType>& v4, const StMatrix<X>& m1)
1677 {
1678  if (m1.numRow() == 4 && m1.numCol() == 4) {
1679  return StLorentzVector<DataType>(m1[0][0]*v4.x()+m1[1][0]*v4.y()+m1[2][0]*v4.z()+m1[0][3]*v4.t(),
1680  m1[0][1]*v4.x()+m1[1][1]*v4.y()+m1[2][1]*v4.z()+m1[1][3]*v4.t(),
1681  m1[0][2]*v4.x()+m1[1][2]*v4.y()+m1[2][2]*v4.z()+m1[2][3]*v4.t(),
1682  m1[0][3]*v4.x()+m1[1][3]*v4.y()+m1[2][3]*v4.z()+m1[2][3]*v4.t());
1683  }
1684  else {
1685 #ifndef ST_NO_EXCEPTIONS
1686  throw out_of_range("operator*(): StLorentzVector<> * StMatrix<>: Matrix Must be 3x3.");
1687 
1688 #else
1689  cerr << "StLorentzVector<> * StMatrix<>: Matrix Must be 3x3" << endl;
1690 #endif
1691  return StLorentzVector<DataType>();
1692  }
1693 }
1694 
1695 template<class DataType, class X>
1696 StMatrix<DataType> operator+(const StMatrix<DataType>& m1,const StMatrix<X>& m2)
1697 {
1698  if(m1.numRow() == m2.numRow() && m1.numCol() == m2.numCol()) {
1699  StMatrix<DataType> mret(m1);
1700  mret +=m2;
1701  return mret;
1702  }
1703  else {
1704 #ifndef ST_NO_EXCEPTIONS
1705  throw out_of_range("operator+(): Matrix Sizes must be the same.");
1706 
1707 #else
1708  cerr << "operator+(): Matrix Sizes must be the same." << endl;
1709 #endif
1710  return StMatrix<DataType>();
1711  }
1712 }
1713 
1714 template<class DataType, class X>
1715 StMatrix<DataType> operator-(const StMatrix<DataType>& m1,const StMatrix<X>& m2)
1716 {
1717  if(m1.numRow() == m2.numRow() && m1.numCol() == m2.numCol()) {
1718  return StMatrix<DataType>(m1) -= m2;
1719  }
1720  else {
1721 #ifndef ST_NO_EXCEPTIONS
1722  throw out_of_range("operator-(): Matrix Sizes must be the same.");
1723 
1724 #else
1725  cerr << "operator-(): Matrix Sizes must be the same." << endl;
1726 #endif
1727  return StMatrix<DataType>();
1728  }
1729 }
1730 
1731 // Print the Matrix.
1732 template<class DataType>
1733 ostream& operator<<(ostream& s, const StMatrix<DataType>& q)
1734 {
1735  s << "\n";
1736  // Fixed format needs 3 extra characters for field
1737  // Scientific format needs 7
1738  unsigned int width;
1739  if(s.flags()&std::ios::fixed)
1740  width = s.precision()+3;
1741  else
1742  width = s.precision()+7;
1743  for(unsigned int irow = 1; irow<= q.numRow(); irow++)
1744  {
1745  for(unsigned int icol=1; icol<=q.numCol(); icol++)
1746  {
1747  s.width(width);
1748  s << q(irow,icol) << " ";
1749  }
1750  s<< endl;
1751  }
1752  return s;
1753 }
1754 
1755 template<class DataType>
1756 DataType norm_infinity(const StMatrix<DataType>& m1)
1757 {
1758  return normInfinity(m1);
1759 }
1760 
1761 template<class DataType>
1762 DataType normInfinity(const StMatrix<DataType>& m1)
1763 {
1764  DataType max=0,sum;
1765  for(unsigned int r=1; r<=m1.numRow(); r++) {
1766  sum=0;
1767  for(unsigned int c=1; c<=m1.numCol(); c++) {
1768  sum+=fabs(m1(r,c));
1769  }
1770  if(sum>max) max=sum;
1771  }
1772  return max;
1773 }
1774 
1775 template<class DataType>
1776 DataType norm1(const StMatrix<DataType>& m1)
1777 {
1778  DataType max=0,sum;
1779  for(unsigned int c=1; c<=m1.numCol(); c++) {
1780  sum=0;
1781  for(unsigned int r=1; r<=m1.num_row(); r++)
1782  sum+=fabs(m1(r,c));
1783  if(sum>max) max=sum;
1784  }
1785  return max;
1786 }
1787 #endif /* __CINT__ */
1788 #endif /* ST_MATRIX_HH */
Generated by   doxygen 1.8.5