| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

LHCbMath

MatrixUtils.h

Go to the documentation of this file.
00001 // $Id: MatrixUtils.h,v 1.8 2008/02/24 19:40:21 ibelyaev Exp $
00002 // ============================================================================
00003 #ifndef LHCBMATH_MATRIXUTILS_H
00004 #define LHCBMATH_MATRIXUTILS_H 1
00005 // ============================================================================
00006 // Include files
00007 // ============================================================================
00008 // STD & STL 
00009 // ============================================================================
00010 #include <algorithm>
00011 #include <functional>
00012 #include <utility>
00013 #include <cmath>
00014 // ============================================================================
00015 // ROOT
00016 // ============================================================================
00017 #include "Math/SMatrix.h"
00018 #include "Math/SVector.h"
00019 #include "Math/Point3D.h"
00020 #include "Math/Vector4D.h"
00021 #include "Math/Vector3D.h"
00022 // ============================================================================
00044 // ============================================================================
00045 namespace Gaudi
00046 {
00047   namespace Math 
00048   {
00049     // ========================================================================
00069     template <class T, unsigned int D>
00070     inline size_t 
00071     setToScalar
00072     ( ROOT::Math::SVector<T,D>& m , const T& value = T() ) 
00073     { 
00074       std::fill ( m.begin() , m.end() , value ) ; 
00075       return D ;
00076     } 
00077     // ========================================================================
00098     template <class T, unsigned int D1, unsigned int D2, class R>
00099     inline size_t 
00100     setToScalar
00101     ( ROOT::Math::SMatrix<T,D1,D2,R>& m , const T& value = T() ) 
00102     { 
00103       std::fill ( m.begin() , m.end() , value ) ; 
00104       return m.end() - m.begin() ;
00105     } 
00106     // ========================================================================
00125     template <class T, unsigned int D, class R>
00126     inline size_t 
00127     setToUnit
00128     ( ROOT::Math::SMatrix<T,D,D,R>& m , const T& value = T(1) ) 
00129     { 
00131       std::fill ( m.begin() , m.end() , T(0.0) ) ; 
00133       for ( unsigned int i = 0 ; i < D ; ++i ) { m(i,i) = value ; }
00134       return m.end() - m.begin() ;
00135     }
00136     // ========================================================================
00158     template <class T, unsigned int D1, unsigned int D2, class R>
00159     inline size_t 
00160     scale 
00161     ( ROOT::Math::SMatrix<T,D1,D2,R>& m , const T& value ) 
00162     { 
00163       typedef typename ROOT::Math::SMatrix<T,D1,D2,R>::iterator iterator ;
00164       for ( iterator it = m.begin() ;  m.end() != it ; ++it ) { (*it) *= value ; }
00165       return m.end() - m.begin() ;
00166     } 
00167     // ========================================================================
00185     template <class T, unsigned int D>
00186     inline size_t 
00187     scale 
00188     ( ROOT::Math::SVector<T,D>& m , const T& value ) 
00189     { 
00190       typedef typename ROOT::Math::SVector<T,D>::iterator iterator ;
00191       for ( iterator it = m.begin() ;  m.end() != it ; ++it ) { (*it) *= value ; }
00192       return D ;
00193     } 
00194     // ========================================================================
00200     template <class T>
00201     struct _AbsCompare : public std::binary_function<T,T,bool>
00202     {
00203       inline bool operator() ( const T v1 , const T v2 ) const 
00204       { return ::fabs( v1 ) < ::fabs( v2 ) ; }
00205     } ;
00206     // ========================================================================
00213     template <class T,unsigned int D1,unsigned int D2, class R>
00214     inline T 
00215     max_element 
00216     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m ) 
00217     { return *std::max_element ( m.begin() , m.end() ) ; }
00218     // ========================================================================
00225     template <class T,unsigned int D1,unsigned int D2, class R>
00226     inline T 
00227     min_element 
00228     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m ) 
00229     { return *std::min_element ( m.begin() , m.end() ) ; }
00230     // ========================================================================
00237     template <class T,unsigned int D>
00238     inline T 
00239     max_element 
00240     ( const ROOT::Math::SVector<T,D>& m ) 
00241     { return *std::max_element ( m.begin() , m.end() ) ; }
00242     // ========================================================================
00249     template <class T,unsigned int D>
00250     inline T 
00251     min_element 
00252     ( const ROOT::Math::SVector<T,D>& m ) 
00253     { return *std::min_element ( m.begin() , m.end() ) ; }
00254     // ========================================================================
00276     template <class T,unsigned int D1,unsigned int D2, class R>
00277     inline T 
00278     maxabs_element 
00279     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m ) 
00280     { return *std::max_element ( m.begin() , m.end()  , _AbsCompare<T>() ) ; }
00281     // ========================================================================
00288     template <class T,unsigned int D1,unsigned int D2, class R>
00289     inline T 
00290     minabs_element 
00291     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m ) 
00292     { return *std::min_element ( m.begin() , m.end()  , _AbsCompare<T>() ) ; }
00293     // ========================================================================
00301     template <class T,unsigned int D1,unsigned int D2, class R, class CMP>
00302     inline std::pair<unsigned int,unsigned int> 
00303     ind_max_element 
00304     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m , CMP cmp ) 
00305     { 
00306       std::pair<unsigned int,unsigned int> result (0,0) ;
00307       T tmp = m(0,0) ;
00308       for ( unsigned int i = 0 ; i < D1 ; ++i ) 
00309       {
00310         for ( unsigned int j = 0 ; j < D2 ; ++j ) 
00311         {
00312           const T val = m(i,j) ;
00313           if ( !cmp( tmp , val ) ) { continue ; }
00314           tmp = val ; 
00315           result.first  = i ; 
00316           result.second = j ;
00317         }
00318       }
00319       return result ;
00320     } 
00321     // ========================================================================
00329     template <class T,unsigned int D, class CMP>
00330     inline std::pair<unsigned int,unsigned int> 
00331     ind_max_element 
00332     ( const ROOT::Math::SMatrix<T,D,D,ROOT::Math::MatRepSym<T,D> >& m , CMP cmp ) 
00333     { 
00334       std::pair<unsigned int,unsigned int> result (0,0) ;
00335       T tmp = m(0,0) ;
00336       for ( unsigned int i = 0 ; i < D ; ++i ) 
00337       {
00338         for ( unsigned int j = i ; j < D ; ++j ) 
00339         {
00340           const T val = m(i,j) ;
00341           if ( !cmp ( tmp , val ) ) { continue ; }
00342           tmp = val ; 
00343           result.first  = i ; 
00344           result.second = j ;
00345         }
00346       }
00347       return result ;
00348     }     
00349     // ========================================================================
00356     template <class T,unsigned int D1,unsigned int D2, class R>
00357     inline std::pair<unsigned int,unsigned int> 
00358     ind_max_element 
00359     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m ) 
00360     { return ind_max_element ( m , std::less<T>() ) ; } ;
00361     // ========================================================================
00369     template <class T,unsigned int D1,unsigned int D2, class R, class CMP>
00370     inline std::pair<unsigned int,unsigned int> 
00371     ind_min_element 
00372     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m , CMP cmp ) 
00373     { 
00374       std::pair<unsigned int,unsigned int> result (0,0) ;
00375       T tmp = m(0,0) ;
00376       for ( unsigned int i = 0 ; i < D1 ; ++i ) 
00377       {
00378         for ( unsigned int j = 0 ; j < D2 ; ++j ) 
00379         {
00380           const T val = m(i,j) ;
00381           if ( !cmp( val , tmp ) ) { continue ; }
00382           tmp = val ; 
00383           result.first  = i ; 
00384           result.second = j ;
00385         }
00386       }
00387       return result ;
00388     } 
00389     // ========================================================================
00397     template <class T,unsigned int D, class CMP>
00398     inline std::pair<unsigned int,unsigned int> 
00399     ind_min_element 
00400     ( const ROOT::Math::SMatrix<T,D,D,ROOT::Math::MatRepSym<T,D> >& m , CMP cmp ) 
00401     { 
00402       std::pair<unsigned int,unsigned int> result (0,0) ;
00403       T tmp = m(0,0) ;
00404       for ( unsigned int i = 0 ; i < D ; ++i ) 
00405       {
00406         for ( unsigned int j = i ; j < D ; ++j ) 
00407         {
00408           const T val = m(i,j) ;
00409           if ( !cmp ( tmp , val ) ) { continue ; }
00410           tmp = val ; 
00411           result.first  = i ; 
00412           result.second = j ;
00413         }
00414       }
00415       return result ;
00416     } 
00417     // ========================================================================
00440     template <class T, unsigned int D1, unsigned int D2, class R>
00441     inline std::pair<unsigned int,unsigned int> 
00442     ind_min_element 
00443     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m ) 
00444     { return ind_min_element( m , std::less<T>() ) ; } 
00445     // ========================================================================
00453     template <class T,unsigned int D, class CMP>
00454     inline unsigned int 
00455     ind_max_element 
00456     ( const ROOT::Math::SVector<T,D>& m , CMP cmp )
00457     { return std::max_element( m.begin() , m.end() , cmp ) - m.begin() ; }    
00458     // ========================================================================
00466     template <class T,unsigned int D, class CMP>
00467     inline unsigned int 
00468     ind_min_element 
00469     ( const ROOT::Math::SVector<T,D>& m , CMP cmp )
00470     { return std::min_element( m.begin() , m.end() , cmp ) - m.begin() ; }
00471     // ========================================================================
00478     template <class T,unsigned int D>
00479     inline unsigned int 
00480     ind_max_element 
00481     ( const ROOT::Math::SVector<T,D>& m )
00482     { return std::max_element( m.begin() , m.end() ) - m.begin() ; }    
00483     // ========================================================================
00490     template <class T,unsigned int D>
00491     inline unsigned int 
00492     ind_min_element 
00493     ( const ROOT::Math::SVector<T,D>& m )
00494     { return std::min_element ( m.begin() , m.end() ) - m.begin() ; }    
00495     // ========================================================================
00519     template <class T, unsigned int D1, unsigned int D2, class R>
00520     inline std::pair<unsigned int,unsigned int> 
00521     ind_maxabs_element 
00522     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m ) 
00523     { return ind_max_element( m , _AbsCompare<T>() ) ; } 
00524     // ========================================================================
00548     template <class T, unsigned int D1, unsigned int D2, class R>
00549     inline std::pair<unsigned int,unsigned int> 
00550     ind_minabs_element 
00551     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m ) 
00552     { return ind_min_element( m , _AbsCompare<T>() ) ; } 
00553     // ========================================================================
00560     template <class T,unsigned int D>
00561     inline unsigned int 
00562     ind_maxabs_element 
00563     ( const ROOT::Math::SVector<T,D>& m )
00564     { return ind_max_element ( m , _AbsCompare<T>() ) ; }    
00565     // ========================================================================
00572     template <class T,unsigned int D>
00573     inline unsigned int 
00574     ind_minabs_element 
00575     ( const ROOT::Math::SVector<T,D>& m )
00576     { return ind_min_element ( m , _AbsCompare<T>() ) ; }
00577     // ========================================================================
00597     template <class T, unsigned int D, class R>
00598     inline T 
00599     trace 
00600     ( const ROOT::Math::SMatrix<T,D,D,R>& m ) 
00601     {
00602       T result = m(0,0) ;
00603       for ( unsigned int i = 1 ; i < D ; ++i ) { result += m(i,i) ; }
00604       return result ;
00605     }
00625     template <class B, class T, unsigned int D, class R>
00626     inline T 
00627     trace 
00628     ( ROOT::Math::Expr<B,T,D,D,R>& m ) 
00629     {
00630       T result = m(0,0) ;
00631       for ( unsigned int i = 1 ; i < D ; ++i ) { result += m(i,i) ; }
00632       return result ;
00633     }
00634     // ========================================================================
00640     template <class T, unsigned int D, class R, class CMP>
00641     inline T 
00642     min_diagonal 
00643     ( const ROOT::Math::SMatrix<T,D,D,R>& m , CMP cmp ) 
00644     {
00645       T result = m(0,0);
00646       for ( unsigned int i = 1 ; i < D ; ++i ) 
00647       { 
00648         const T value = m(i,i) ;  
00649         if ( cmp ( value , result ) ) { result = value ; }
00650       }
00651       return result ;
00652     }
00653     // ========================================================================
00659     template <class T, unsigned int D, class R, class CMP>
00660     inline T 
00661     max_diagonal 
00662     ( const ROOT::Math::SMatrix<T,D,D,R>& m , CMP cmp ) 
00663     {
00664       T result = m(0,0);
00665       for ( unsigned int i = 1 ; i < D ; ++i ) 
00666       { 
00667         const T value = m(i,i) ;  
00668         if ( cmp ( result , value ) ) { result = value ; }
00669       }
00670       return result ;
00671     } 
00672     // ========================================================================
00693     template <class T, unsigned int D, class R>
00694     inline T 
00695     max_diagonal 
00696     ( const ROOT::Math::SMatrix<T,D,D,R>& m ) 
00697     { return max_diagonal( m , std::less<T>() ) ; } ;    
00698     // ========================================================================
00723     template <class T, unsigned int D, class R>
00724     inline T 
00725     min_diagonal 
00726     ( const ROOT::Math::SMatrix<T,D,D,R>& m ) 
00727     { return min_diagonal( m , std::less<T>() ) ; } ;    
00728     // ========================================================================
00749     template <class T, unsigned int D, class R>
00750     inline T 
00751     maxabs_diagonal 
00752     ( const ROOT::Math::SMatrix<T,D,D,R>& m ) 
00753     { return max_diagonal( m , _AbsCompare<T>() ) ; } ;
00754     // ========================================================================
00779     template <class T, unsigned int D, class R>
00780     inline T 
00781     minabs_diagonal 
00782     ( const ROOT::Math::SMatrix<T,D,D,R>& m ) 
00783     { return min_diagonal( m , _AbsCompare<T>() ) ; } ;
00784     // ========================================================================
00815     template <class T, unsigned int D1, unsigned int D2, class R, class P>
00816     inline size_t 
00817     count_if 
00818     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m , P pred )
00819     { return std::count_if ( m.begin() , m.end() , pred ) ; }
00820     // ========================================================================
00855     template <class T, unsigned int D, class P>
00856     inline size_t 
00857     count_if 
00858     ( const ROOT::Math::SMatrix<T,D,D,ROOT::Math::MatRepSym<T,D> >& m , P pred )
00859     { 
00860       size_t result = 0 ;
00861       for ( unsigned int i = 0 ; i < D ; ++i ) 
00862       {
00863         if ( pred ( m ( i , i ) ) ) { result += 1 ; }
00864         for ( unsigned int j = i + 1 ; j < D ; ++j ) 
00865         {
00866           if ( pred (  m ( i , j ) ) ) { result +=2 ; }  // ATTENTION! 
00867         }
00868       }
00869       return result ;
00870     } 
00871     // ========================================================================
00897     template <class T, unsigned int D, class R, class P>
00898     inline size_t 
00899     count_diagonal
00900     ( const ROOT::Math::SMatrix<T,D,D,R>& m , P pred )
00901     {
00902       size_t result = 0 ;
00903       for ( unsigned int i = 0 ; i < D ; ++i ) 
00904       { if ( pred ( m ( i , i ) ) ) { result += 1 ; } }
00905       return result ;
00906     } 
00907     // ========================================================================
00937     template <class T, unsigned int D1, unsigned int D2, class R, class P>
00938     inline bool
00939     check_if 
00940     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m , P pred )
00941     { return m.end() != std::find_if ( m.begin() , m.end() , pred ) ; }
00942     // ========================================================================
00968     template <class T, unsigned int D, class R, class P>
00969     inline bool
00970     check_diagonal 
00971     ( const ROOT::Math::SMatrix<T,D,D,R>& m , P pred )
00972     { 
00973       for ( unsigned int i = 0 ; i < D ; ++i ) 
00974       { if ( pred ( m ( i , i ) ) ) { return true ; } }
00975       return false ;
00976     } 
00977     // ========================================================================
01015     template <class T1 , class T2 ,
01016               unsigned int D1 , unsigned int D2 ,
01017               class R1 , class R2 , class P>
01018     inline bool 
01019     equal_if  
01020     ( const ROOT::Math::SMatrix<T1,D1,D2,R1>& m1 , 
01021       const ROOT::Math::SMatrix<T2,D1,D2,R2>& m2 , P pred )
01022     { 
01023       for ( unsigned int i = 0 ; i < D1 ; ++i ) 
01024       {
01025         for ( unsigned int j = 0 ; j < D2 ; ++j ) 
01026         {
01027           if ( !pred( m1(i,j) , m2(i,j) ) ) { return false ; } // RETURN 
01028         }
01029       }
01030       return true ;
01031     } 
01032     // ========================================================================
01068     template <class T,unsigned int D1, unsigned int D2,class R,class P>
01069     inline bool 
01070     equal_if  
01071     ( const ROOT::Math::SMatrix<T,D1,D2,R>& m1 , 
01072       const ROOT::Math::SMatrix<T,D1,D2,R>& m2 , P pred )
01073     { return std::equal ( m1.begin() , m1.end() , m2.begin() , pred ) ; } 
01074     // ========================================================================
01075 
01076 
01077 
01078     // =========================================================================
01087     template <class T, class T2, unsigned int D>
01088     void update 
01089     ( ROOT::Math::SMatrix<T,D,D,ROOT::Math::MatRepSym<T,D> > & left  , 
01090       const ROOT::Math::SVector<T2,D>&                         vect  , 
01091       const double                                             scale ) 
01092     {
01093       for ( unsigned int i = 0 ; i < D ; ++i ) 
01094       {
01095         for ( unsigned int j = i ; j < D ; ++j ) 
01096         { left ( i , j ) += scale * vect(i) * vect(j) ; }
01097       }
01098     }
01099     // =========================================================================
01108     template <class T, class B, class T2, unsigned int D>
01109     void update 
01110     ( ROOT::Math::SMatrix<T,D,D,ROOT::Math::MatRepSym<T,D> > & left  , 
01111       const ROOT::Math::VecExpr<B,T2,D>&                       vect  , 
01112       const double                                             scale ) 
01113     {
01114       for ( unsigned int i = 0 ; i < D ; ++i ) 
01115       {
01116         for ( unsigned int j = i ; j < D ; ++j ) 
01117         { left ( i , j ) += scale * vect(i) * vect(j) ; }
01118       }
01119     }
01120     // =========================================================================
01130     template <class T, class R, class T2, class T3, unsigned int D1,unsigned int D2>
01131     void update 
01132     ( ROOT::Math::SMatrix<T,D1,D2,R>  & left        , 
01133       const ROOT::Math::SVector<T2,D1>& vct1        , 
01134       const ROOT::Math::SVector<T3,D2>& vct2        , 
01135       const double                      scale = 1.0 ) 
01136     {
01137       for ( unsigned int i = 0 ; i < D1 ; ++i ) 
01138       {
01139         for ( unsigned int j = 0 ; j < D2 ; ++j ) 
01140         { left ( i , j ) += scale * vct1(i) * vct2(j) ; }
01141       }
01142     }
01143     // =========================================================================
01152     template<class T, class T1, class T2, class R, unsigned int D1, unsigned int D2>
01153     T mult 
01154     ( const ROOT::Math::SVector<T1,D1>&     vct1  ,
01155       const ROOT::Math::SMatrix<T,D1,D2,R>& mtrx  , 
01156       const ROOT::Math::SVector<T2,D2>&     vct2  ) 
01157     { 
01158       return ROOT::Math::Dot ( vct1 , mtrx * vct2 ) ; 
01159     }
01160     // =========================================================================
01169     template <class T, class T2, class R, unsigned int D>
01170     void update 
01171     ( ROOT::Math::SMatrix<T,D,D,ROOT::Math::MatRepSym<T,D> >&   left        , 
01172       const ROOT::Math::SMatrix<T2,D,D,R>&                      right       , 
01173       const double                                              scale = 1.0 ) 
01174     {
01175       for ( unsigned int i = 0 ; i < D ; ++i ) 
01176       {
01177         for ( unsigned int j = i ; j < D ; ++j ) 
01178         { left ( i , j ) += scale * ( right ( i , j ) + right ( j , i ) ) ; }
01179       }
01180     }
01181     // =========================================================================
01190     template <class T, class T2, class B, class R, unsigned int D>
01191     void update 
01192     ( ROOT::Math::SMatrix<T,D,D,ROOT::Math::MatRepSym<T,D> >&   left        , 
01193       const ROOT::Math::Expr<B,T2,D,D,R>&                       right       , 
01194       const double                                              scale = 1.0 ) 
01195     {
01196       for ( unsigned int i = 0 ; i < D ; ++i ) 
01197       {
01198         for ( unsigned int j = i ; j < D ; ++j ) 
01199         { left ( i , j ) += scale * ( right ( i , j ) + right ( j , i ) ) ; }
01200       }
01201     }
01202     // ========================================================================
01203   } // end of namespace Math  
01204 } // end of namespace Gaudi
01205 
01206 // ============================================================================
01207 // The END 
01208 // ============================================================================
01209 #endif // LHCBMATH_MATRIXUTILS_H
01210 // ============================================================================
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |
Generated on Mon Apr 11 20:02:57 2011 for LHCbMath by doxygen 1.4.7