StRoot: StarClassLibrary/StMatrix.hh Source File

00001 /***************************************************************************
00002  *
00003  * $Id: StMatrix.hh,v 1.19 2006/04/07 22:02:34 ullrich Exp $
00004  *
00005  * Author: Original code from CLHEP by Mike Smyth
00006  *         Modified April 17, 1998 Brian Lasiuk (templated version)
00007  *
00008  ***************************************************************************
00009  *
00010  * Description:
00011  * StMatrix class, does basic matrix operations
00012  *
00013  * Remarks:   Since not all compilers support member templates
00014  *            we have to specialize the templated member on these
00015  *            platforms. If member templates are not supported the
00016  *            ST_NO_MEMBER_TEMPLATES flag has to be set. tu.
00017  *
00018  ***************************************************************************
00019  *
00020  * $Log: StMatrix.hh,v $
00021  * Revision 1.19  2006/04/07 22:02:34  ullrich
00022  * Fixed bug in dfinv and dfact that were affecting
00023  * determinant(), invert(), and inverse().
00024  *
00025  * Revision 1.18  2006/01/09 23:47:27  fisyak
00026  * Add missing methods (found by Zhangbu) to Cint dictionary
00027  *
00028  * Revision 1.17  2005/09/22 20:09:20  fisyak
00029  * Make StLorentzVector persistent
00030  *
00031  * Revision 1.16  2005/07/06 18:49:56  fisyak
00032  * Replace StHelixD, StLorentzVectorD,StLorentzVectorF,StMatrixD,StMatrixF,StPhysicalHelixD,StThreeVectorD,StThreeVectorF by templated version
00033  *
00034  * Revision 1.15  2004/01/14 22:37:27  fisyak
00035  *  unsigned int => size_t to make alpha happy
00036  *
00037  * Revision 1.14  2003/09/02 17:59:35  perev
00038  * gcc 3.2 updates + WarnOff
00039  *
00040  * Revision 1.13  2002/07/02 22:20:53  ullrich
00041  * Add dummy return statements to get rid of warnings under Linux.
00042  *
00043  * Revision 1.12  2001/10/31 15:11:36  ullrich
00044  * Rewrote swap() to work as non-friend to StMatrix.
00045  *
00046  * Revision 1.11  2001/10/31 00:33:34  ullrich
00047  * Remove macro ifdef for GCC which is not needed anymore.
00048  *
00049  * Revision 1.10  2000/02/02 18:31:05  lasiuk
00050  * restore files
00051  *
00052  * Revision 1.9  2000/02/01 16:02:59  lasiuk
00053  * namespace std is different on SUN CC5 and KCC.  Redefine macros!
00054  *
00055  * Revision 1.8  2000/01/31 20:53:45  lasiuk
00056  * using std::swap
00057  *
00058  * Revision 1.7  1999/12/21 15:14:13  ullrich
00059  * Modified to cope with new compiler version on Sun (CC5.0).
00060  *
00061  * Revision 1.6  1999/12/07 23:43:04  ullrich
00062  * Modified to get rid of warnings on Linux.
00063  *
00064  * Revision 1.5  1999/04/14 23:12:07  fisyak
00065  * Add __CINT__ to handle references
00066  *
00067  * Revision 1.4  1999/03/04 18:12:24  ullrich
00068  * Added namespace 'std'.
00069  *
00070  * Revision 1.3  1999/02/17 11:38:54  ullrich
00071  * Removed specialization for 'long double'.
00072  *
00073  * Revision 1.2  1999/02/14 23:11:43  fisyak
00074  * Fixes for Rootcint
00075  *
00076  * Revision 1.1  1999/01/30 03:59:03  fisyak
00077  * Root Version of StarClassLibrary
00078  *
00079  * Revision 1.1  1999/01/23 00:27:56  ullrich
00080  * Initial Revision
00081  *
00082  **************************************************************************/
00083 // This code was modified from the CLHEP files:
00084 // GenMatrixD.h/.icc/.cc MatrixD.h/.icc/.cc
00085 // Comments that were kept from the original code are preceding
00086 // by ///.
00087 
00107 //  The Matrix class provides the tools to do matrix arithmetic.
00108 //  The Matrices are templates with the default type being "double"
00109 //  A Matrix is declared by:
00110 //
00111 //   StMatrix<double> m1(n,m);    // n x m matrix
00112 //   StMatrix<> m2(n,m);          // default type is double
00113 //
00114 //  A Matrix may be initialized by giving it a third integer
00115 //  argument, either 0 or 1:
00116 //  - Zero (0) means initialize all elements to 0.
00117 //  - One (1) means initialize to the identity matrix if the matrix is
00118 //    square (NxN); otherwise an exception is thrown.
00119 //  - The default behavior is to initialize all elements to zero (0)
00120 //
00121 //   StMatrix<double> m2(n, m, 0);  // all zeros
00122 //   StMatrix<double> m1(n, n, 1);  // diagnol elements are 1 
00123 //   StMatrix<double> m1(n, m, 1);  // error if m!=n
00124 //
00125 //  Copying and Assignment:
00126 //
00127 //   StMatrix<double> m1(m2);      // m2 copied to m1
00128 //   StMatrix<double> m1 = m2;     // m1 assigned value of m1
00129 //
00130 //  Element Access
00131 //  Single elements of the matrix can be accessed in two differernt ways:
00132 //  - C array style operators []
00133 //
00134 //  For a pXq matrix:
00135 //  --> m1[i][j]:
00136 //       0 <= i <= (p-1)
00137 //       0 <= j <= (q-1)
00138 //
00139 //  - Fortran array style operators ()
00140 //  --> m1(i,j):
00141 //       1 <= i <= p
00142 //       1 <= j <= q
00143 //
00150 //
00152 //   #define MATRIX_BOUND_CHECK
00153 //   before including StMatrix.hh
00155 //  Access Functions:
00156 //  Properties of the matrices can be extracted via the following
00157 //  access functions:
00158 //
00159 //   m1.num_row();   // number of rows
00160 //   m1.numRow();    // "
00161 //   m1.num_col();   // number of columns
00162 //   m1.numCol();    // "
00163 //
00165 //  - Scaler operations:
00166 //    A matrix may be multiplied or divided by a scaler
00167 //    via the /= and *= operators:
00168 //    m1 /= 5;      // divides all elements in matrix m1 by 5
00169 //
00170 //  - Matrix operations:
00171 //    Two matrices can be added, subtracted, and multiplied with the
00172 //    operators +=, -=, *=
00173 //
00174 //    m1 += m2;   // same as m1 = m1+m2;
00175 //
00180 //   mt = m.transpose() ;    // "
00181 //
00182 //   unsigned int ierr;      // unsigned int == unsigned int
00184 //                           // ierr will be non-zero if the matrix is singular.
00186 //                           // faster (avoids copy operator)
00189 //                                                    //  m2(row_min:row_max, col_min:col_max)
00191 //                           //  is replaced with m1.
00194 //                //  applies f to every element of m1 and places it in m.
00195 //
00198 //
00199 //  There also exists functionality to multiply a
00200 //  1) StThreeVector
00201 //  2) StLorentzVector
00202 //  3) Std C++ Library vector<>
00203 //  by a matrix.
00204 //
00205 //  StThreeVector<> a(1,1,1);
00206 //  StMatrix<> M1(3,3,1);
00207 //  StThreeVector<> b = a*M1;   // multiplication by identity matrix
00208 //
00209 //  This can be used for simple rotations.  The same syntax is also used
00210 //  for the LorentzVector and Std C++ vector.  In order to use the
00211 //  StThreeVector and StLorentzVector capablilties simply define:
00212 //
00213 //  #define WITH_ST_THREEVECTOR     // for use with StThreeVector
00214 //  #define WITH_ST_LORENTZVECTOR   // for use with StLorentzVector
00223 
00224 #ifndef ST_MATRIX_HH
00225 #define ST_MATRIX_HH
00226 
00227 #include <Stiostream.h>
00228 #include <math.h>
00229 #include <vector>
00230 #include <float.h>
00231 #if !defined(ST_NO_NAMESPACES)
00232 using std::vector;
00233 #endif
00234 
00235 #ifndef ST_NO_EXCEPTIONS
00236 #   include <stdexcept>
00237 #   if !defined(ST_NO_NAMESPACES)
00238         using std::out_of_range;
00239         using std::domain_error;
00240 #   endif
00241 #endif
00242 
00243 #include "StThreeVector.hh"
00244 #include "StLorentzVector.hh"
00245 
00246 // in types.h assumed:
00247 // typedef   unsigned int  size_t;
00248 
00249 template<class DataType> class StMatrix {
00250 public:
00251     StMatrix();
00252     StMatrix(size_t p, size_t q, size_t init=0);
00253 
00254     // Copy constructor.
00255 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
00256     template<class X>
00257     StMatrix(const StMatrix<X>&);
00258     StMatrix(const StMatrix<DataType>&);
00259 #else
00260     StMatrix(const StMatrix<float>&);
00261     StMatrix(const StMatrix<double>&);
00262 #endif
00263     
00264     // Assignment operators.
00265 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
00266     template<class X>
00267     StMatrix<DataType>& operator=(const StMatrix<X>&);
00268     StMatrix<DataType>& operator=(const StMatrix<DataType>&);
00269 #else
00270     StMatrix<DataType>& operator=(const StMatrix<float>&);
00271     StMatrix<DataType>& operator=(const StMatrix<double>&);
00272 #endif
00273     
00274     // Destructor. --
00275     //Problem with LINUX (virtual table error when exception thrown)
00276     virtual ~StMatrix();
00277 
00278     //
00279     // access functions
00280     //
00281     unsigned int numRow()  const;
00282     unsigned int numCol()  const;
00283     unsigned int numSize() const;
00284     unsigned int num_row()  const;
00285     unsigned int num_col()  const;
00286     unsigned int num_size() const;
00287     //
00288     // Element access:
00289     // ** Note that the indexing starts from:
00290     // For (,)   --> (1,1). **
00291     // For [][] --> [0][0]
00292     const DataType& operator()(size_t row, size_t col) const;
00293     DataType& operator()(size_t row, size_t col);
00294 
00295     // classes for implementing m[i][j]
00296     class StMatrixRow {
00297     public:
00298         StMatrixRow(StMatrix<DataType>&, size_t);
00299         DataType& operator[](size_t);
00300     private:
00301 #ifndef __CINT__
00302         StMatrix<DataType>& _a;
00303 #else
00304         StMatrix<DataType>* _a;
00305 #endif
00306         size_t _r;
00307     };
00308     class StMatrixRowConst {
00309     public:
00310         StMatrixRowConst (const StMatrix<DataType>&, size_t);
00311         const DataType & operator[](size_t) const;
00312     private:
00313 #ifndef __CINT__
00314         const StMatrix<DataType>& _a;
00315 #else
00316         const StMatrix<DataType>* _a;
00317 #endif
00318         size_t _r;
00319     };
00320 
00321     StMatrixRow operator[] (size_t r);
00322     const StMatrixRowConst operator[] (size_t) const;
00323     
00324     // Scaler Operations
00325     StMatrix<DataType>& operator*=(double t);
00326     StMatrix<DataType>& operator/=(double t);
00327 
00328     // Matrix Operations
00329 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
00330     template<class X> StMatrix<DataType>& operator+=(const StMatrix<X>&);
00331     template<class X> StMatrix<DataType>& operator-=(const StMatrix<X>&);
00332     template<class X> StMatrix<DataType>  dot(const StMatrix<X>&);
00333 #else
00334     StMatrix<DataType>& operator+=(const StMatrix<float>&);
00335     StMatrix<DataType>& operator+=(const StMatrix<double>&);
00336     
00337     StMatrix<DataType>& operator-=(const StMatrix<float>&);
00338     StMatrix<DataType>& operator-=(const StMatrix<double>&);
00339 
00340     StMatrix<DataType>  dot(const StMatrix<float>&);
00341     StMatrix<DataType>  dot(const StMatrix<double>&);
00342 #endif
00343     
00344     // unary Operations
00345     StMatrix<DataType> operator+ () const;
00346     StMatrix<DataType> operator- () const;
00347 
00348     // Equality
00349     bool operator==(const StMatrix<DataType>&) const;
00350     bool operator!=(const StMatrix<DataType>&) const;
00351     
00352     // Apply a function to all elements of the matrix.
00353     StMatrix<DataType> apply(DataType (*f)(DataType, size_t, size_t)) const;
00354 
00355     StMatrix<DataType> T() const;
00356     StMatrix<DataType> transpose() const;
00357 
00358     StMatrix<DataType> sub(size_t min_row, size_t max_row, size_t min_col, size_t max_col) const;
00359     
00360     void sub(size_t row, size_t col, const StMatrix<DataType> &m1);
00361 
00362     StMatrix<DataType> inverse(size_t& ierr) const;
00363     void invert(size_t& ierr);
00364 
00365     DataType determinant() const;
00366     
00367     // Must be of same type to swap
00368     static void   swap(unsigned int&, unsigned int&);
00369     static void   swap(DataType *&, DataType *&);
00370     
00371 private:
00372     // Friend classes.
00373     friend class  StMatrixRow;
00374     friend class  StMatrixRowConst;
00375 
00376     // If successful, the return code is 0.
00377     // On return, det is the determinant and ir[] is row-unsigned interchange
00378     // matrix. See CERNLIB's DFACT routine.
00379     int dfact(DataType &det, int *ir); // factorize the matrix. 
00380 
00381     // invert the matrix. See CERNLIB DFINV.
00382     int dfinv(int *ir);
00383 
00384 protected:
00385     DataType *mElement;
00386     unsigned int mRow, mCol;
00387     unsigned int mSize;
00388 #ifdef __ROOT__
00389   ClassDef(StMatrix,3)
00390 #endif /* __ROOT__ */
00391 };
00392 
00393 #ifndef __CINT__
00394 // Constructors. 
00395 
00396 template<class DataType>
00397 inline StMatrix<DataType>::StMatrix()
00398     : mElement(0), mRow(0), mCol(0), mSize(0) {}
00399 
00400 template<class DataType>
00401 StMatrix<DataType>::StMatrix(size_t p,size_t q, size_t init)
00402     : mRow(p), mCol(q)
00403 {
00404     mSize    = mRow*mCol;
00405     mElement = new DataType[mSize];
00406 
00407     DataType *a = mElement;
00408     DataType *b = mElement + mSize;
00409     for( ; a<b; a++)
00410         *a = 0;
00411     
00412     if (mSize > 0) {
00413         switch(init)
00414             {
00415             case 0:
00416                 break;
00417             case 1:
00418                 if ( mCol == mRow ) {
00419                     DataType *a = mElement;
00420                     DataType *b = mElement + mSize;
00421                     for( ; a<b; a+=(mCol+1)) *a = 1.0;
00422                 }
00423                 else {
00424 #ifndef ST_NO_EXCEPTIONS
00425                     throw domain_error("StMatrix<T>::StMatrix(): Matrix must be NxN");
00426 #else
00427                     cerr << "StMatrix<T>::StMatrix(): Matrix must be NxN" << endl;
00428 #endif
00429                 }
00430                 break;
00431             default:
00432 #ifndef ST_NO_EXCEPTIONS
00433                 throw domain_error("StMatrix<T>::StMatrix(p,q,init): init must be 0 or 1");
00434 #else
00435                 cerr << "StMatrix<T>::StMatrix(p,q,init): init must be 0 or 1" << endl;
00436 #endif
00437                 break;
00438             }
00439     }
00440 }
00441 
00442 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
00443 template<class DataType>
00444 template<class X>
00445 StMatrix<DataType>::StMatrix(const StMatrix<X>& m1)
00446     : mRow(m1.numRow()), mCol(m1.numCol()), mSize(m1.numSize())
00447 {
00448     mElement = new DataType[mSize];
00449 
00450     for(unsigned int ii=0; ii<mRow; ii++)
00451         for(unsigned int jj=0; jj<mCol; jj++)
00452             *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
00453 }
00454 
00455 template<class DataType>
00456 StMatrix<DataType>::StMatrix(const StMatrix<DataType>& m1)
00457     : mRow(m1.numRow()), mCol(m1.numCol()), mSize(m1.numSize())
00458 {
00459     mElement = new DataType[mSize];
00460 
00461     for(unsigned int ii=0; ii<mRow; ii++)
00462         for(unsigned int jj=0; jj<mCol; jj++)
00463             *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
00464 }
00465 #else
00466 template<class DataType>
00467 StMatrix<DataType>::StMatrix(const StMatrix<float>& m1)
00468     : mRow(m1.numRow()), mCol(m1.numCol()), mSize(m1.numSize())
00469 {
00470     mElement = new DataType[mSize];
00471 
00472     for(unsigned int ii=0; ii<mRow; ii++)
00473         for(unsigned int jj=0; jj<mCol; jj++)
00474             *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
00475 }
00476 
00477 template<class DataType>
00478 StMatrix<DataType>::StMatrix(const StMatrix<double>& m1)
00479     : mRow(m1.numRow()), mCol(m1.numCol()), mSize(m1.numSize())
00480 {
00481     mElement = new DataType[mSize];
00482 
00483     for(unsigned int ii=0; ii<mRow; ii++)
00484         for(unsigned int jj=0; jj<mCol; jj++)
00485             *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
00486 }
00487 
00488 #endif
00489 
00490 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
00491 template<class DataType>
00492 template<class X>
00493 StMatrix<DataType>& StMatrix<DataType>::operator=(const StMatrix<X>& m1)
00494 {   
00495     if ((void *)&m1 == (void *)this) {
00496         return *this;
00497     }
00498     else {
00499         delete [] mElement;
00500         mSize    = m1.numRow()*m1.numCol();
00501         mElement = new DataType[mSize];
00502 
00503         mRow = m1.numRow();
00504         mCol = m1.numCol();
00505         
00506         for(unsigned int ii=0; ii<mRow; ii++)
00507             for(unsigned int jj=0; jj<mCol; jj++) {
00508                 *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
00509             }
00510         return (*this);
00511     }
00512 }
00513 
00514 template<class DataType>
00515 StMatrix<DataType>& StMatrix<DataType>::operator=(const StMatrix<DataType>& m1)
00516 {   
00517     if ((void *)&m1 == (void *)this) {
00518         return *this;
00519     }
00520     else {
00521         delete [] mElement;
00522         mSize    = m1.numRow()*m1.numCol();
00523         mElement = new DataType[mSize];
00524 
00525         mRow = m1.numRow();
00526         mCol = m1.numCol();
00527         
00528         for(unsigned int ii=0; ii<mRow; ii++)
00529             for(unsigned int jj=0; jj<mCol; jj++) {
00530                 *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
00531             }
00532         return (*this);
00533     }
00534 }
00535 #else // ST_NO_MEMBER_TEMPLATES
00536 template<class DataType>
00537 StMatrix<DataType>& StMatrix<DataType>::operator=(const StMatrix<float>& m1)
00538 {   
00539     if ((void *)&m1 == (void *)this) {
00540         return *this;
00541     }
00542     else {
00543         delete [] mElement;
00544         mSize    = m1.numRow()*m1.numCol();
00545         mElement = new DataType[mSize];
00546 
00547         mRow = m1.numRow();
00548         mCol = m1.numCol();
00549         
00550         for(unsigned int ii=0; ii<mRow; ii++)
00551             for(unsigned int jj=0; jj<mCol; jj++) {
00552                 *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
00553             }
00554         return (*this);
00555     }
00556 }
00557 template<class DataType>
00558 StMatrix<DataType>& StMatrix<DataType>::operator=(const StMatrix<double>& m1)
00559 {   
00560     if ((void *)&m1 == (void *)this) {
00561         return *this;
00562     }
00563     else {
00564         delete [] mElement;
00565         mSize    = m1.numRow()*m1.numCol();
00566         mElement = new DataType[mSize];
00567 
00568         mRow = m1.numRow();
00569         mCol = m1.numCol();
00570         
00571         for(unsigned int ii=0; ii<mRow; ii++)
00572             for(unsigned int jj=0; jj<mCol; jj++) {
00573                 *(mElement+(ii)*mCol+jj) = m1(ii+1,jj+1);
00574             }
00575         return (*this);
00576     }
00577 }
00578 
00579 #endif
00580 
00581 // Destructor
00582 template<class DataType>
00583 StMatrix<DataType>::~StMatrix() {
00584     delete [] mElement;
00585 }
00586 
00587 //
00588 // access functions
00589 //
00590 template<class DataType>
00591 unsigned int StMatrix<DataType>::numRow() const { return mRow;}
00592 
00593 template<class DataType>
00594 unsigned int StMatrix<DataType>::numCol() const  { return mCol;}
00595 
00596 template<class DataType>
00597 unsigned int StMatrix<DataType>::numSize() const { return mSize;}
00598 
00599 // backward compatibility
00600 template<class DataType>
00601 unsigned int StMatrix<DataType>::num_row() const { return mRow;}
00602 
00603 template<class DataType>
00604 unsigned int StMatrix<DataType>::num_col() const  { return mCol;}
00605 
00606 template<class DataType>
00607 unsigned int StMatrix<DataType>::num_size() const { return mSize;}
00608 
00609 //
00610 // Element Access
00611 //
00612 template<class DataType>
00613 const DataType& StMatrix<DataType>::operator()(size_t row, size_t col) const 
00614 {
00615 #ifdef MATRIX_BOUND_CHECK
00616     if(row<1 || row>numRow() || col<1 || col>numCol()) {
00617 #ifndef ST_NO_EXCEPTIONS
00618         throw out_of_range("StMatrix<DataType>::operator(): const Bad Index");
00619 #else
00620         cerr << "StMatrix<DataType>::operator(): const Bad Index" << endl;
00621 #endif
00622     }
00623 #endif
00624     return *(mElement+(row-1)*mCol+col-1);
00625 }
00626 
00627 template<class DataType>
00628 DataType& StMatrix<DataType>::operator()(size_t row, size_t col)
00629 {
00630 #ifdef MATRIX_BOUND_CHECK
00631     if(row<1 || row>mRow || col<1 || col>mCol) {
00632 #ifndef ST_NO_EXCEPTIONS
00633         throw out_of_range("StMatrix<DataType>::operator(): Bad Index");
00634 #else
00635         cerr << "StMatrix<DataType>::operator(): Bad Index" << endl;
00636 #endif
00637     }
00638 #endif
00639     return *(mElement+(row-1)*mCol+col-1);
00640 }
00641 
00642 /* -----------------------------------------------------------------------
00643    This section contains the assignment and inplace operators =,+=,-=,*=,/=.
00644    ----------------------------------------------------------------------- */
00645 
00646 template<class DataType>
00647 StMatrix<DataType>& StMatrix<DataType>::operator*=(double fact)
00648 {
00649     for(unsigned int ii=0; ii<mRow; ii++)
00650         for(unsigned int jj=0; jj<mCol; jj++)
00651             *(mElement+(ii)*mCol+jj) *= fact;
00652     
00653     return (*this);
00654 }
00655 
00656 template<class DataType>
00657 StMatrix<DataType> & StMatrix<DataType>::operator/=(double fact)
00658 {
00659     if(fact == 0) {
00660 #ifndef ST_NO_EXCEPTIONS
00661         throw out_of_range("StMatrix<T>::operator/=(): Cannot divide by zero!");
00662 #else
00663         cerr << "StMatrix<T>::operator/=(): Cannot divide by zero!" << endl;
00664 #endif
00665     }
00666     for(unsigned int ii=0; ii<mCol; ii++)
00667         for(unsigned int jj=0; jj<mRow; jj++)
00668             *(mElement+(ii)*mCol+jj) /= fact;
00669     
00670     return (*this);
00671 }
00672 
00673 #if !defined(ST_NO_MEMBER_TEMPLATES) && !defined(__CINT__)
00674 template<class DataType>
00675 template<class X>
00676 StMatrix<DataType>& StMatrix<DataType>::operator+=(const StMatrix<X>& m2)
00677 {
00678     if(mRow == m2.numRow() && mCol == m2.numCol()) {
00679         for(unsigned int ii=0; ii<mRow; ii++)
00680             for(unsigned int jj=0; jj<mCol; jj++)
00681                 *(mElement+(ii)*mCol+jj) += m2(ii+1,jj+1);
00682     }
00683     else {
00684 #ifndef ST_NO_EXCEPTIONS
00685         throw out_of_range("StMatrix<T>::operator+=(): Matrices are not same size!");
00686 #else
00687         cerr << "StMatrix<T>::operator+=(): Matrices are not same size!" << endl;
00688 #endif
00689     }   
00690     return (*this);
00691 }
00692 
00693 template<class DataType>
00694 template<class X>
00695 StMatrix<DataType>& StMatrix<DataType>::operator-=(const StMatrix<X>& m2)
00696 {
00697     if(mRow == m2.numRow() && mCol == m2.numCol()) {
00698         for(unsigned int ii=0; ii<mRow; ii++)
00699             for(unsigned int jj=0; jj<mCol; jj++)
00700                 *(mElement+(ii)*mCol+jj) -= m2(ii+1,jj+1);
00701     }
00702     else {
00703 #ifndef ST_NO_EXCEPTIONS
00704         throw out_of_range("StMatrix<T>::operator-=(): Matrices are not same size!");
00705 #else
00706         cerr << "StMatrix<T>::operator-=(): Matrices are not same size!" << endl;
00707 #endif
00708     }
00709     return (*this);
00710 }
00711 
00712 template<class DataType>
00713 template<class X>
00714 StMatrix<DataType> StMatrix<DataType>::dot(const StMatrix<X>& m2)
00715 {
00716     if(mCol == m2.numRow() ) {
00717         StMatrix<DataType> mret(mRow, m2.numCol(), 0);  
00718         for(unsigned int i=0; i<mRow; i++)
00719             for(unsigned int j=0; j<m2.numCol(); j++) {
00720                 for(unsigned int kk=0; kk<mCol; kk++)
00721                     mret(i+1, j+1) += (*(mElement+(i)*mCol+kk))*m2(kk+1,j+1);
00722             }
00723         return mret;
00724     }
00725     else {
00726 #ifndef ST_NO_EXCEPTIONS
00727         throw out_of_range("StMatrix<T>::dot(): incompatible matrix sizes");
00728 #else
00729         cerr << "StMatrix<T>::dot(): incompatible matrix sizes" << endl;
00730 #endif
00731         return StMatrix();
00732     }
00733 }
00734 #else  // ST_NO_MEMBER_TEMPLATES
00735 // operator+=
00736 template<class DataType>
00737 StMatrix<DataType>& StMatrix<DataType>::operator+=(const StMatrix<float>& m2)
00738 {
00739     if(mRow == m2.numRow() && mCol == m2.numCol()) {
00740         for(unsigned int ii=0; ii<mRow; ii++)
00741             for(unsigned int jj=0; jj<mCol; jj++)
00742                 *(mElement+(ii)*mCol+jj) += m2(ii+1,jj+1);
00743     }
00744     else {
00745 #ifndef ST_NO_EXCEPTIONS
00746         throw out_of_range("StMatrix<float>::operator+=(): Matrices are not same size!");
00747 #else
00748         cerr << "StMatrix<float>::operator+=(): Matrices are not same size!" << endl;
00749 #endif
00750     }   
00751     return (*this);
00752 }
00753 
00754 template<class DataType>
00755 StMatrix<DataType>& StMatrix<DataType>::operator+=(const StMatrix<double>& m2)
00756 {
00757     if(mRow == m2.numRow() && mCol == m2.numCol()) {
00758         for(unsigned int ii=0; ii<mRow; ii++)
00759             for(unsigned int jj=0; jj<mCol; jj++)
00760                 *(mElement+(ii)*mCol+jj) += m2(ii+1,jj+1);
00761     }
00762     else {
00763 #ifndef ST_NO_EXCEPTIONS
00764         throw out_of_range("StMatrix<double>::operator+=(): Matrices are not same size!");
00765 #else
00766         cerr << "StMatrix<double>::operator+= Matrices are not same size!" << endl;
00767 #endif
00768     }   
00769     return (*this);
00770 }
00771 
00772 // operator -=
00773 template<class DataType>
00774 StMatrix<DataType>& StMatrix<DataType>::operator-=(const StMatrix<float>& m2)
00775 {
00776     if(mRow == m2.numRow() && mCol == m2.numCol()) {
00777         for(unsigned int ii=0; ii<mRow; ii++)
00778             for(unsigned int jj=0; jj<mCol; jj++)
00779                 *(mElement+(ii)*mCol+jj) -= m2(ii+1,jj+1);
00780     }
00781     else {
00782 #ifndef ST_NO_EXCEPTIONS
00783         throw out_of_range("StMatrix<float>::operator-=(): Matrices are not same size!");
00784 #else
00785         cerr << "StMatrix<float>::operator-=(): Matrices are not same size!" << endl;
00786 #endif
00787     }
00788     return (*this);
00789 }
00790 
00791 template<class DataType>
00792 StMatrix<DataType>& StMatrix<DataType>::operator-=(const StMatrix<double>& m2)
00793 {
00794     if(mRow == m2.numRow() && mCol == m2.numCol()) {
00795         for(unsigned int ii=0; ii<mRow; ii++)
00796             for(unsigned int jj=0; jj<mCol; jj++)
00797                 *(mElement+(ii)*mCol+jj) -= m2(ii+1,jj+1);
00798     }
00799     else {
00800 #ifndef ST_NO_EXCEPTIONS
00801         throw out_of_range("StMatrix<double>::operator-=(): Matrices are not same size!");
00802 #else
00803         cerr << "StMatrix<double>::operator-=(): Matrices are not same size!" << endl;
00804 #endif
00805     }
00806     return (*this);
00807 }
00808 
00809 // operator::dot
00810 template<class DataType>
00811 StMatrix<DataType> StMatrix<DataType>::dot(const StMatrix<float>& m2)
00812 {
00813     if(mCol == m2.numRow() ) {
00814         StMatrix<DataType> mret(mRow, m2.numCol(), 0);
00815         
00816         for(unsigned int i=0; i<mRow; i++)
00817             for(unsigned int j=0; j<m2.numCol(); j++) {
00818                 for(unsigned int kk=0; kk<mCol; kk++)
00819                     mret(i+1, j+1) += (*(mElement+(i)*mCol+kk))*m2(kk+1,j+1);
00820             }
00821         return mret;
00822     }
00823     else {
00824 #ifndef ST_NO_EXCEPTIONS
00825         throw out_of_range("StMatrix<float>::dot(): Incompatible matrix sizes");
00826 #else
00827         cerr << "StMatrix<float>::dot(): Incompatible matrix sizes" << endl;
00828 #endif
00829         return StMatrix();
00830     }
00831 }
00832 
00833 template<class DataType>
00834 StMatrix<DataType> StMatrix<DataType>::dot(const StMatrix<double>& m2)
00835 {
00836     if(mCol == m2.numRow() ) {
00837         StMatrix<DataType> mret(mRow, m2.numCol(), 0);
00838         
00839         for(unsigned int i=0; i<mRow; i++)
00840             for(unsigned int j=0; j<m2.numCol(); j++) {
00841                 for(unsigned int kk=0; kk<mCol; kk++)
00842                     mret(i+1, j+1) += (*(mElement+(i)*mCol+kk))*m2(kk+1,j+1);
00843             }
00844         return mret;
00845     }
00846     else {
00847 #ifndef ST_NO_EXCEPTIONS
00848         throw out_of_range("StMatrix<double>::dot(): Incompatible matrix sizes");
00849 #else
00850         cerr << "StMatrix<double>::dot(): Incompatible matrix sizes" << endl;
00851 #endif
00852         return StMatrix();
00853     }
00854 }
00855 
00857 #endif // ST_NO_MEMBER_TEMPLATES
00858 
00859 template<class DataType>
00860 StMatrix<DataType> StMatrix<DataType>::operator+ () const 
00861 {
00862     return *this;
00863 }
00864 
00865 template<class DataType>
00866 StMatrix<DataType> StMatrix<DataType>::operator- () const 
00867 {
00868     // this can be streamlined
00869     StMatrix<DataType> mret(*this);
00870     return mret*=-1;
00871 }
00872 
00873 template<class DataType>
00874 bool StMatrix<DataType>::operator== (const StMatrix<DataType>& m1) const
00875 {
00876     if (mCol == m1.numCol() && mRow == m1.numRow()) {
00877         for(unsigned int ii=0; ii<mRow; ii++)
00878             for(unsigned int jj=0; jj<mCol; jj++) {
00879                 if(*(mElement+(ii)*mCol+jj) != m1(ii+1,jj+1))
00880                     return false;
00881             }
00882         return true;
00883     }
00884     else
00885         return false;  // not even right size
00886 }
00887 
00888 template<class DataType>
00889 bool StMatrix<DataType>::operator!= (const StMatrix<DataType>& m1) const
00890 {
00891     return !(*this == m1);
00892 }
00893 
00894 //Apply a function to all elements
00895 template<class DataType>
00896 StMatrix<DataType> StMatrix<DataType>::apply(DataType (*f)(DataType, size_t, size_t)) const
00897 {
00898     StMatrix<DataType> mret(mRow, mCol);
00899     DataType *a = mElement;
00900     for(unsigned int ir=1; ir<=mRow; ir++) {
00901         for(unsigned int ic=1; ic<=mCol; ic++) {
00902             mret(ir,ic) = (*f)(*(a++), ir, ic);
00903         }
00904     }
00905     return mret;
00906 }
00907 
00908 template<class DataType>
00909 StMatrix<DataType> StMatrix<DataType>::T() const
00910 {
00911     StMatrix<DataType> mret(mCol,mRow);
00912     register DataType *pl = mElement + mSize;
00913     register DataType *pme = mElement;
00914     register DataType *pt = mret.mElement;
00915     register DataType *ptl = mret.mElement + mSize;
00916     for (; pme < pl; pme++, pt+=mRow)
00917         {
00918             if (pt >= ptl)
00919                 pt -= (mSize-1);
00920             (*pt) = (*pme);
00921         }
00922     return mret;
00923 }
00924 template<class DataType>
00925 StMatrix<DataType> StMatrix<DataType>::transpose() const
00926 {
00927     return this->T();
00928 }
00929 
00930 //
00931 // Sub matrix
00932 //
00933 template<class DataType>
00934 StMatrix<DataType> StMatrix<DataType>::sub(size_t min_row, size_t max_row,
00935                                                  size_t min_col,size_t max_col) const
00936 {
00937     StMatrix<DataType> mret(max_row-min_row+1,max_col-min_col+1);
00938     if(max_row > mRow || max_col > mCol) {
00939 #ifndef ST_NO_EXCEPTIONS
00940         throw out_of_range("StMatrix<DataType>::sub(): Index out of range");
00941 #else
00942         cerr << "StMatrix<DataType>::sub(): Index out of range" << endl;
00943 #endif
00944     }
00945     int nrows = mret.numRow();
00946     int ncols = mret.numCol();
00947     for(int ii=0; ii<nrows; ii++)
00948         for(int jj=0; jj<ncols; jj++) {
00949             mret(ii+1, jj+1) =
00950                 *(mElement+(min_row+ii-1)*mCol+(min_col+jj-1));
00951         }
00952     return mret;
00953 }
00954 
00955 template<class DataType>
00956 void StMatrix<DataType>::sub(size_t row,size_t col,const StMatrix<DataType> &m1)
00957 {
00958     if((row <1)                   ||
00959        (row+m1.numRow()-1 > mRow) ||
00960        (col <1)                   ||
00961        (col+m1.numCol()-1 > mCol)) {
00962 #ifndef ST_NO_EXCEPTIONS
00963         throw out_of_range("StMatrix<DataType>::sub(): Index out of range");
00964 #else
00965         cerr << "StMatrix<DataType>::sub(): Index out of range" << endl;
00966 #endif
00967     }
00968     DataType *a = m1.mElement;
00969     unsigned int nc = numCol();
00970     DataType *b1 = mElement + (row - 1) * nc + col - 1;
00971     
00972     for(unsigned int irow=1; irow<=m1.numRow(); irow++) {
00973         DataType *brc = b1;
00974         for(unsigned int icol=1; icol<=m1.numCol(); icol++) {
00975             *(brc++) = *(a++);
00976         }
00977         b1 += nc;
00978     }
00979 }
00980 
00981 // Contents of the matrix are replaced by the inverse--less overhead then inverse().
00982 template<class DataType>
00983 StMatrix<DataType> StMatrix<DataType>::inverse(size_t &ierr) const
00984 {
00985     StMatrix<DataType> tmp(*this);
00986     tmp.invert(ierr);
00987     return tmp;
00988 }
00989 
00990 template<class DataType>
00991 void StMatrix<DataType>::invert(size_t &ierr) {
00992     if(mCol != mRow) {
00993 #ifndef ST_NO_EXCEPTIONS
00994         throw domain_error("StMatrix<DataType>::invert(): not a NxN matrix");
00995 #else
00996         cerr << "StMatrix<DataType>::invert(): not a NxN matrix" << endl;
00997 #endif
00998     }
00999     static unsigned int max_array = 20;
01000     static int *ir = new int [max_array+1];
01001     
01002     if (mCol > max_array) {
01003         delete [] ir;
01004         max_array = mRow;
01005         ir = new int [max_array+1];
01006     }
01007     DataType t1, t2, t3;
01008     DataType det, temp, s;
01009     int ifail;
01010     switch(mRow) {
01011     case 3:
01012         {
01013             DataType c11,c12,c13,c21,c22,c23,c31,c32,c33;
01014             ifail = 0;
01015             c11 = (*(mElement+4)) * (*(mElement+8)) - (*(mElement+5)) * (*(mElement+7));
01016             c12 = (*(mElement+5)) * (*(mElement+6)) - (*(mElement+3)) * (*(mElement+8));
01017             c13 = (*(mElement+3)) * (*(mElement+7)) - (*(mElement+4)) * (*(mElement+6));
01018             c21 = (*(mElement+7)) * (*(mElement+2)) - (*(mElement+8)) * (*(mElement+1));
01019             c22 = (*(mElement+8)) * (*mElement)     - (*(mElement+6)) * (*(mElement+2));
01020             c23 = (*(mElement+6)) * (*(mElement+1)) - (*(mElement+7)) * (*mElement);
01021             c31 = (*(mElement+1)) * (*(mElement+5)) - (*(mElement+2)) * (*(mElement+4));
01022             c32 = (*(mElement+2)) * (*(mElement+3)) - (*mElement) * (*(mElement+5));
01023             c33 = (*mElement)     * (*(mElement+4)) - (*(mElement+1)) * (*(mElement+3));
01024             t1  = fabs(*mElement);
01025             t2  = fabs(*(mElement+3));
01026             t3  = fabs(*(mElement+6));
01027             if (t1 >= t2) {
01028                 if (t3 >= t1) {
01029                     temp = *(mElement+6);
01030                     det = c23*c12-c22*c13;
01031                 } else {
01032                     temp = *mElement;
01033                     det = c22*c33-c23*c32;
01034                 }
01035             } else if (t3 >= t2) {
01036                 temp = *(mElement+6);
01037                 det = c23*c12-c22*c13;
01038             } else {
01039                 temp = *(mElement+3);
01040                 det = c13*c32-c12*c33;
01041             }
01042             if (det==0) {
01043                 ierr = 1;
01044                 return;
01045             }
01046             {
01047                 DataType s = temp/det;
01048                 DataType *tmp = mElement;
01049                 *(tmp++) = s*c11;
01050                 *(tmp++) = s*c21;
01051                 *(tmp++) = s*c31;
01052                 *(tmp++) = s*c12;
01053                 *(tmp++) = s*c22;
01054                 *(tmp++) = s*c32;
01055                 *(tmp++) = s*c13;
01056                 *(tmp++) = s*c23;
01057                 *(tmp) = s*c33;
01058             }
01059         }
01060         break;
01061     case 2:
01062         ifail = 0;
01063         det = (*mElement)*(*(mElement+3)) - (*(mElement+1))*(*(mElement+2));
01064         if (det==0) {
01065             ierr = 1;
01066             return;
01067         }
01068         s = 1.0/det;
01069         temp = s*(*(mElement+3));
01070         *(mElement+1) *= -s;
01071         *(mElement+2) *= -s;
01072         *(mElement+3) = s*(*mElement);
01073         *mElement = temp;
01074         break;
01075     case 1:
01076         ifail = 0;
01077         if ((*mElement)==0) {
01078             ierr = 0;
01079             return;
01080         }
01081         *mElement = 1.0/(*mElement);
01082         break;
01083     default:
01084         ifail = dfact(det, ir);
01085         if(ifail) {
01086             ierr = 1;
01087             return;
01088         }
01089         dfinv(ir);
01090         break;
01091     }
01092     ierr = 0;
01093     return;
01094 }
01095 
01096 template<class DataType>
01097 DataType StMatrix<DataType>::determinant() const {
01098     static unsigned int max_array = 20;
01099     static int *ir = new int [max_array+1];
01100     if(mCol != mRow) {
01101 #ifndef ST_NO_EXCEPTIONS
01102                 throw out_of_range("StMatrix<DataType>::determinant(): not a NxN matrix");
01103 #else
01104                 cerr << "StMatrix<DataType>::determinant(): not a NxN matrix" << endl;
01105 #endif
01106     }
01107     if (mCol > max_array) {
01108         delete [] ir;
01109         max_array = mRow;
01110         ir = new int [max_array+1];
01111     }
01112     DataType det;
01113     StMatrix<DataType> mt(*this);
01114     int i = mt.dfact(det, ir);
01115     if(i==0) return det;
01116     return 0;
01117 }
01118 
01119 //
01120 // Access operators
01121 //
01122 template<class DataType>
01123 StMatrix<DataType>::StMatrixRow::StMatrixRow(StMatrix<DataType>& a, size_t r) 
01124     : _a(a) {
01125         _r = r;
01126 }
01127 
01128 template<class DataType>
01129 DataType& StMatrix<DataType>::StMatrixRow::operator[](size_t c) {
01130 #ifdef MATRIX_BOUND_CHECK
01131     if (_r<0 || _r>=_a.numRow() || c<0 || c>=_a.numCol()) {
01132 #ifndef ST_NO_EXCEPTIONS
01133         throw out_of_range("StMatrix<DataType>::operator[]: index out of range");
01134 #else
01135         cerr << "StMatrix<DataType>::operator[]: index out of range" << endl;
01136 #endif
01137     }
01138 #endif
01139     return *(_a.mElement+_r*_a.mCol+c);
01140 }
01141 
01142 template<class DataType>
01143 StMatrix<DataType>::StMatrixRowConst::StMatrixRowConst(const StMatrix<DataType>& a, size_t r) 
01144     : _a(a) {
01145     _r = r;
01146 }
01147 
01148 template<class DataType>
01149 const DataType& StMatrix<DataType>::StMatrixRowConst::operator[](size_t c) const
01150 {
01151 #ifdef MATRIX_BOUND_CHECK
01152     if (_r<0 || _r>=_a.numRow() || c<0 || c>=_a.numCol()) {
01153 #ifndef ST_NO_EXCEPTIONS
01154         throw out_of_range("StMatrix<DataType>::operator[]: const index out of range");
01155 #else
01156         cerr << "StMatrix<DataType>::operator[]: const index out of range" << endl;
01157 #endif
01158     }
01159 #endif
01160     return *(_a.mElement+_r*_a.mCol+c);
01161 }
01162 
01163 template<class DataType> 
01164 typename StMatrix<DataType>::StMatrixRow StMatrix<DataType>::operator[] (size_t r)
01165 {
01166     StMatrixRow b(*this,r);
01167     return b;
01168 }
01169 
01170 template<class DataType>
01171 const typename StMatrix<DataType>::StMatrixRowConst StMatrix<DataType>::operator[] (size_t r) const
01172 {
01173     StMatrixRowConst b(*this,r);
01174     return b;
01175 }
01176 
01177 template<class DataType>
01178 void StMatrix<DataType>::swap(unsigned int &i, unsigned int &j)
01179 {
01180     unsigned int tmp=i;
01181     i=j;
01182     j=tmp;
01183 }
01184 
01185 template<class DataType>
01186 void StMatrix<DataType>::swap(DataType *&i, DataType *&j)
01187 {
01188     DataType *tmp=i;
01189     i=j;
01190     j=tmp;
01191 }
01192 
01193 // This function swaps two Matrices doing a full copy
01194 template<class DataType>
01195 void swap(StMatrix<DataType>& m1, StMatrix<DataType>& m2) {
01196     StMatrix<DataType> tmp(m1);
01197     m1 = m2;
01198     m2 = tmp;
01199 }
01200 
01201 // This function swaps two Matrices without doing a full copy.
01202 // The version below had some problems with the new GCC compiler
01203 // (friendship with StMatrix caused trouble) and was replaced
01204 // by the version above. tu
01205 // template<class DataType>
01206 // void swap(StMatrix<DataType>& m1, StMatrix<DataType>& m2) {
01207 //     StMatrix<DataType>::swap(m1.mElement, m2.mElement);
01208 //     StMatrix<DataType>::swap(m1.mRow,     m2.mRow);
01209 //     StMatrix<DataType>::swap(m1.mCol,     m2.mCol);
01210 //     StMatrix<DataType>::swap(m1.mSize,    m2.mSize);
01211 // }
01212 
01213 
01214 template<class DataType>
01215 int StMatrix<DataType>::dfact(DataType &det, int *ir)
01216 {
01217     if (mCol!=mRow) {
01218 #ifndef ST_NO_EXCEPTIONS
01219         throw domain_error("StMatrix<DataType>::dfact(): Matrix not NxN");
01220 #else
01221         cerr << "StMatrix<DataType>::dfact(): Matrix not NxN" << endl;
01222 #endif
01223     }
01224  
01225     int ifail, jfail;
01226     register int n = mCol;
01227  
01228     DataType tf;
01229     DataType g1 = 1.0e-19;
01230     DataType g2 = 1.0e19;
01231  
01232     DataType p, q, t;
01233     DataType s11, s12;
01234 
01235     
01236     DataType epsilon;
01237     if (sizeof(DataType) == sizeof(double))
01238         epsilon = 8*DBL_EPSILON;
01239     else
01240         epsilon = 8*FLT_EPSILON;
01241     // could be set to zero (like it was before)
01242     // but then the algorithm often doesn't detect
01243     // that a matrix is singular
01244     
01245     int normal = 0, imposs = -1;
01246     int jrange = 0, jover = 1, junder = -1;
01247     ifail = normal;
01248     jfail = jrange;
01249     int nxch = 0;
01250     det = 1.0;
01251     DataType *mj = mElement;
01252     DataType *mjj = mj;
01253     for (int j=1;j<=n;j++) {
01254         int k = j;
01255         p = (fabs(*mjj));
01256         if (j!=n) {
01257             DataType *mij = mj + n + j - 1;
01258             for (int i=j+1;i<=n;i++) {
01259                 q = (fabs(*(mij)));
01260                 if (q > p) {
01261                     k = i;
01262                     p = q;
01263                 }
01264                 mij += n;
01265             }
01266             if (k==j) {
01267                 if (p <= epsilon) {
01268                     det = 0;
01269                     ifail = imposs;
01270                     jfail = jrange;
01271                     return ifail;
01272                 }
01273                 det = -det; // in this case the sign of the determinant
01274                 // must not change. So I change it twice.
01275             }
01276             DataType *mjl = mj;
01277             DataType *mkl = mElement + (k-1)*n;
01278             for (int l=1;l<=n;l++) {
01279                 tf = *mjl;
01280                 *(mjl++) = *mkl;
01281                 *(mkl++) = tf;
01282             }
01283             nxch = nxch + 1;  // this makes the determinant change its sign
01284             ir[nxch] = (((j)<<12)+(k));
01285         } else {
01286             if (p <= epsilon) {
01287                 det = 0.0;
01288                 ifail = imposs;
01289                 jfail = jrange;
01290                 return ifail;
01291             }
01292         }
01293         det *= *mjj;
01294         *mjj = 1.0 / *mjj;
01295         t = (fabs(det));
01296         if (t < g1) {
01297             det = 0.0;
01298             if (jfail == jrange) jfail = junder;
01299         }
01300         else if (t > g2) {
01301             det = 1.0;
01302             if (jfail==jrange) jfail = jover;
01303         }
01304         if (j!=n) {
01305             DataType *mk = mj + n;
01306             DataType *mkjp = mk + j;
01307             DataType *mjk = mj + j;
01308             for (k=j+1;k<=n;k++) {
01309                 s11 = - (*mjk);
01310                 s12 = - (*mkjp);
01311                 if (j!=1) {
01312                     DataType *mik = mElement + k - 1;
01313                     DataType *mijp = mElement + j;
01314                     DataType *mki = mk;
01315                     DataType *mji = mj;
01316                     for (int i=1;i<j;i++) {
01317                         s11 += (*mik) * (*(mji++));
01318                         s12 += (*mijp) * (*(mki++));
01319                         mik += n;
01320                         mijp += n;
01321                     }
01322                 }
01323                 *(mjk++) = -s11 * (*mjj);
01324                 *(mkjp) = -(((*(mjj+1)))*((*(mkjp-1)))+(s12));
01325                 mk += n;
01326                 mkjp += n;
01327             }
01328         }
01329         mj += n;
01330         mjj += (n+1);
01331     }
01332     if (nxch%2==1) det = -det;
01333     if (jfail !=jrange) det = 0.0;
01334     ir[n] = nxch;
01335     return 0;
01336 }
01337 
01338 template<class DataType>
01339 int StMatrix<DataType>::dfinv(int *ir)
01340 {
01341     if (mCol != mRow) {
01342 #ifndef ST_NO_EXCEPTIONS
01343         throw domain_error("StMatrix<DataType>::dfinv(): Matrix not NxN");
01344 #else
01345         cerr << "StMatrix<DataType>::dfinv(): Matrix not NxN" << endl;
01346 #endif
01347     }
01348     
01349     register int n = mCol;
01350     if (n==1) return 0;
01351     
01352     DataType s31, s32;
01353     register DataType s33, s34;
01354     
01355     DataType *m11 = mElement;
01356     DataType *m12 = m11 + 1;
01357     DataType *m21 = m11 + n;
01358     DataType *m22 = m12 + n;
01359     *m21 = -(*m22) * (*m11) * (*m21);
01360     *m12 = -(*m12);
01361     if (n>2) {
01362         DataType *mi = mElement + 2 * n;
01363         DataType *mii= mElement + 2 * n + 2;
01364         DataType *mimim = mElement + n + 1;
01365         for (int i=3;i<=n;i++) {
01366             int im2 = i - 2;
01367             DataType *mj = mElement;
01368             DataType *mji = mj + i - 1;
01369             DataType *mij = mi;
01370             for (int j=1;j<=im2;j++) {
01371                 s31 = 0.0;
01372                 s32 = *mji;
01373                 DataType *mkj = mj + j - 1;
01374                 DataType *mik = mi + j - 1;
01375                 DataType *mjkp = mj + j;
01376                 DataType *mkpi = mj + n + i - 1;
01377                 for (int k=j;k<=im2;k++) {
01378                     s31 += (*mkj) * (*(mik++));
01379                     s32 += (*(mjkp++)) * (*mkpi);
01380                     mkj += n;
01381                     mkpi += n;
01382                 }
01383                 *mij = -(*mii) * (((*(mij-n)))*( (*(mii-1)))+(s31));
01384                 *mji = -s32;
01385                 mj += n;
01386                 mji += n;
01387                 mij++;
01388             }
01389             *(mii-1) = -(*mii) * (*mimim) * (*(mii-1));
01390             *(mimim+1) = -(*(mimim+1));
01391             mi += n;
01392             mimim += (n+1);
01393             mii += (n+1);
01394         }
01395     }
01396     DataType *mi = mElement;
01397     DataType *mii = mElement;
01398     for (int i=1;i<n;i++) {
01399         int ni = n - i;
01400         DataType *mij = mi;
01401         int j;
01402         for (j=1; j<=i;j++) {
01403             s33 = *mij;
01404             register DataType *mikj = mi + n + j - 1;
01405             register DataType *miik = mii + 1;
01406             DataType *min_end = mi + n;
01407             for (;miik<min_end;) {
01408                 s33 += (*mikj) * (*(miik++));
01409                 mikj += n;
01410             }
01411             *(mij++) = s33;
01412         }
01413         for (j=1;j<=ni;j++) {
01414             s34 = 0.0;
01415             DataType *miik = mii + j;
01416             DataType *mikij = mii + j * n + j;
01417             for (int k=j;k<=ni;k++) {
01418                 s34 += *mikij * (*(miik++));
01419                 mikij += n;
01420             }
01421             *(mii+j) = s34;
01422         }
01423         mi += n;
01424         mii += (n+1);
01425     }
01426     int nxch = ir[n];
01427     if (nxch==0) return 0;
01428     for (int mm=1;mm<=nxch;mm++) {
01429         int k = nxch - mm + 1;
01430         int ij = ir[k];
01431         int i = ij >> 12;
01432         int j = ij%4096;
01433         DataType *mki = mElement + i - 1;
01434         DataType *mkj =mElement + j - 1;
01435         for (k=1; k<=n;k++) {
01436             DataType ti = *mki;
01437             *mki = *mkj;
01438             *mkj = ti;
01439             mki += n;
01440             mkj += n;
01441         }
01442     }
01443     return 0;
01444 }
01445  
01446 
01447 
01448 template<class DataType>
01449 StMatrix<DataType> operator/(const StMatrix<DataType>& m1, double fact)
01450 {
01451     StMatrix<DataType> mret(m1);
01452     mret /= fact;
01453     return mret;
01454 }
01455 
01456 template<class DataType>
01457 StMatrix<DataType> operator*(const StMatrix<DataType>& m1, double fact)
01458 {
01459     StMatrix<DataType> mret(m1);
01460     mret *= fact;
01461     return mret;
01462 }
01463 
01464 template<class DataType>
01465 StMatrix<DataType> operator*(double fact,const StMatrix<DataType>& m1)
01466 {
01467     StMatrix<DataType> mret(m1);
01468     mret *= fact;
01469     return mret;
01470 }
01471 
01472 // Private Members
01473 // Direct sum of two matricies
01474 template<class DataType>
01475 StMatrix<DataType> dsum(const StMatrix<DataType> &m1, const StMatrix<DataType> &m2)
01476 {
01477     StMatrix<DataType> mret(m1.numRow() + m2.numRow(), m1.numCol() + m2.numCol());
01478     mret.sub(1,1,m1);
01479     mret.sub(m1.numRow()+1, m1.numCol()+1, m2);
01480     return mret;
01481 }
01482 
01483 #endif /* ! __CINT__ */
01484 #ifdef __CINT__
01485 template<> StMatrix<double> operator*(const StMatrix<double>& m1,const StMatrix<double>& m2);
01486 template<> StMatrix<double> operator*(const StMatrix<double>& m1,const StMatrix<float>&  m2);
01487 template<> StMatrix<double> operator*(const StMatrix<float>&  m1,const StMatrix<double>& m2);
01488 template<> StMatrix<float>  operator*(const StMatrix<float>&  m1,const StMatrix<float>&  m2);
01489 
01490 template<> vector<double> operator*(const StMatrix<double>& m1, const vector<double>& v);
01491 template<> vector<double> operator*(const vector<double>&    v, const StMatrix<double>& m1);
01492 template<> vector<double> operator*(const StMatrix<double>& m1, const vector<float>& v);
01493 template<> vector<double> operator*(const vector<float>&     v, const StMatrix<double>& m1);
01494 template<> vector<double> operator*(const StMatrix<float>&  m1, const vector<double>& v);
01495 template<> vector<double> operator*(const vector<double>&    v, const StMatrix<float>& m1);
01496 template<> vector<float>  operator*(const StMatrix<float>&  m1, const vector<float>& v);
01497 template<> vector<float>  operator*(const vector<float>&     v, const StMatrix<float>& m1);
01498 
01499 template<> StThreeVector<double> operator*(const StMatrix<double>& m1, const StThreeVector<double>& v);
01500 template<> StThreeVector<double> operator*(const StThreeVector<double>&    v, const StMatrix<double>& m1);
01501 template<> StThreeVector<double> operator*(const StMatrix<double>& m1, const StThreeVector<float>& v);
01502 template<> StThreeVector<double> operator*(const StThreeVector<float>&     v, const StMatrix<double>& m1);
01503 template<> StThreeVector<double> operator*(const StMatrix<float>&  m1, const StThreeVector<double>& v);
01504 template<> StThreeVector<double> operator*(const StThreeVector<double>&    v, const StMatrix<float>& m1);
01505 template<> StThreeVector<float>  operator*(const StMatrix<float>&  m1, const StThreeVector<float>& v);
01506 template<> StThreeVector<float>  operator*(const StThreeVector<float>&     v, const StMatrix<float>& m1);
01507 
01508 template<> StLorentzVector<double> operator*(const StMatrix<double>& m1, const StLorentzVector<double>& v);
01509 template<> StLorentzVector<double> operator*(const StLorentzVector<double>&    v, const StMatrix<double>& m1);
01510 template<> StLorentzVector<double> operator*(const StMatrix<double>& m1, const StLorentzVector<float>& v);
01511 template<> StLorentzVector<double> operator*(const StLorentzVector<float>&     v, const StMatrix<double>& m1);
01512 template<> StLorentzVector<double> operator*(const StMatrix<float>&  m1, const StLorentzVector<double>& v);
01513 template<> StLorentzVector<double> operator*(const StLorentzVector<double>&    v, const StMatrix<float>& m1);
01514 template<> StLorentzVector<float>  operator*(const StMatrix<float>&  m1, const StLorentzVector<float>& v);
01515 template<> StLorentzVector<float>  operator*(const StLorentzVector<float>&     v, const StMatrix<float>& m1);
01516 
01517 template<> StMatrix<double> operator+(const StMatrix<double>& m1,const StMatrix<double>& m2);
01518 template<> StMatrix<double> operator+(const StMatrix<float>&  m1,const StMatrix<double>& m2);
01519 template<> StMatrix<double> operator+(const StMatrix<double>& m1,const StMatrix<float>& m2);
01520 template<> StMatrix<float>  operator+(const StMatrix<float>&  m1,const StMatrix<float>& m2);
01521 template<> StMatrix<double> operator-(const StMatrix<double>& m1,const StMatrix<double>& m2);
01522 template<> StMatrix<double> operator-(const StMatrix<float>&  m1,const StMatrix<double>& m2);
01523 template<> StMatrix<double> operator-(const StMatrix<double>& m1,const StMatrix<float>& m2);
01524 template<> StMatrix<float>  operator-(const StMatrix<float>&  m1,const StMatrix<float>& m2);
01525 
01526 template<> ostream& operator<<(ostream& s, const StMatrix<double>& q);
01527 template<> ostream& operator<<(ostream& s, const StMatrix<float>& q);
01528 
01529 template<> double norm_infinity(const StMatrix<double>& m1);
01530 template<> double normInfinity(const StMatrix<double>& m1);
01531 template<> double norm1(const StMatrix<double>& m1);
01532 template<> float norm_infinity(const StMatrix<float>& m1);
01533 template<> float normInfinity(const StMatrix<float>& m1);
01534 template<> float norm1(const StMatrix<float>& m1);
01535 #else /* ! __CINT__ */
01536 // Non-Member
01537 template<class DataType, class X>
01538 StMatrix<DataType> operator*(const StMatrix<DataType>& m1,const StMatrix<X>& m2)
01539 {
01540     return StMatrix<DataType>(m1).dot(m2);
01541 }
01542 
01543 // STL vector
01544 template <class DataType, class X>
01545 #ifdef ST_NO_TEMPLATE_DEF_ARGS
01546 vector<DataType, allocator<DataType> >
01547 operator*(const StMatrix<X>& m1, const vector<DataType, allocator<DataType> >& v)
01548 #else
01549 vector<DataType> operator*(const StMatrix<X>& m1, const vector<DataType>& v)
01550 #endif
01551 {
01552     vector<DataType> mret(m1.numRow());
01553     if(m1.numCol() == v.size() ) {
01554         for (unsigned int i=0; i<v.size(); i++) {
01555             register DataType temp = 0;
01556             for(unsigned int j=0; j<m1.numCol(); j++)
01557                 temp += m1(i+1,j+1)*v[j];
01558             mret[i]=temp;
01559         }
01560     }
01561     else {
01562 #ifndef ST_NO_EXCEPTIONS
01563         throw out_of_range("operator*(): StMatrix * STL vector Sizes are incompatible");
01564 #else
01565         cerr << "operator*(): StMatrix * STL vector Sizes are incompatible" << endl;
01566 #endif
01567     }
01568     return mret;
01569 }
01570 
01571 template <class DataType, class X>
01572 #ifdef ST_NO_TEMPLATE_DEF_ARGS
01573 vector<DataType, allocator<DataType> >
01574 operator*(const vector<DataType, allocator<DataType> >& v, const StMatrix<X>& m1)
01575 #else
01576 vector<DataType> operator*(const vector<DataType>& v, const StMatrix<X>& m1)
01577 #endif
01578 {
01579 
01580     if(v.size() == m1.numRow()) {
01581         vector<DataType> mret(m1.numCol());
01582         for (unsigned int i=0; i<m1.numCol(); i++) {
01583             register DataType temp = 0;
01584             for(unsigned int j=0; j<m1.numRow(); j++)
01585                 temp+=v[j]*m1(j+1,i+1);
01586             mret[i]=temp;
01587         }
01588         return mret;
01589     }
01590     else {
01591 #ifndef ST_NO_EXCEPTIONS
01592         throw out_of_range("operator*(): STL vector * StMatrix : Matrix Sizes are incompatible");
01593 #else
01594         cerr << "operator*(): STL vector * StMatrix : Matrix Sizes are incompatible" << endl;
01595 #endif
01596         return vector<DataType>();
01597     }
01598 }
01599 
01600 template <class DataType, class X>
01601 StThreeVector<DataType> operator*(const StMatrix<X>& m1, const StThreeVector<DataType>& v3)
01602 {
01603     if(m1.numRow() == 3 && m1.numCol() == 3) {
01604         return StThreeVector<DataType>(m1[0][0]*v3.x()+m1[0][1]*v3.y()+m1[0][2]*v3.z(),
01605                                        m1[1][0]*v3.x()+m1[1][1]*v3.y()+m1[1][2]*v3.z(),
01606                                        m1[2][0]*v3.x()+m1[2][1]*v3.y()+m1[2][2]*v3.z());
01607     }
01608     else {
01609 #ifndef ST_NO_EXCEPTIONS
01610         throw out_of_range("operator*(): StMatrix<> * StThreeVector<> : Matrix Must be 3x3.");
01611         
01612 #else
01613         cerr << "StMatrix<> * StThreeVector<>: Matrix Must be 3x3" << endl;
01614 #endif
01615         return StThreeVector<DataType>();
01616     }
01617 }
01618 
01619 
01620 template <class DataType, class X>
01621 StThreeVector<DataType> operator*(const StThreeVector<DataType>& v3, const StMatrix<X>& m1)
01622 {
01623     if(m1.numRow() == 3 && m1.numCol() == 3) {
01624         return StThreeVector<DataType>(m1[0][0]*v3.x()+m1[1][0]*v3.y()+m1[2][0]*v3.z(),
01625                                        m1[0][1]*v3.x()+m1[1][1]*v3.y()+m1[2][1]*v3.z(),
01626                                        m1[0][2]*v3.x()+m1[1][2]*v3.y()+m1[2][2]*v3.z());
01627     }
01628     else {
01629 #ifndef ST_NO_EXCEPTIONS
01630         throw out_of_range("operator*(): StThreeVector<> * StMatrix<>: Matrix Must be 3x3.");
01631         
01632 #else
01633         cerr << "operator*(): StThreeVector<> * StMatrix<>: Matrix Must be 3x3" << endl;
01634 #endif
01635         return StThreeVector<DataType>();
01636     }
01637 }
01638 
01639 template <class DataType, class X>
01640 StLorentzVector<DataType> operator*(const StMatrix<X>& m1, const StLorentzVector<DataType>& v4)
01641 {
01642     if (m1.numRow() == 4 && m1.numCol() == 4) {
01643         return StLorentzVector<DataType>(m1[0][0]*v4.x()+m1[0][1]*v4.y()+m1[0][2]*v4.z()+m1[0][3]*v4.t(),
01644                                          m1[1][0]*v4.x()+m1[1][1]*v4.y()+m1[1][2]*v4.z()+m1[1][3]*v4.t(),
01645                                          m1[2][0]*v4.x()+m1[2][1]*v4.y()+m1[2][2]*v4.z()+m1[2][3]*v4.t(),
01646                                          m1[3][0]*v4.x()+m1[3][1]*v4.y()+m1[3][2]*v4.z()+m1[3][3]*v4.t());
01647     }
01648     else {
01649 #ifndef ST_NO_EXCEPTIONS
01650         throw out_of_range("operator*(): StMatrix<> * StLorentzVector<> : Matrix Must be 4x4.");
01651         
01652 #else
01653         cerr << "operator*(): StMatrix<> * StLorentzVector<>: Matrix Must be 4x4" << endl;
01654 #endif
01655         return StLorentzVector<DataType>();
01656     }   
01657 }
01658 
01659 template <class DataType, class X>
01660 StLorentzVector<DataType> operator*(const StLorentzVector<DataType>& v4, const StMatrix<X>& m1)
01661 {
01662     if (m1.numRow() == 4 && m1.numCol() == 4) {
01663         return StLorentzVector<DataType>(m1[0][0]*v4.x()+m1[1][0]*v4.y()+m1[2][0]*v4.z()+m1[0][3]*v4.t(),
01664                                          m1[0][1]*v4.x()+m1[1][1]*v4.y()+m1[2][1]*v4.z()+m1[1][3]*v4.t(),
01665                                          m1[0][2]*v4.x()+m1[1][2]*v4.y()+m1[2][2]*v4.z()+m1[2][3]*v4.t(),
01666                                          m1[0][3]*v4.x()+m1[1][3]*v4.y()+m1[2][3]*v4.z()+m1[2][3]*v4.t());
01667     }
01668     else {
01669 #ifndef ST_NO_EXCEPTIONS
01670         throw out_of_range("operator*(): StLorentzVector<> * StMatrix<>: Matrix Must be 3x3.");
01671         
01672 #else
01673         cerr << "StLorentzVector<> * StMatrix<>: Matrix Must be 3x3" << endl;
01674 #endif
01675         return StLorentzVector<DataType>();
01676     }
01677 }
01678 
01679 template<class DataType, class X>
01680 StMatrix<DataType> operator+(const StMatrix<DataType>& m1,const StMatrix<X>& m2)
01681 {
01682     if(m1.numRow() == m2.numRow() && m1.numCol() == m2.numCol()) {
01683         StMatrix<DataType> mret(m1);
01684         mret +=m2;
01685         return mret;
01686     }
01687     else {
01688 #ifndef ST_NO_EXCEPTIONS
01689         throw out_of_range("operator+(): Matrix Sizes must be the same.");
01690         
01691 #else
01692         cerr << "operator+(): Matrix Sizes must be the same." << endl;
01693 #endif
01694         return StMatrix<DataType>();
01695     }   
01696 }
01697 
01698 template<class DataType, class X>
01699 StMatrix<DataType> operator-(const StMatrix<DataType>& m1,const StMatrix<X>& m2)
01700 {
01701     if(m1.numRow() == m2.numRow() && m1.numCol() == m2.numCol()) {
01702         return StMatrix<DataType>(m1) -= m2;
01703     }
01704     else {
01705 #ifndef ST_NO_EXCEPTIONS
01706         throw out_of_range("operator-(): Matrix Sizes must be the same.");
01707         
01708 #else
01709         cerr << "operator-(): Matrix Sizes must be the same." << endl;
01710 #endif
01711         return StMatrix<DataType>();
01712     }
01713 }
01714 
01715 // Print the Matrix.
01716 template<class DataType>
01717 ostream& operator<<(ostream& s, const StMatrix<DataType>& q)
01718 {
01719     s << "\n";
01720     // Fixed format needs 3 extra characters for field
01721     // Scientific format needs 7
01722     unsigned int width;
01723     if(s.flags()&ios::fixed)
01724         width = s.precision()+3;
01725     else
01726         width = s.precision()+7;
01727     for(unsigned int irow = 1; irow<= q.numRow(); irow++)
01728         {
01729             for(unsigned int icol=1; icol<=q.numCol(); icol++)
01730                 {
01731                     s.width(width);
01732                     s << q(irow,icol) << " ";
01733                 }
01734             s<< endl;
01735         }
01736     return s;
01737 }
01738 
01739 template<class DataType>
01740 DataType norm_infinity(const StMatrix<DataType>& m1)
01741 {
01742     return normInfinity(m1);
01743 }
01744 
01745 template<class DataType>
01746 DataType normInfinity(const StMatrix<DataType>& m1)
01747 {
01748     DataType max=0,sum;
01749     for(unsigned int r=1; r<=m1.numRow(); r++) {
01750         sum=0;
01751         for(unsigned int c=1; c<=m1.numCol(); c++) {
01752             sum+=fabs(m1(r,c));
01753         }
01754         if(sum>max) max=sum;
01755     }
01756     return max;
01757 }
01758 
01759 template<class DataType>
01760 DataType norm1(const StMatrix<DataType>& m1)
01761 {
01762     DataType max=0,sum;
01763     for(unsigned int c=1; c<=m1.numCol(); c++) {
01764         sum=0;
01765         for(unsigned int r=1; r<=m1.num_row(); r++)
01766             sum+=fabs(m1(r,c));
01767         if(sum>max) max=sum;
01768     }
01769     return max;
01770 }
01771 #endif /* __CINT__ */
01772 #endif /* ST_MATRIX_HH */