ROOT_528-00b_version: include/Math/MatrixFunctions.h Source File

00001 // @(#)root/smatrix:$Id: MatrixFunctions.h 30749 2009-10-15 16:33:04Z brun $
00002 // Authors: T. Glebe, L. Moneta    2005  
00003 
00004 #ifndef ROOT_Math_MatrixFunctions
00005 #define ROOT_Math_MatrixFunctions
00006 // ********************************************************************
00007 //
00008 // source:
00009 //
00010 // type:      source code
00011 //
00012 // created:   20. Mar 2001
00013 //
00014 // author:    Thorsten Glebe
00015 //            HERA-B Collaboration
00016 //            Max-Planck-Institut fuer Kernphysik
00017 //            Saupfercheckweg 1
00018 //            69117 Heidelberg
00019 //            Germany
00020 //            E-mail: T.Glebe@mpi-hd.mpg.de
00021 //
00022 // Description: Functions/Operators special to Matrix
00023 //
00024 // changes:
00025 // 20 Mar 2001 (TG) creation
00026 // 20 Mar 2001 (TG) Added Matrix * Vector multiplication
00027 // 21 Mar 2001 (TG) Added transpose, product
00028 // 11 Apr 2001 (TG) transpose() speed improvment by removing rows(), cols()
00029 //                  through static members of Matrix and Expr class
00030 //
00031 // ********************************************************************
00032 
00033 //doxygen tag
00034 /**
00035    @defgroup MatrixFunctions Matrix Template Functions 
00036    @ingroup SMatrixGroup
00037 
00038    These function apply to matrices (and also Matrix expression) and can return a 
00039    matrix expression of a particular defined type, like in the matrix multiplication or 
00040    a vector, like in the matrix-vector product or a scalar like in the Similarity vector-matrix product.
00041 */
00042 
00043 #ifndef ROOT_Math_BinaryOpPolicy
00044 #include "Math/BinaryOpPolicy.h"
00045 #endif
00046 
00047 namespace ROOT { 
00048 
00049   namespace Math { 
00050 
00051 
00052 
00053 #ifdef XXX
00054 //==============================================================================
00055 // SMatrix * SVector
00056 //==============================================================================
00057 template <class T, unsigned int D1, unsigned int D2, class R>
00058 SVector<T,D1> operator*(const SMatrix<T,D1,D2,R>& rhs, const SVector<T,D2>& lhs)
00059 {
00060   SVector<T,D1> tmp;
00061   for(unsigned int i=0; i<D1; ++i) {
00062     const unsigned int rpos = i*D2;
00063     for(unsigned int j=0; j<D2; ++j) {
00064       tmp[i] += rhs.apply(rpos+j) * lhs.apply(j);
00065     }
00066   }
00067   return tmp;
00068 }
00069 #endif
00070 
00071 
00072 // matrix-vector product: 
00073 // use apply(i) function for matrices. Tested  (11/05/06) with using (i,j) but 
00074 // performances are slightly worse (not clear  why)
00075 
00076 //==============================================================================
00077 // meta_row_dot
00078 //==============================================================================
00079 template <unsigned int I>
00080 struct meta_row_dot {
00081   template <class A, class B>
00082   static inline typename A::value_type f(const A& lhs, const B& rhs,
00083                                          const unsigned int offset) {
00084     return lhs.apply(offset+I) * rhs.apply(I) + meta_row_dot<I-1>::f(lhs,rhs,offset);
00085   }
00086 };
00087 
00088 
00089 //==============================================================================
00090 // meta_row_dot<0>
00091 //==============================================================================
00092 template <>
00093 struct meta_row_dot<0> {
00094   template <class A, class B>
00095   static inline typename A::value_type f(const A& lhs, const B& rhs,
00096                                          const unsigned int offset) {
00097     return lhs.apply(offset) * rhs.apply(0);
00098   }
00099 };
00100 
00101 //==============================================================================
00102 // VectorMatrixRowOp
00103 //==============================================================================
00104 template <class Matrix, class Vector, unsigned int D2>
00105 class VectorMatrixRowOp {
00106 public:
00107   ///
00108   VectorMatrixRowOp(const Matrix& lhs, const Vector& rhs) :
00109     lhs_(lhs), rhs_(rhs) {}
00110 
00111   ///
00112   ~VectorMatrixRowOp() {}
00113 
00114   /// calc \f$ \sum_{j} a_{ij} * v_j \f$
00115   inline typename Matrix::value_type apply(unsigned int i) const {
00116     return meta_row_dot<D2-1>::f(lhs_, rhs_, i*D2);
00117   }
00118 
00119 protected:
00120   const Matrix& lhs_;
00121   const Vector& rhs_;
00122 };
00123 
00124 
00125 //==============================================================================
00126 // meta_col_dot
00127 //==============================================================================
00128 template <unsigned int I>
00129 struct meta_col_dot {
00130   template <class Matrix, class Vector>
00131   static inline typename Matrix::value_type f(const Matrix& lhs, const Vector& rhs,
00132                                               const unsigned int offset) {
00133     return lhs.apply(Matrix::kCols*I+offset) * rhs.apply(I) + 
00134            meta_col_dot<I-1>::f(lhs,rhs,offset);
00135   }
00136 };
00137 
00138 
00139 //==============================================================================
00140 // meta_col_dot<0>
00141 //==============================================================================
00142 template <>
00143 struct meta_col_dot<0> {
00144   template <class Matrix, class Vector>
00145   static inline typename Matrix::value_type f(const Matrix& lhs, const Vector& rhs,
00146                                               const unsigned int offset) {
00147     return lhs.apply(offset) * rhs.apply(0);
00148   }
00149 };
00150 
00151 //==============================================================================
00152 // VectorMatrixColOp
00153 //==============================================================================
00154 /**
00155    Class for Vector-Matrix multiplication
00156 
00157    @ingroup Expression
00158  */
00159 template <class Vector, class Matrix, unsigned int D1>
00160 class VectorMatrixColOp {
00161 public:
00162   ///
00163   VectorMatrixColOp(const Vector& lhs, const Matrix& rhs) :
00164     lhs_(lhs), rhs_(rhs) {}
00165 
00166   ///
00167   ~VectorMatrixColOp() {}
00168 
00169   /// calc \f$ \sum_{j} a_{ij} * v_j \f$
00170   inline typename Matrix::value_type apply(unsigned int i) const {
00171     return meta_col_dot<D1-1>::f(rhs_, lhs_, i);
00172   }
00173 
00174 protected:
00175   const Vector&    lhs_;
00176   const Matrix&    rhs_;
00177 };
00178 
00179 /**
00180    Matrix *  Vector multiplication   \f$ a(i) = \sum_{j} M(i,j) * b(j) \f$
00181    returning a vector expression 
00182 
00183    @ingroup MatrixFunctions
00184  */
00185 //==============================================================================
00186 // operator*: SMatrix * SVector
00187 //==============================================================================
00188 template <class T, unsigned int D1, unsigned int D2, class R>
00189 inline VecExpr<VectorMatrixRowOp<SMatrix<T,D1,D2,R>,SVector<T,D2>, D2>, T, D1>
00190  operator*(const SMatrix<T,D1,D2,R>& lhs, const SVector<T,D2>& rhs) {
00191   typedef VectorMatrixRowOp<SMatrix<T,D1,D2,R>,SVector<T,D2>, D2> VMOp;
00192   return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
00193 }
00194 
00195 //==============================================================================
00196 // operator*: SMatrix * Expr<A,T,D2>
00197 //==============================================================================
00198 template <class A, class T, unsigned int D1, unsigned int D2, class R>
00199 inline VecExpr<VectorMatrixRowOp<SMatrix<T,D1,D2,R>, VecExpr<A,T,D2>, D2>, T, D1>
00200  operator*(const SMatrix<T,D1,D2,R>& lhs, const VecExpr<A,T,D2>& rhs) {
00201   typedef VectorMatrixRowOp<SMatrix<T,D1,D2,R>,VecExpr<A,T,D2>, D2> VMOp;
00202   return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
00203 }
00204 
00205 //==============================================================================
00206 // operator*: Expr<A,T,D1,D2> * SVector
00207 //==============================================================================
00208 template <class A, class T, unsigned int D1, unsigned int D2, class R>
00209 inline VecExpr<VectorMatrixRowOp<Expr<A,T,D1,D2,R>, SVector<T,D2>, D2>, T, D1>
00210  operator*(const Expr<A,T,D1,D2,R>& lhs, const SVector<T,D2>& rhs) {
00211   typedef VectorMatrixRowOp<Expr<A,T,D1,D2,R>,SVector<T,D2>, D2> VMOp;
00212   return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
00213 }
00214 
00215 //==============================================================================
00216 // operator*: Expr<A,T,D1,D2> * VecExpr<B,T,D2>
00217 //==============================================================================
00218 template <class A, class B, class T, unsigned int D1, unsigned int D2, class R>
00219 inline VecExpr<VectorMatrixRowOp<Expr<A,T,D1,D2,R>, VecExpr<B,T,D2>, D2>, T, D1>
00220  operator*(const Expr<A,T,D1,D2,R>& lhs, const VecExpr<B,T,D2>& rhs) {
00221   typedef VectorMatrixRowOp<Expr<A,T,D1,D2,R>,VecExpr<B,T,D2>, D2> VMOp;
00222   return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
00223 }
00224 
00225 //==============================================================================
00226 // operator*: SVector * SMatrix
00227 //==============================================================================
00228 template <class T, unsigned int D1, unsigned int D2, class R>
00229 inline VecExpr<VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2,R>, D1>, T, D2>
00230  operator*(const SVector<T,D1>& lhs, const SMatrix<T,D1,D2,R>& rhs) {
00231   typedef VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2,R>, D1> VMOp;
00232   return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
00233 }
00234 
00235 //==============================================================================
00236 // operator*: SVector * Expr<A,T,D1,D2>
00237 //==============================================================================
00238 template <class A, class T, unsigned int D1, unsigned int D2, class R>
00239 inline VecExpr<VectorMatrixColOp<SVector<T,D1>, Expr<A,T,D1,D2,R>, D1>, T, D2>
00240  operator*(const SVector<T,D1>& lhs, const Expr<A,T,D1,D2,R>& rhs) {
00241   typedef VectorMatrixColOp<SVector<T,D1>, Expr<A,T,D1,D2,R>, D1> VMOp;
00242   return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
00243 }
00244 
00245 //==============================================================================
00246 // operator*: VecExpr<A,T,D1> * SMatrix
00247 //==============================================================================
00248 template <class A, class T, unsigned int D1, unsigned int D2, class R>
00249 inline VecExpr<VectorMatrixColOp<VecExpr<A,T,D1>, SMatrix<T,D1,D2,R>, D1>, T, D2>
00250  operator*(const VecExpr<A,T,D1>& lhs, const SMatrix<T,D1,D2,R>& rhs) {
00251   typedef VectorMatrixColOp<VecExpr<A,T,D1>, SMatrix<T,D1,D2,R>, D1> VMOp;
00252   return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
00253 }
00254 
00255 //==============================================================================
00256 // operator*: VecExpr<A,T,D1> * Expr<B,T,D1,D2>
00257 //==============================================================================
00258 template <class A, class B, class T, unsigned int D1, unsigned int D2, class R>
00259 inline VecExpr<VectorMatrixColOp<VecExpr<A,T,D1>, Expr<B,T,D1,D2,R>, D1>, T, D2>
00260  operator*(const VecExpr<A,T,D1>& lhs, const Expr<B,T,D1,D2,R>& rhs) {
00261   typedef VectorMatrixColOp<VecExpr<A,T,D1>, Expr<B,T,D1,D2,R>, D1> VMOp;
00262   return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
00263 }
00264 
00265 //==============================================================================
00266 // meta_matrix_dot
00267 //==============================================================================
00268 template <unsigned int I>
00269 struct meta_matrix_dot {
00270 
00271   template <class MatrixA, class MatrixB>
00272   static inline typename MatrixA::value_type f(const MatrixA& lhs, 
00273                                                const MatrixB& rhs,
00274                                                const unsigned int offset) {
00275     return lhs.apply(offset/MatrixB::kCols*MatrixA::kCols + I) *
00276            rhs.apply(MatrixB::kCols*I + offset%MatrixB::kCols) + 
00277            meta_matrix_dot<I-1>::f(lhs,rhs,offset);
00278   }
00279 
00280   // multiplication using i and j indeces
00281   template <class MatrixA, class MatrixB>
00282   static inline typename MatrixA::value_type g(const MatrixA& lhs, 
00283                                                const MatrixB& rhs,
00284                                                unsigned int i, 
00285                                                unsigned int j) {
00286     return lhs(i, I) * rhs(I , j) + 
00287            meta_matrix_dot<I-1>::g(lhs,rhs,i,j);
00288   }
00289 };
00290 
00291 
00292 //==============================================================================
00293 // meta_matrix_dot<0>
00294 //==============================================================================
00295 template <>
00296 struct meta_matrix_dot<0> {
00297 
00298   template <class MatrixA, class MatrixB>
00299   static inline typename MatrixA::value_type f(const MatrixA& lhs, 
00300                                                const MatrixB& rhs,
00301                                                const unsigned int offset) {
00302     return lhs.apply(offset/MatrixB::kCols*MatrixA::kCols) *
00303            rhs.apply(offset%MatrixB::kCols);
00304   }
00305 
00306   // multiplication using i and j 
00307   template <class MatrixA, class MatrixB>
00308   static inline typename MatrixA::value_type g(const MatrixA& lhs, 
00309                                                const MatrixB& rhs,
00310                                                unsigned int i, unsigned int j) {
00311     return lhs(i,0) * rhs(0,j);
00312   }
00313 
00314 };
00315 
00316 //==============================================================================
00317 // MatrixMulOp
00318 //==============================================================================
00319 /**
00320    Class for Matrix-Matrix multiplication 
00321 
00322    @ingroup Expression
00323  */
00324 template <class MatrixA, class MatrixB, class T, unsigned int D>
00325 class MatrixMulOp {
00326 public:
00327   ///
00328   MatrixMulOp(const MatrixA& lhs, const MatrixB& rhs) :
00329     lhs_(lhs), rhs_(rhs) {}
00330 
00331   ///
00332   ~MatrixMulOp() {}
00333 
00334   /// calc \f$\sum_{j} a_{ik} * b_{kj}\f$
00335   inline T apply(unsigned int i) const {
00336     return meta_matrix_dot<D-1>::f(lhs_, rhs_, i);
00337   }
00338 
00339   inline T operator() (unsigned int i, unsigned j) const {
00340     return meta_matrix_dot<D-1>::g(lhs_, rhs_, i, j);
00341   }
00342 
00343   inline bool IsInUse (const T * p) const { 
00344     return lhs_.IsInUse(p) || rhs_.IsInUse(p);
00345   }
00346 
00347 
00348 protected:
00349   const MatrixA&    lhs_;
00350   const MatrixB&    rhs_;
00351 };
00352 
00353 
00354 /**
00355    Matrix *  Matrix multiplication , \f$ C(i,j) = \sum_{k} A(i,k) * B(k,j)\f$ 
00356    returning a matrix expression 
00357 
00358    @ingroup MatrixFunctions
00359  */
00360 //==============================================================================
00361 // operator* (SMatrix * SMatrix, binary)
00362 //==============================================================================
00363 template <  class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
00364 inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
00365  operator*(const SMatrix<T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
00366   typedef MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, T,D> MatMulOp;
00367   return Expr<MatMulOp,T,D1,D2,
00368     typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
00369 }
00370 
00371 //==============================================================================
00372 // operator* (SMatrix * Expr, binary)
00373 //==============================================================================
00374 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
00375 inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
00376  operator*(const SMatrix<T,D1,D,R1>& lhs, const Expr<A,T,D,D2,R2>& rhs) {
00377   typedef MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>,T,D> MatMulOp;
00378   return Expr<MatMulOp,T,D1,D2,
00379     typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
00380 }
00381 
00382 //==============================================================================
00383 // operator* (Expr * SMatrix, binary)
00384 //==============================================================================
00385 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
00386 inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
00387  operator*(const Expr<A,T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
00388   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>,T,D> MatMulOp;
00389   return Expr<MatMulOp,T,D1,D2,
00390     typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
00391 }
00392 
00393 //==============================================================================
00394 // operator* (Expr * Expr, binary)
00395 //==============================================================================
00396 template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
00397 inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
00398  operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
00399   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, T,D> MatMulOp;
00400   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
00401 }
00402 
00403 
00404 
00405 #ifdef XXX
00406 //==============================================================================
00407 // MatrixMulOp
00408 //==============================================================================
00409 template <class MatrixA, class MatrixB, unsigned int D>
00410 class MatrixMulOp {
00411 public:
00412   ///
00413   MatrixMulOp(const MatrixA& lhs, const MatrixB& rhs) :
00414     lhs_(lhs), rhs_(rhs) {}
00415 
00416   ///
00417   ~MatrixMulOp() {}
00418 
00419   /// calc $\sum_{j} a_{ik} * b_{kj}$
00420   inline typename MatrixA::value_type apply(unsigned int i) const {
00421     return meta_matrix_dot<D-1>::f(lhs_, rhs_, i);
00422   }
00423 
00424 protected:
00425   const MatrixA&    lhs_;
00426   const MatrixB&    rhs_;
00427 };
00428 
00429 
00430 //==============================================================================
00431 // operator* (SMatrix * SMatrix, binary)
00432 //==============================================================================
00433 template <  class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
00434 inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
00435  operator*(const SMatrix<T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
00436   typedef MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, D> MatMulOp;
00437   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
00438 }
00439 
00440 //==============================================================================
00441 // operator* (SMatrix * Expr, binary)
00442 //==============================================================================
00443 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
00444 inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
00445  operator*(const SMatrix<T,D1,D,R1>& lhs, const Expr<A,T,D,D2,R2>& rhs) {
00446   typedef MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>, D> MatMulOp;
00447   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
00448 }
00449 
00450 //==============================================================================
00451 // operator* (Expr * SMatrix, binary)
00452 //==============================================================================
00453 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
00454 inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
00455  operator*(const Expr<A,T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
00456   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>, D> MatMulOp;
00457   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
00458 }
00459 
00460 //=============================================================================
00461 // operator* (Expr * Expr, binary)
00462 //=============================================================================
00463 template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
00464 inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
00465  operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
00466   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, D> MatMulOp;
00467   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
00468 }
00469 #endif
00470 
00471 //==============================================================================
00472 // TransposeOp
00473 //==============================================================================
00474 /**
00475    Class for Transpose Operations
00476 
00477    @ingroup Expression
00478  */
00479 template <class Matrix, class T, unsigned int D1, unsigned int D2=D1>
00480 class TransposeOp {
00481 public:
00482   ///
00483   TransposeOp( const Matrix& rhs) :
00484     rhs_(rhs) {}
00485 
00486   ///
00487   ~TransposeOp() {}
00488 
00489   ///
00490   inline T apply(unsigned int i) const {
00491     return rhs_.apply( (i%D1)*D2 + i/D1);
00492   }
00493   inline T operator() (unsigned int i, unsigned j) const {
00494     return rhs_( j, i);
00495   }
00496 
00497   inline bool IsInUse (const T * p) const { 
00498     return rhs_.IsInUse(p); 
00499   }
00500 
00501 protected:
00502   const Matrix& rhs_;
00503 };
00504 
00505 
00506 /**
00507    Matrix Transpose   B(i,j) = A(j,i) 
00508    returning a matrix expression
00509 
00510    @ingroup MatrixFunctions
00511  */
00512 //==============================================================================
00513 // transpose
00514 //==============================================================================
00515 template <class T, unsigned int D1, unsigned int D2, class R>
00516 inline Expr<TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2>, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
00517  Transpose(const SMatrix<T,D1,D2, R>& rhs) {
00518   typedef TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2> MatTrOp;
00519 
00520   return Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs));
00521 }
00522 
00523 //==============================================================================
00524 // transpose
00525 //==============================================================================
00526 template <class A, class T, unsigned int D1, unsigned int D2, class R>
00527 inline Expr<TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2>, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
00528  Transpose(const Expr<A,T,D1,D2,R>& rhs) {
00529   typedef TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2> MatTrOp;
00530 
00531   return Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs));
00532 }
00533 
00534 
00535 #ifdef ENABLE_TEMPORARIES_TRANSPOSE 
00536 // sometimes is faster to create a temp, not clear why
00537 
00538 //==============================================================================
00539 // transpose
00540 //==============================================================================
00541 template <class T, unsigned int D1, unsigned int D2, class R>
00542 inline SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
00543  Transpose(const SMatrix<T,D1,D2, R>& rhs) {
00544   typedef TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2> MatTrOp;
00545 
00546   return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType> 
00547     ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
00548 }
00549 
00550 //==============================================================================
00551 // transpose
00552 //==============================================================================
00553 template <class A, class T, unsigned int D1, unsigned int D2, class R>
00554 inline  SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
00555  Transpose(const Expr<A,T,D1,D2,R>& rhs) {
00556   typedef TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2> MatTrOp;
00557 
00558   return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType> 
00559     ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
00560 }
00561 
00562 #endif 
00563 
00564 
00565 #ifdef OLD
00566 //==============================================================================
00567 // product: SMatrix/SVector calculate v^T * A * v
00568 //==============================================================================
00569 template <class T, unsigned int D, class R>
00570 inline T Product(const SMatrix<T,D,D,R>& lhs, const SVector<T,D>& rhs) {
00571   return Dot(rhs, lhs * rhs);
00572 }
00573 
00574 //==============================================================================
00575 // product: SVector/SMatrix calculate v^T * A * v
00576 //==============================================================================
00577 template <class T, unsigned int D, class R>
00578 inline T Product(const SVector<T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
00579   return Dot(lhs, rhs * lhs);
00580 }
00581 
00582 //==============================================================================
00583 // product: SMatrix/Expr calculate v^T * A * v
00584 //==============================================================================
00585 template <class A, class T, unsigned int D, class R>
00586 inline T Product(const SMatrix<T,D,D,R>& lhs, const VecExpr<A,T,D>& rhs) {
00587   return Dot(rhs, lhs * rhs);
00588 }
00589 
00590 //==============================================================================
00591 // product: Expr/SMatrix calculate v^T * A * v
00592 //==============================================================================
00593 template <class A, class T, unsigned int D, class R>
00594 inline T Product(const VecExpr<A,T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
00595   return Dot(lhs, rhs * lhs);
00596 }
00597 
00598 //==============================================================================
00599 // product: SVector/Expr calculate v^T * A * v
00600 //==============================================================================
00601 template <class A, class T, unsigned int D, class R>
00602 inline T Product(const SVector<T,D>& lhs, const Expr<A,T,D,D,R>& rhs) {
00603   return Dot(lhs, rhs * lhs);
00604 }
00605 
00606 //==============================================================================
00607 // product: Expr/SVector calculate v^T * A * v
00608 //==============================================================================
00609 template <class A, class T, unsigned int D, class R>
00610 inline T Product(const Expr<A,T,D,D,R>& lhs, const SVector<T,D>& rhs) {
00611   return Dot(rhs, lhs * rhs);
00612 }
00613 
00614 //==============================================================================
00615 // product: Expr/Expr calculate v^T * A * v
00616 //==============================================================================
00617 template <class A, class B, class T, unsigned int D, class R>
00618 inline T Product(const Expr<A,T,D,D,R>& lhs, const VecExpr<B,T,D>& rhs) {
00619   return Dot(rhs, lhs * rhs);
00620 }
00621 
00622 //==============================================================================
00623 // product: Expr/Expr calculate v^T * A * v
00624 //==============================================================================
00625 template <class A, class B, class T, unsigned int D, class R>
00626 inline T Product(const VecExpr<A,T,D>& lhs, const Expr<B,T,D,D,R>& rhs) {
00627   return Dot(lhs, rhs * lhs);
00628 }
00629 #endif
00630 
00631 /**
00632    Similarity Vector - Matrix Product:  v^T * A * v
00633    returning a scalar value of type T   \f$ s = \sum_{i,j} v(i) * A(i,j) * v(j)\f$ 
00634 
00635    @ingroup MatrixFunctions
00636  */
00637 
00638 //==============================================================================
00639 // product: SMatrix/SVector calculate v^T * A * v
00640 //==============================================================================
00641 template <class T, unsigned int D, class R>
00642 inline T Similarity(const SMatrix<T,D,D,R>& lhs, const SVector<T,D>& rhs) {
00643   return Dot(rhs, lhs * rhs);
00644 }
00645 
00646 //==============================================================================
00647 // product: SVector/SMatrix calculate v^T * A * v
00648 //==============================================================================
00649 template <class T, unsigned int D, class R>
00650 inline T Similarity(const SVector<T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
00651   return Dot(lhs, rhs * lhs);
00652 }
00653 
00654 //==============================================================================
00655 // product: SMatrix/Expr calculate v^T * A * v
00656 //==============================================================================
00657 template <class A, class T, unsigned int D, class R>
00658 inline T Similarity(const SMatrix<T,D,D,R>& lhs, const VecExpr<A,T,D>& rhs) {
00659   return Dot(rhs, lhs * rhs);
00660 }
00661 
00662 //==============================================================================
00663 // product: Expr/SMatrix calculate v^T * A * v
00664 //==============================================================================
00665 template <class A, class T, unsigned int D, class R>
00666 inline T Similarity(const VecExpr<A,T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
00667   return Dot(lhs, rhs * lhs);
00668 }
00669 
00670 //==============================================================================
00671 // product: SVector/Expr calculate v^T * A * v
00672 //==============================================================================
00673 template <class A, class T, unsigned int D, class R>
00674 inline T Similarity(const SVector<T,D>& lhs, const Expr<A,T,D,D,R>& rhs) {
00675   return Dot(lhs, rhs * lhs);
00676 }
00677 
00678 //==============================================================================
00679 // product: Expr/SVector calculate v^T * A * v
00680 //==============================================================================
00681 template <class A, class T, unsigned int D, class R>
00682 inline T Similarity(const Expr<A,T,D,D,R>& lhs, const SVector<T,D>& rhs) {
00683   return Dot(rhs, lhs * rhs);
00684 }
00685 
00686 //==============================================================================
00687 // product: Expr/Expr calculate v^T * A * v
00688 //==============================================================================
00689 template <class A, class B, class T, unsigned int D, class R>
00690 inline T Similarity(const Expr<A,T,D,D,R>& lhs, const VecExpr<B,T,D>& rhs) {
00691   return Dot(rhs, lhs * rhs);
00692 }
00693 
00694 //==============================================================================
00695 // product: Expr/Expr calculate v^T * A * v
00696 //==============================================================================
00697 template <class A, class B, class T, unsigned int D, class R>
00698 inline T Similarity(const VecExpr<A,T,D>& lhs, const Expr<B,T,D,D,R>& rhs) {
00699   return Dot(lhs, rhs * lhs);
00700 }
00701 
00702 
00703 /**
00704    Similarity Matrix Product :  B = U * A * U^T for A symmetric 
00705    returning a symmetric matrix expression: 
00706    \f$ B(i,j) = \sum_{k,l} U(i,k) * A(k,l) * U(j,l) \f$ 
00707 
00708    @ingroup MatrixFunctions
00709  */
00710 //==============================================================================
00711 // product: SMatrix/SMatrix calculate M * A * M^T where A is a symmetric matrix
00712 // return matrix will be nrows M x nrows M
00713 //==============================================================================
00714 template <class T, unsigned int D1, unsigned int D2, class R>
00715 inline SMatrix<T,D1,D1,MatRepSym<T,D1> > Similarity(const SMatrix<T,D1,D2,R>& lhs, const SMatrix<T,D2,D2,MatRepSym<T,D2> >& rhs) {
00716   SMatrix<T,D1,D2, MatRepStd<T,D1,D2> > tmp = lhs * rhs;
00717   typedef  SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym; 
00718   SMatrixSym mret; 
00719   AssignSym::Evaluate(mret,  tmp * Transpose(lhs)  ); 
00720   return mret; 
00721 }
00722 
00723 //==============================================================================
00724 // product: SMatrix/SMatrix calculate M * A * M^T where A is a symmetric matrix
00725 // return matrix will be nrowsM x nrows M
00726 // M is a matrix expression
00727 //==============================================================================
00728 template <class A, class T, unsigned int D1, unsigned int D2, class R>
00729 inline SMatrix<T,D1,D1,MatRepSym<T,D1> > Similarity(const Expr<A,T,D1,D2,R>& lhs, const SMatrix<T,D2,D2,MatRepSym<T,D2> >& rhs) {
00730   SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = lhs * rhs;
00731   typedef  SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym; 
00732   SMatrixSym mret; 
00733   AssignSym::Evaluate(mret,  tmp * Transpose(lhs)  ); 
00734   return mret; 
00735 }
00736 
00737 #ifdef XXX
00738     // not needed (
00739 //==============================================================================
00740 // product: SMatrix/SMatrix calculate M * A * M where A and M are symmetric matrices
00741 // return matrix will be nrows M x nrows M
00742 //==============================================================================
00743 template <class T, unsigned int D1>
00744 inline SMatrix<T,D1,D1,MatRepSym<T,D1> > Similarity(const SMatrix<T,D1,D1,MatRepSym<T,D1> >& lhs, const SMatrix<T,D1,D1,MatRepSym<T,D1> >& rhs) {
00745   SMatrix<T,D1,D1, MatRepStd<T,D1,D1> > tmp = lhs * rhs;
00746   typedef  SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym; 
00747   SMatrixSym mret; 
00748   AssignSym::Evaluate(mret,  tmp * lhs  ); 
00749   return mret; 
00750 }
00751 #endif
00752 
00753 
00754 /**
00755    Transpose Similarity Matrix Product :  B = U^T * A * U for A symmetric
00756    returning a symmetric matrix expression: \f$ B(i,j) = \sum_{k,l} U(k,i) * A(k,l) * U(l,j) \f$ 
00757 
00758    @ingroup MatrixFunctions
00759  */
00760 //==============================================================================
00761 // product: SMatrix/SMatrix calculate M^T * A * M where A is a symmetric matrix
00762 // return matrix will be ncolsM x ncols M
00763 //==============================================================================
00764 template <class T, unsigned int D1, unsigned int D2, class R>
00765 inline SMatrix<T,D2,D2,MatRepSym<T,D2> > SimilarityT(const SMatrix<T,D1,D2,R>& lhs, const SMatrix<T,D1,D1,MatRepSym<T,D1> >& rhs) {
00766   SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = rhs * lhs;
00767   typedef  SMatrix<T,D2,D2,MatRepSym<T,D2> > SMatrixSym; 
00768   SMatrixSym mret; 
00769   AssignSym::Evaluate(mret,  Transpose(lhs) * tmp ); 
00770   return mret; 
00771 }
00772 
00773 //==============================================================================
00774 // product: SMatrix/SMatrix calculate M^T * A * M where A is a symmetric matrix
00775 // return matrix will be ncolsM x ncols M
00776 // M is a matrix expression
00777 //==============================================================================
00778 template <class A, class T, unsigned int D1, unsigned int D2, class R>
00779 inline SMatrix<T,D2,D2,MatRepSym<T,D2> > SimilarityT(const Expr<A,T,D1,D2,R>& lhs, const SMatrix<T,D1,D1,MatRepSym<T,D1> >& rhs) {
00780   SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = rhs * lhs;
00781   typedef  SMatrix<T,D2,D2,MatRepSym<T,D2> > SMatrixSym; 
00782   SMatrixSym mret; 
00783   AssignSym::Evaluate(mret,  Transpose(lhs) * tmp ); 
00784   return mret; 
00785 }
00786 
00787 
00788 
00789 
00790 
00791 // //==============================================================================
00792 // // Mult * (Expr * Expr, binary) with a symmetric result
00793 // // the operation is done only for half
00794 // //==============================================================================
00795 // template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
00796 // inline Expr<MatrixMulOp<Expr<A,T,D,D,MatRepSym<T,D> >, Expr<B,T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
00797 //  operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
00798 //   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, T,D> MatMulOp;
00799 //   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
00800 // }
00801 
00802 
00803 
00804 //==============================================================================
00805 // TensorMulOp
00806 //==============================================================================
00807 /**
00808    Class for Tensor Multiplication (outer product) of two vectors   
00809    giving a matrix 
00810 
00811    @ingroup Expression
00812  */
00813 template <class Vector1, class Vector2>
00814 class TensorMulOp {
00815 public:
00816   ///
00817   TensorMulOp( const Vector1 & lhs, const Vector2 & rhs) :
00818     lhs_(lhs),
00819     rhs_(rhs) {}
00820 
00821   ///
00822   ~TensorMulOp() {}
00823 
00824   /// Vector2::kSize is the number of columns in the resulting matrix
00825   inline typename Vector1::value_type apply(unsigned int i) const {
00826     return lhs_.apply( i/ Vector2::kSize) * rhs_.apply( i % Vector2::kSize );
00827   }
00828   inline typename Vector1::value_type operator() (unsigned int i, unsigned j) const {
00829     return lhs_.apply(i) * rhs_.apply(j);
00830   }
00831 
00832   inline bool IsInUse (const typename Vector1::value_type * ) const { 
00833     return false; 
00834   }
00835 
00836 
00837 protected:
00838 
00839   const Vector1 & lhs_;
00840   const Vector2 & rhs_;
00841 
00842 };
00843 
00844 
00845 
00846 /**
00847    Tensor Vector Product : M(i,j) = v(i) * v(j) 
00848    returning a matrix expression 
00849 
00850    @ingroup VectFunction
00851  */
00852 
00853 #ifndef _WIN32
00854 
00855     // Tensor Prod (use default MatRepStd for the returned expression
00856     // cannot make a symmetric matrix
00857 //==============================================================================
00858 // TensorProd (SVector x SVector)
00859 //==============================================================================
00860 template <class T, unsigned int D1, unsigned int D2>
00861 inline Expr<TensorMulOp<SVector<T,D1>, SVector<T,D2>  >, T, D1, D2 >
00862  TensorProd(const SVector<T,D1>& lhs, const SVector<T,D2>& rhs) {
00863   typedef TensorMulOp<SVector<T,D1>, SVector<T,D2> > TVMulOp;
00864   return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
00865 }
00866 
00867 //==============================================================================
00868 // TensorProd (VecExpr x SVector)
00869 //==============================================================================
00870  template <class T, unsigned int D1, unsigned int D2, class A>
00871 inline Expr<TensorMulOp<VecExpr<A,T,D1>, SVector<T,D2>  >, T, D1, D2 >
00872  TensorProd(const VecExpr<A,T,D1>& lhs, const SVector<T,D2>& rhs) {
00873   typedef TensorMulOp<VecExpr<A,T,D1>, SVector<T,D2> > TVMulOp;
00874   return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
00875 }
00876 
00877 //==============================================================================
00878 // TensorProd (SVector x VecExpr)
00879 //==============================================================================
00880  template <class T, unsigned int D1, unsigned int D2, class A>
00881 inline Expr<TensorMulOp<SVector<T,D1>, VecExpr<A,T,D2>  >, T, D1, D2 >
00882  TensorProd(const SVector<T,D1>& lhs, const VecExpr<A,T,D2>& rhs) {
00883   typedef TensorMulOp<SVector<T,D1>, VecExpr<A,T,D2> > TVMulOp;
00884   return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
00885 }
00886 
00887 
00888 //==============================================================================
00889 // TensorProd (VecExpr x VecExpr)
00890 //==============================================================================
00891  template <class T, unsigned int D1, unsigned int D2, class A, class B>
00892 inline Expr<TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2>  >, T, D1, D2 >
00893  TensorProd(const VecExpr<A,T,D1>& lhs, const VecExpr<B,T,D2>& rhs) {
00894   typedef TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2> > TVMulOp;
00895   return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
00896 }
00897 
00898 #endif
00899 #ifdef _WIN32 
00900 /// case of WINDOWS - problem using Expression (  C1001: INTERNAL COMPILER ERROR )
00901 
00902 //==============================================================================
00903 // TensorProd (SVector x SVector)
00904 //==============================================================================
00905 template <class T, unsigned int D1, unsigned int D2>
00906 inline SMatrix<T,D1,D2>  TensorProd(const SVector<T,D1>& lhs, const SVector<T,D2>& rhs) {
00907   SMatrix<T,D1,D2> tmp;
00908   for (unsigned int i=0; i< D1; ++i) 
00909     for (unsigned int j=0; j< D2; ++j) {
00910       tmp(i,j) = lhs[i]*rhs[j];
00911     }
00912 
00913   return tmp;
00914 }
00915 //==============================================================================
00916 // TensorProd (VecExpr x SVector)
00917 //==============================================================================
00918  template <class T, unsigned int D1, unsigned int D2, class A>
00919 inline SMatrix<T,D1,D2>  TensorProd(const VecExpr<A,T,D1>& lhs, const SVector<T,D2>& rhs) {
00920   SMatrix<T,D1,D2> tmp;
00921   for (unsigned int i=0; i< D1; ++i)
00922     for (unsigned int j=0; j< D2; ++j)
00923       tmp(i,j) = lhs.apply(i) * rhs.apply(j);
00924 
00925   return tmp;
00926 }
00927 //==============================================================================
00928 // TensorProd (SVector x VecExpr)
00929 //==============================================================================
00930  template <class T, unsigned int D1, unsigned int D2, class A>
00931 inline SMatrix<T,D1,D2> TensorProd(const SVector<T,D1>& lhs, const VecExpr<A,T,D2>& rhs) {
00932   SMatrix<T,D1,D2> tmp;
00933   for (unsigned int i=0; i< D1; ++i)
00934     for (unsigned int j=0; j< D2; ++j)
00935       tmp(i,j) = lhs.apply(i) * rhs.apply(j);
00936 
00937   return tmp;
00938 }
00939 
00940 //==============================================================================
00941 // TensorProd (VecExpr x VecExpr)
00942 //==============================================================================
00943 
00944  template <class T, unsigned int D1, unsigned int D2, class A, class B>
00945 inline SMatrix<T,D1,D2  > TensorProd(const VecExpr<A,T,D1>& lhs, const VecExpr<B,T,D2>& rhs) {
00946   SMatrix<T,D1,D2> tmp;
00947   for (unsigned int i=0; i< D1; ++i)
00948     for (unsigned int j=0; j< D2; ++j)
00949       tmp(i,j) = lhs.apply(i) * rhs.apply(j);
00950 
00951   return tmp;
00952 }
00953 
00954 
00955 #endif
00956 
00957 // solving a positive defined symmetric linear system using Choleski decompositions
00958 // matrix will be decomposed and the returned vector will be overwritten in vec
00959 // If the user wants to pass const objects need to copy the matrices
00960 // It will work only for symmetric matrices
00961 template <class T, unsigned int D>
00962 bool SolveChol( SMatrix<T, D, D, MatRepSym<T, D>  > & mat,  SVector<T, D> & vec ) { 
00963    CholeskyDecomp<T, D> decomp(mat);
00964    return decomp.Solve(vec); 
00965 }
00966 
00967 /// same function as before but not overwriting the matrix and returning a copy of the vector
00968 /// (this is the slow version)
00969 template <class T, unsigned int D>
00970 SVector<T,D> SolveChol( const SMatrix<T, D, D, MatRepSym<T, D>  > & mat,  const SVector<T, D> & vec, int & ifail  ) { 
00971    SMatrix<T, D, D, MatRepSym<T, D> > atmp(mat); 
00972    SVector<T,D> vret(vec); 
00973    bool ok = SolveChol( atmp, vret); 
00974    ifail =  (ok) ? 0 : -1;
00975    return vret; 
00976 }
00977 
00978 
00979 
00980   }  // namespace Math
00981 
00982 }  // namespace ROOT
00983           
00984 
00985 #endif  /* ROOT_Math_MatrixFunctions */