root/Math/MatrixFunctions.h

0001 // @(#)root/smatrix:$Id$
0002 // Authors: T. Glebe, L. Moneta    2005
0003
0004 #ifndef ROOT_Math_MatrixFunctions
0005 #define ROOT_Math_MatrixFunctions
0006 // ********************************************************************
0007 //
0008 // source:
0009 //
0010 // type:      source code
0011 //
0012 // created:   20. Mar 2001
0013 //
0014 // author:    Thorsten Glebe
0015 //            HERA-B Collaboration
0016 //            Max-Planck-Institut fuer Kernphysik
0017 //            Saupfercheckweg 1
0018 //            69117 Heidelberg
0019 //            Germany
0020 //            E-mail: T.Glebe@mpi-hd.mpg.de
0021 //
0022 // Description: Functions/Operators special to Matrix
0023 //
0024 // changes:
0025 // 20 Mar 2001 (TG) creation
0026 // 20 Mar 2001 (TG) Added Matrix * Vector multiplication
0027 // 21 Mar 2001 (TG) Added transpose, product
0028 // 11 Apr 2001 (TG) transpose() speed improvement by removing rows(), cols()
0029 //                  through static members of Matrix and Expr class
0030 //
0031 // ********************************************************************
0032
0033
0034 /**
0035 \defgroup MatrixFunctions Matrix Template Functions
0036 \ingroup SMatrixGroup
0037
0038 These function apply to matrices (and also Matrix expression) and can return a
0039 matrix expression of a particular defined type, like in the matrix multiplication or
0040 a vector, like in the matrix-vector product or a scalar like in the Similarity
0041 vector-matrix product.
0042
0043 */
0044
0045 #include "Math/BinaryOpPolicy.h"
0046 #include "Math/Expression.h"
0047 #include "Math/HelperOps.h"
0048 #include "Math/CholeskyDecomp.h"
0049
0050 namespace ROOT {
0051
0052   namespace Math {
0053
0054     template <class T, unsigned int D>
0055     class SVector;
0056
0057 #ifdef XXX
0058 //==============================================================================
0059 // SMatrix * SVector
0060 //==============================================================================
0061 template <class T, unsigned int D1, unsigned int D2, class R>
0062 SVector<T,D1> operator*(const SMatrix<T,D1,D2,R>& rhs, const SVector<T,D2>& lhs)
0063 {
0064   SVector<T,D1> tmp;
0065   for(unsigned int i=0; i<D1; ++i) {
0066     const unsigned int rpos = i*D2;
0067     for(unsigned int j=0; j<D2; ++j) {
0068       tmp[i] += rhs.apply(rpos+j) * lhs.apply(j);
0069     }
0070   }
0071   return tmp;
0072 }
0073 #endif
0074
0075
0076 // matrix-vector product:
0077 // use apply(i) function for matrices. Tested  (11/05/06) with using (i,j) but
0078 // performances are slightly worse (not clear  why)
0079
0080 //==============================================================================
0081 // meta_row_dot
0082 //==============================================================================
0083 template <unsigned int I>
0084 struct meta_row_dot {
0085   template <class A, class B>
0086   static inline typename A::value_type f(const A& lhs, const B& rhs,
0087                                          const unsigned int offset) {
0088     return lhs.apply(offset+I) * rhs.apply(I) + meta_row_dot<I-1>::f(lhs,rhs,offset);
0089   }
0090 };
0091
0092
0093 //==============================================================================
0094 // meta_row_dot<0>
0095 //==============================================================================
0096 template <>
0097 struct meta_row_dot<0> {
0098   template <class A, class B>
0099   static inline typename A::value_type f(const A& lhs, const B& rhs,
0100                                          const unsigned int offset) {
0101     return lhs.apply(offset) * rhs.apply(0);
0102   }
0103 };
0104
0105 //==============================================================================
0106 // VectorMatrixRowOp
0107 //==============================================================================
0108 template <class Matrix, class Vector, unsigned int D2>
0109 class VectorMatrixRowOp {
0110 public:
0111
0112   typedef typename Vector::value_type T;
0113
0114   ///
0115   VectorMatrixRowOp(const Matrix& lhs, const Vector& rhs) :
0116     lhs_(lhs), rhs_(rhs) {}
0117
0118   ///
0119   ~VectorMatrixRowOp() {}
0120
0121   /// calc \f$ \sum_{j} a_{ij} * v_j \f$
0122   inline typename Matrix::value_type apply(unsigned int i) const {
0123     return meta_row_dot<D2-1>::f(lhs_, rhs_, i*D2);
0124   }
0125
0126    // check if passed pointer is in use
0127    // check only the vector since this is a vector expression
0128   inline bool IsInUse (const T * p) const {
0129     return rhs_.IsInUse(p);
0130   }
0131
0132
0133 protected:
0134   const Matrix& lhs_;
0135   const Vector& rhs_;
0136 };
0137
0138
0139 //==============================================================================
0140 // meta_col_dot
0141 //==============================================================================
0142 template <unsigned int I>
0143 struct meta_col_dot {
0144   template <class Matrix, class Vector>
0145   static inline typename Matrix::value_type f(const Matrix& lhs, const Vector& rhs,
0146                                               const unsigned int offset) {
0147     return lhs.apply(Matrix::kCols*I+offset) * rhs.apply(I) +
0148            meta_col_dot<I-1>::f(lhs,rhs,offset);
0149   }
0150 };
0151
0152
0153 //==============================================================================
0154 // meta_col_dot<0>
0155 //==============================================================================
0156 template <>
0157 struct meta_col_dot<0> {
0158   template <class Matrix, class Vector>
0159   static inline typename Matrix::value_type f(const Matrix& lhs, const Vector& rhs,
0160                                               const unsigned int offset) {
0161     return lhs.apply(offset) * rhs.apply(0);
0162   }
0163 };
0164
0165 //==============================================================================
0166 // VectorMatrixColOp
0167 //==============================================================================
0168 /**
0169    Class for Vector-Matrix multiplication
0170
0171    @ingroup Expression
0172  */
0173 template <class Vector, class Matrix, unsigned int D1>
0174 class VectorMatrixColOp {
0175 public:
0176
0177   typedef typename Vector::value_type T;
0178   ///
0179   VectorMatrixColOp(const Vector& lhs, const Matrix& rhs) :
0180     lhs_(lhs), rhs_(rhs) {}
0181
0182   ///
0183   ~VectorMatrixColOp() {}
0184
0185   /// calc \f$ \sum_{j} a_{ij} * v_j \f$
0186   inline typename Matrix::value_type apply(unsigned int i) const {
0187     return meta_col_dot<D1-1>::f(rhs_, lhs_, i);
0188   }
0189
0190    // check if passed pointer is in use
0191    // check only the vector since this is a vector expression
0192    inline bool IsInUse (const T * p) const {
0193     return lhs_.IsInUse(p);
0194   }
0195
0196
0197 protected:
0198   const Vector&    lhs_;
0199   const Matrix&    rhs_;
0200 };
0201
0202 /**
0203    Matrix *  Vector multiplication   \f$ a(i) = \sum_{j} M(i,j) * b(j) \f$
0204    returning a vector expression
0205
0206    @ingroup MatrixFunctions
0207  */
0208 //==============================================================================
0209 // operator*: SMatrix * SVector
0210 //==============================================================================
0211 template <class T, unsigned int D1, unsigned int D2, class R>
0212 inline VecExpr<VectorMatrixRowOp<SMatrix<T,D1,D2,R>,SVector<T,D2>, D2>, T, D1>
0213  operator*(const SMatrix<T,D1,D2,R>& lhs, const SVector<T,D2>& rhs) {
0214   typedef VectorMatrixRowOp<SMatrix<T,D1,D2,R>,SVector<T,D2>, D2> VMOp;
0215   return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
0216 }
0217
0218 //==============================================================================
0219 // operator*: SMatrix * Expr<A,T,D2>
0220 //==============================================================================
0221 template <class A, class T, unsigned int D1, unsigned int D2, class R>
0222 inline VecExpr<VectorMatrixRowOp<SMatrix<T,D1,D2,R>, VecExpr<A,T,D2>, D2>, T, D1>
0223  operator*(const SMatrix<T,D1,D2,R>& lhs, const VecExpr<A,T,D2>& rhs) {
0224   typedef VectorMatrixRowOp<SMatrix<T,D1,D2,R>,VecExpr<A,T,D2>, D2> VMOp;
0225   return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
0226 }
0227
0228 //==============================================================================
0229 // operator*: Expr<A,T,D1,D2> * SVector
0230 //==============================================================================
0231 template <class A, class T, unsigned int D1, unsigned int D2, class R>
0232 inline VecExpr<VectorMatrixRowOp<Expr<A,T,D1,D2,R>, SVector<T,D2>, D2>, T, D1>
0233  operator*(const Expr<A,T,D1,D2,R>& lhs, const SVector<T,D2>& rhs) {
0234   typedef VectorMatrixRowOp<Expr<A,T,D1,D2,R>,SVector<T,D2>, D2> VMOp;
0235   return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
0236 }
0237
0238 //==============================================================================
0239 // operator*: Expr<A,T,D1,D2> * VecExpr<B,T,D2>
0240 //==============================================================================
0241 template <class A, class B, class T, unsigned int D1, unsigned int D2, class R>
0242 inline VecExpr<VectorMatrixRowOp<Expr<A,T,D1,D2,R>, VecExpr<B,T,D2>, D2>, T, D1>
0243  operator*(const Expr<A,T,D1,D2,R>& lhs, const VecExpr<B,T,D2>& rhs) {
0244   typedef VectorMatrixRowOp<Expr<A,T,D1,D2,R>,VecExpr<B,T,D2>, D2> VMOp;
0245   return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
0246 }
0247
0248 //==============================================================================
0249 // operator*: SVector * SMatrix
0250 //==============================================================================
0251 template <class T, unsigned int D1, unsigned int D2, class R>
0252 inline VecExpr<VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2,R>, D1>, T, D2>
0253  operator*(const SVector<T,D1>& lhs, const SMatrix<T,D1,D2,R>& rhs) {
0254   typedef VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2,R>, D1> VMOp;
0255   return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
0256 }
0257
0258 //==============================================================================
0259 // operator*: SVector * Expr<A,T,D1,D2>
0260 //==============================================================================
0261 template <class A, class T, unsigned int D1, unsigned int D2, class R>
0262 inline VecExpr<VectorMatrixColOp<SVector<T,D1>, Expr<A,T,D1,D2,R>, D1>, T, D2>
0263  operator*(const SVector<T,D1>& lhs, const Expr<A,T,D1,D2,R>& rhs) {
0264   typedef VectorMatrixColOp<SVector<T,D1>, Expr<A,T,D1,D2,R>, D1> VMOp;
0265   return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
0266 }
0267
0268 //==============================================================================
0269 // operator*: VecExpr<A,T,D1> * SMatrix
0270 //==============================================================================
0271 template <class A, class T, unsigned int D1, unsigned int D2, class R>
0272 inline VecExpr<VectorMatrixColOp<VecExpr<A,T,D1>, SMatrix<T,D1,D2,R>, D1>, T, D2>
0273  operator*(const VecExpr<A,T,D1>& lhs, const SMatrix<T,D1,D2,R>& rhs) {
0274   typedef VectorMatrixColOp<VecExpr<A,T,D1>, SMatrix<T,D1,D2,R>, D1> VMOp;
0275   return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
0276 }
0277
0278 //==============================================================================
0279 // operator*: VecExpr<A,T,D1> * Expr<B,T,D1,D2>
0280 //==============================================================================
0281 template <class A, class B, class T, unsigned int D1, unsigned int D2, class R>
0282 inline VecExpr<VectorMatrixColOp<VecExpr<A,T,D1>, Expr<B,T,D1,D2,R>, D1>, T, D2>
0283  operator*(const VecExpr<A,T,D1>& lhs, const Expr<B,T,D1,D2,R>& rhs) {
0284   typedef VectorMatrixColOp<VecExpr<A,T,D1>, Expr<B,T,D1,D2,R>, D1> VMOp;
0285   return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
0286 }
0287
0288 //==============================================================================
0289 // meta_matrix_dot
0290 //==============================================================================
0291 template <unsigned int I>
0292 struct meta_matrix_dot {
0293
0294   template <class MatrixA, class MatrixB>
0295   static inline typename MatrixA::value_type f(const MatrixA& lhs,
0296                                                const MatrixB& rhs,
0297                                                const unsigned int offset) {
0298     return lhs.apply(offset/MatrixB::kCols*MatrixA::kCols + I) *
0299            rhs.apply(MatrixB::kCols*I + offset%MatrixB::kCols) +
0300            meta_matrix_dot<I-1>::f(lhs,rhs,offset);
0301   }
0302
0303   // multiplication using i and j indices
0304   template <class MatrixA, class MatrixB>
0305   static inline typename MatrixA::value_type g(const MatrixA& lhs,
0306                                                const MatrixB& rhs,
0307                                                unsigned int i,
0308                                                unsigned int j) {
0309     return lhs(i, I) * rhs(I , j) +
0310            meta_matrix_dot<I-1>::g(lhs,rhs,i,j);
0311   }
0312 };
0313
0314
0315 //==============================================================================
0316 // meta_matrix_dot<0>
0317 //==============================================================================
0318 template <>
0319 struct meta_matrix_dot<0> {
0320
0321   template <class MatrixA, class MatrixB>
0322   static inline typename MatrixA::value_type f(const MatrixA& lhs,
0323                                                const MatrixB& rhs,
0324                                                const unsigned int offset) {
0325     return lhs.apply(offset/MatrixB::kCols*MatrixA::kCols) *
0326            rhs.apply(offset%MatrixB::kCols);
0327   }
0328
0329   // multiplication using i and j
0330   template <class MatrixA, class MatrixB>
0331   static inline typename MatrixA::value_type g(const MatrixA& lhs,
0332                                                const MatrixB& rhs,
0333                                                unsigned int i, unsigned int j) {
0334     return lhs(i,0) * rhs(0,j);
0335   }
0336
0337 };
0338
0339 //==============================================================================
0340 // MatrixMulOp
0341 //==============================================================================
0342 /**
0343    Class for Matrix-Matrix multiplication
0344
0345    @ingroup Expression
0346  */
0347 template <class MatrixA, class MatrixB, class T, unsigned int D>
0348 class MatrixMulOp {
0349 public:
0350   ///
0351   MatrixMulOp(const MatrixA& lhs, const MatrixB& rhs) :
0352     lhs_(lhs), rhs_(rhs) {}
0353
0354   ///
0355   ~MatrixMulOp() {}
0356
0357   /// calc \f$\sum_{j} a_{ik} * b_{kj}\f$
0358   inline T apply(unsigned int i) const {
0359     return meta_matrix_dot<D-1>::f(lhs_, rhs_, i);
0360   }
0361
0362   inline T operator() (unsigned int i, unsigned j) const {
0363     return meta_matrix_dot<D-1>::g(lhs_, rhs_, i, j);
0364   }
0365
0366   inline bool IsInUse (const T * p) const {
0367     return lhs_.IsInUse(p) || rhs_.IsInUse(p);
0368   }
0369
0370
0371 protected:
0372   const MatrixA&    lhs_;
0373   const MatrixB&    rhs_;
0374 };
0375
0376
0377 /**
0378    Matrix *  Matrix multiplication , \f$ C(i,j) = \sum_{k} A(i,k) * B(k,j)\f$
0379    returning a matrix expression
0380
0381    @ingroup MatrixFunctions
0382  */
0383 //==============================================================================
0384 // operator* (SMatrix * SMatrix, binary)
0385 //==============================================================================
0386 template <  class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
0387 inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
0388  operator*(const SMatrix<T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
0389   typedef MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, T,D> MatMulOp;
0390   return Expr<MatMulOp,T,D1,D2,
0391     typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
0392 }
0393
0394 //==============================================================================
0395 // operator* (SMatrix * Expr, binary)
0396 //==============================================================================
0397 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
0398 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>
0399  operator*(const SMatrix<T,D1,D,R1>& lhs, const Expr<A,T,D,D2,R2>& rhs) {
0400   typedef MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>,T,D> MatMulOp;
0401   return Expr<MatMulOp,T,D1,D2,
0402     typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
0403 }
0404
0405 //==============================================================================
0406 // operator* (Expr * SMatrix, binary)
0407 //==============================================================================
0408 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
0409 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>
0410  operator*(const Expr<A,T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
0411   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>,T,D> MatMulOp;
0412   return Expr<MatMulOp,T,D1,D2,
0413     typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
0414 }
0415
0416 //==============================================================================
0417 // operator* (Expr * Expr, binary)
0418 //==============================================================================
0419 template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
0420 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>
0421  operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
0422   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, T,D> MatMulOp;
0423   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
0424 }
0425
0426
0427
0428 #ifdef XXX
0429 //==============================================================================
0430 // MatrixMulOp
0431 //==============================================================================
0432 template <class MatrixA, class MatrixB, unsigned int D>
0433 class MatrixMulOp {
0434 public:
0435   ///
0436   MatrixMulOp(const MatrixA& lhs, const MatrixB& rhs) :
0437     lhs_(lhs), rhs_(rhs) {}
0438
0439   ///
0440   ~MatrixMulOp() {}
0441
0442   /// calc \f$\sum_{j} a_{ik} * b_{kj}\f$
0443   inline typename MatrixA::value_type apply(unsigned int i) const {
0444     return meta_matrix_dot<D-1>::f(lhs_, rhs_, i);
0445   }
0446
0447 protected:
0448   const MatrixA&    lhs_;
0449   const MatrixB&    rhs_;
0450 };
0451
0452
0453 //==============================================================================
0454 // operator* (SMatrix * SMatrix, binary)
0455 //==============================================================================
0456 template <  class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
0457 inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
0458  operator*(const SMatrix<T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
0459   typedef MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, D> MatMulOp;
0460   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
0461 }
0462
0463 //==============================================================================
0464 // operator* (SMatrix * Expr, binary)
0465 //==============================================================================
0466 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
0467 inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
0468  operator*(const SMatrix<T,D1,D,R1>& lhs, const Expr<A,T,D,D2,R2>& rhs) {
0469   typedef MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>, D> MatMulOp;
0470   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
0471 }
0472
0473 //==============================================================================
0474 // operator* (Expr * SMatrix, binary)
0475 //==============================================================================
0476 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
0477 inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
0478  operator*(const Expr<A,T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
0479   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>, D> MatMulOp;
0480   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
0481 }
0482
0483 //=============================================================================
0484 // operator* (Expr * Expr, binary)
0485 //=============================================================================
0486 template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
0487 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>
0488  operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
0489   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, D> MatMulOp;
0490   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
0491 }
0492 #endif
0493
0494 //==============================================================================
0495 // TransposeOp
0496 //==============================================================================
0497 /**
0498    Class for Transpose Operations
0499
0500    @ingroup Expression
0501  */
0502 template <class Matrix, class T, unsigned int D1, unsigned int D2=D1>
0503 class TransposeOp {
0504 public:
0505   ///
0506   TransposeOp( const Matrix& rhs) :
0507     rhs_(rhs) {}
0508
0509   ///
0510   ~TransposeOp() {}
0511
0512   ///
0513   inline T apply(unsigned int i) const {
0514     return rhs_.apply( (i%D1)*D2 + i/D1);
0515   }
0516   inline T operator() (unsigned int i, unsigned j) const {
0517     return rhs_( j, i);
0518   }
0519
0520   inline bool IsInUse (const T * p) const {
0521     return rhs_.IsInUse(p);
0522   }
0523
0524 protected:
0525   const Matrix& rhs_;
0526 };
0527
0528
0529 /**
0530    Matrix Transpose   B(i,j) = A(j,i)
0531    returning a matrix expression
0532
0533    @ingroup MatrixFunctions
0534  */
0535 //==============================================================================
0536 // transpose
0537 //==============================================================================
0538 template <class T, unsigned int D1, unsigned int D2, class R>
0539 inline Expr<TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2>, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
0540  Transpose(const SMatrix<T,D1,D2, R>& rhs) {
0541   typedef TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2> MatTrOp;
0542
0543   return Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs));
0544 }
0545
0546 //==============================================================================
0547 // transpose
0548 //==============================================================================
0549 template <class A, class T, unsigned int D1, unsigned int D2, class R>
0550 inline Expr<TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2>, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
0551  Transpose(const Expr<A,T,D1,D2,R>& rhs) {
0552   typedef TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2> MatTrOp;
0553
0554   return Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs));
0555 }
0556
0557
0558 #ifdef ENABLE_TEMPORARIES_TRANSPOSE
0559 // sometimes is faster to create a temp, not clear why
0560
0561 //==============================================================================
0562 // transpose
0563 //==============================================================================
0564 template <class T, unsigned int D1, unsigned int D2, class R>
0565 inline SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
0566  Transpose(const SMatrix<T,D1,D2, R>& rhs) {
0567   typedef TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2> MatTrOp;
0568
0569   return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
0570     ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
0571 }
0572
0573 //==============================================================================
0574 // transpose
0575 //==============================================================================
0576 template <class A, class T, unsigned int D1, unsigned int D2, class R>
0577 inline  SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
0578  Transpose(const Expr<A,T,D1,D2,R>& rhs) {
0579   typedef TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2> MatTrOp;
0580
0581   return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
0582     ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
0583 }
0584
0585 #endif
0586
0587
0588 #ifdef OLD
0589 //==============================================================================
0590 // product: SMatrix/SVector calculate v^T * A * v
0591 //==============================================================================
0592 template <class T, unsigned int D, class R>
0593 inline T Product(const SMatrix<T,D,D,R>& lhs, const SVector<T,D>& rhs) {
0594   return Dot(rhs, lhs * rhs);
0595 }
0596
0597 //==============================================================================
0598 // product: SVector/SMatrix calculate v^T * A * v
0599 //==============================================================================
0600 template <class T, unsigned int D, class R>
0601 inline T Product(const SVector<T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
0602   return Dot(lhs, rhs * lhs);
0603 }
0604
0605 //==============================================================================
0606 // product: SMatrix/Expr calculate v^T * A * v
0607 //==============================================================================
0608 template <class A, class T, unsigned int D, class R>
0609 inline T Product(const SMatrix<T,D,D,R>& lhs, const VecExpr<A,T,D>& rhs) {
0610   return Dot(rhs, lhs * rhs);
0611 }
0612
0613 //==============================================================================
0614 // product: Expr/SMatrix calculate v^T * A * v
0615 //==============================================================================
0616 template <class A, class T, unsigned int D, class R>
0617 inline T Product(const VecExpr<A,T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
0618   return Dot(lhs, rhs * lhs);
0619 }
0620
0621 //==============================================================================
0622 // product: SVector/Expr calculate v^T * A * v
0623 //==============================================================================
0624 template <class A, class T, unsigned int D, class R>
0625 inline T Product(const SVector<T,D>& lhs, const Expr<A,T,D,D,R>& rhs) {
0626   return Dot(lhs, rhs * lhs);
0627 }
0628
0629 //==============================================================================
0630 // product: Expr/SVector calculate v^T * A * v
0631 //==============================================================================
0632 template <class A, class T, unsigned int D, class R>
0633 inline T Product(const Expr<A,T,D,D,R>& lhs, const SVector<T,D>& rhs) {
0634   return Dot(rhs, lhs * rhs);
0635 }
0636
0637 //==============================================================================
0638 // product: Expr/Expr calculate v^T * A * v
0639 //==============================================================================
0640 template <class A, class B, class T, unsigned int D, class R>
0641 inline T Product(const Expr<A,T,D,D,R>& lhs, const VecExpr<B,T,D>& rhs) {
0642   return Dot(rhs, lhs * rhs);
0643 }
0644
0645 //==============================================================================
0646 // product: Expr/Expr calculate v^T * A * v
0647 //==============================================================================
0648 template <class A, class B, class T, unsigned int D, class R>
0649 inline T Product(const VecExpr<A,T,D>& lhs, const Expr<B,T,D,D,R>& rhs) {
0650   return Dot(lhs, rhs * lhs);
0651 }
0652 #endif
0653
0654 /**
0655    Similarity Vector - Matrix Product:  v^T * A * v
0656    returning a scalar value of type T   \f$ s = \sum_{i,j} v(i) * A(i,j) * v(j)\f$
0657
0658    @ingroup MatrixFunctions
0659  */
0660
0661 //==============================================================================
0662 // product: SMatrix/SVector calculate v^T * A * v
0663 //==============================================================================
0664 template <class T, unsigned int D, class R>
0665 inline T Similarity(const SMatrix<T,D,D,R>& lhs, const SVector<T,D>& rhs) {
0666   return Dot(rhs, lhs * rhs);
0667 }
0668
0669 //==============================================================================
0670 // product: SVector/SMatrix calculate v^T * A * v
0671 //==============================================================================
0672 template <class T, unsigned int D, class R>
0673 inline T Similarity(const SVector<T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
0674   return Dot(lhs, rhs * lhs);
0675 }
0676
0677 //==============================================================================
0678 // product: SMatrix/Expr calculate v^T * A * v
0679 //==============================================================================
0680 template <class A, class T, unsigned int D, class R>
0681 inline T Similarity(const SMatrix<T,D,D,R>& lhs, const VecExpr<A,T,D>& rhs) {
0682   return Dot(rhs, lhs * rhs);
0683 }
0684
0685 //==============================================================================
0686 // product: Expr/SMatrix calculate v^T * A * v
0687 //==============================================================================
0688 template <class A, class T, unsigned int D, class R>
0689 inline T Similarity(const VecExpr<A,T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
0690   return Dot(lhs, rhs * lhs);
0691 }
0692
0693 //==============================================================================
0694 // product: SVector/Expr calculate v^T * A * v
0695 //==============================================================================
0696 template <class A, class T, unsigned int D, class R>
0697 inline T Similarity(const SVector<T,D>& lhs, const Expr<A,T,D,D,R>& rhs) {
0698   return Dot(lhs, rhs * lhs);
0699 }
0700
0701 //==============================================================================
0702 // product: Expr/SVector calculate v^T * A * v
0703 //==============================================================================
0704 template <class A, class T, unsigned int D, class R>
0705 inline T Similarity(const Expr<A,T,D,D,R>& lhs, const SVector<T,D>& rhs) {
0706   return Dot(rhs, lhs * rhs);
0707 }
0708
0709 //==============================================================================
0710 // product: Expr/Expr calculate v^T * A * v
0711 //==============================================================================
0712 template <class A, class B, class T, unsigned int D, class R>
0713 inline T Similarity(const Expr<A,T,D,D,R>& lhs, const VecExpr<B,T,D>& rhs) {
0714   return Dot(rhs, lhs * rhs);
0715 }
0716
0717 //==============================================================================
0718 // product: Expr/Expr calculate v^T * A * v
0719 //==============================================================================
0720 template <class A, class B, class T, unsigned int D, class R>
0721 inline T Similarity(const VecExpr<A,T,D>& lhs, const Expr<B,T,D,D,R>& rhs) {
0722   return Dot(lhs, rhs * lhs);
0723 }
0724
0725
0726 /**
0727    Similarity Matrix Product :  B = U * A * U^T for A symmetric
0728    returning a symmetric matrix expression:
0729    \f$ B(i,j) = \sum_{k,l} U(i,k) * A(k,l) * U(j,l) \f$
0730
0731    @ingroup MatrixFunctions
0732  */
0733 //==============================================================================
0734 // product: SMatrix/SMatrix calculate M * A * M^T where A is a symmetric matrix
0735 // return matrix will be nrows M x nrows M
0736 //==============================================================================
0737 template <class T, unsigned int D1, unsigned int D2, class R>
0738 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) {
0739   SMatrix<T,D1,D2, MatRepStd<T,D1,D2> > tmp = lhs * rhs;
0740   typedef  SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym;
0741   SMatrixSym mret;
0742   AssignSym::Evaluate(mret,  tmp * Transpose(lhs)  );
0743   return mret;
0744 }
0745
0746 //==============================================================================
0747 // product: SMatrix/SMatrix calculate M * A * M^T where A is a symmetric matrix
0748 // return matrix will be nrowsM x nrows M
0749 // M is a matrix expression
0750 //==============================================================================
0751 template <class A, class T, unsigned int D1, unsigned int D2, class R>
0752 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) {
0753   SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = lhs * rhs;
0754   typedef  SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym;
0755   SMatrixSym mret;
0756   AssignSym::Evaluate(mret,  tmp * Transpose(lhs)  );
0757   return mret;
0758 }
0759
0760 #ifdef XXX
0761     // not needed (
0762 //==============================================================================
0763 // product: SMatrix/SMatrix calculate M * A * M where A and M are symmetric matrices
0764 // return matrix will be nrows M x nrows M
0765 //==============================================================================
0766 template <class T, unsigned int D1>
0767 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) {
0768   SMatrix<T,D1,D1, MatRepStd<T,D1,D1> > tmp = lhs * rhs;
0769   typedef  SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym;
0770   SMatrixSym mret;
0771   AssignSym::Evaluate(mret,  tmp * lhs  );
0772   return mret;
0773 }
0774 #endif
0775
0776
0777 /**
0778    Transpose Similarity Matrix Product :  B = U^T * A * U for A symmetric
0779    returning a symmetric matrix expression: \f$ B(i,j) = \sum_{k,l} U(k,i) * A(k,l) * U(l,j) \f$
0780
0781    @ingroup MatrixFunctions
0782  */
0783 //==============================================================================
0784 // product: SMatrix/SMatrix calculate M^T * A * M where A is a symmetric matrix
0785 // return matrix will be ncolsM x ncols M
0786 //==============================================================================
0787 template <class T, unsigned int D1, unsigned int D2, class R>
0788 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) {
0789   SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = rhs * lhs;
0790   typedef  SMatrix<T,D2,D2,MatRepSym<T,D2> > SMatrixSym;
0791   SMatrixSym mret;
0792   AssignSym::Evaluate(mret,  Transpose(lhs) * tmp );
0793   return mret;
0794 }
0795
0796 //==============================================================================
0797 // product: SMatrix/SMatrix calculate M^T * A * M where A is a symmetric matrix
0798 // return matrix will be ncolsM x ncols M
0799 // M is a matrix expression
0800 //==============================================================================
0801 template <class A, class T, unsigned int D1, unsigned int D2, class R>
0802 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) {
0803   SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = rhs * lhs;
0804   typedef  SMatrix<T,D2,D2,MatRepSym<T,D2> > SMatrixSym;
0805   SMatrixSym mret;
0806   AssignSym::Evaluate(mret,  Transpose(lhs) * tmp );
0807   return mret;
0808 }
0809
0810
0811
0812
0813
0814 // //==============================================================================
0815 // // Mult * (Expr * Expr, binary) with a symmetric result
0816 // // the operation is done only for half
0817 // //==============================================================================
0818 // template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
0819 // 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>
0820 //  operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
0821 //   typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, T,D> MatMulOp;
0822 //   return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
0823 // }
0824
0825
0826
0827 //==============================================================================
0828 // TensorMulOp
0829 //==============================================================================
0830 /**
0831    Class for Tensor Multiplication (outer product) of two vectors
0832    giving a matrix
0833
0834    @ingroup Expression
0835  */
0836 template <class Vector1, class Vector2>
0837 class TensorMulOp {
0838 public:
0839   ///
0840   TensorMulOp( const Vector1 & lhs, const Vector2 & rhs) :
0841     lhs_(lhs),
0842     rhs_(rhs) {}
0843
0844   ///
0845   ~TensorMulOp() {}
0846
0847   /// Vector2::kSize is the number of columns in the resulting matrix
0848   inline typename Vector1::value_type apply(unsigned int i) const {
0849     return lhs_.apply( i/ Vector2::kSize) * rhs_.apply( i % Vector2::kSize );
0850   }
0851   inline typename Vector1::value_type operator() (unsigned int i, unsigned j) const {
0852     return lhs_.apply(i) * rhs_.apply(j);
0853   }
0854
0855   inline bool IsInUse (const typename Vector1::value_type * ) const {
0856     return false;
0857   }
0858
0859
0860 protected:
0861
0862   const Vector1 & lhs_;
0863   const Vector2 & rhs_;
0864
0865 };
0866
0867
0868
0869 /**
0870    Tensor Vector Product : M(i,j) = v(i) * v(j)
0871    returning a matrix expression
0872
0873    @ingroup VectFunction
0874  */
0875
0876 #ifndef _WIN32
0877
0878     // Tensor Prod (use default MatRepStd for the returned expression
0879     // cannot make a symmetric matrix
0880 //==============================================================================
0881 // TensorProd (SVector x SVector)
0882 //==============================================================================
0883 template <class T, unsigned int D1, unsigned int D2>
0884 inline Expr<TensorMulOp<SVector<T,D1>, SVector<T,D2>  >, T, D1, D2 >
0885  TensorProd(const SVector<T,D1>& lhs, const SVector<T,D2>& rhs) {
0886   typedef TensorMulOp<SVector<T,D1>, SVector<T,D2> > TVMulOp;
0887   return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
0888 }
0889
0890 //==============================================================================
0891 // TensorProd (VecExpr x SVector)
0892 //==============================================================================
0893  template <class T, unsigned int D1, unsigned int D2, class A>
0894 inline Expr<TensorMulOp<VecExpr<A,T,D1>, SVector<T,D2>  >, T, D1, D2 >
0895  TensorProd(const VecExpr<A,T,D1>& lhs, const SVector<T,D2>& rhs) {
0896   typedef TensorMulOp<VecExpr<A,T,D1>, SVector<T,D2> > TVMulOp;
0897   return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
0898 }
0899
0900 //==============================================================================
0901 // TensorProd (SVector x VecExpr)
0902 //==============================================================================
0903  template <class T, unsigned int D1, unsigned int D2, class A>
0904 inline Expr<TensorMulOp<SVector<T,D1>, VecExpr<A,T,D2>  >, T, D1, D2 >
0905  TensorProd(const SVector<T,D1>& lhs, const VecExpr<A,T,D2>& rhs) {
0906   typedef TensorMulOp<SVector<T,D1>, VecExpr<A,T,D2> > TVMulOp;
0907   return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
0908 }
0909
0910
0911 //==============================================================================
0912 // TensorProd (VecExpr x VecExpr)
0913 //==============================================================================
0914  template <class T, unsigned int D1, unsigned int D2, class A, class B>
0915 inline Expr<TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2>  >, T, D1, D2 >
0916  TensorProd(const VecExpr<A,T,D1>& lhs, const VecExpr<B,T,D2>& rhs) {
0917   typedef TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2> > TVMulOp;
0918   return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
0919 }
0920
0921 #endif
0922 #ifdef _WIN32
0923 /// case of WINDOWS - problem using Expression (  C1001: INTERNAL COMPILER ERROR )
0924
0925 //==============================================================================
0926 // TensorProd (SVector x SVector)
0927 //==============================================================================
0928 template <class T, unsigned int D1, unsigned int D2>
0929 inline SMatrix<T,D1,D2,MatRepStd<T, D1, D2>> TensorProd(const SVector<T,D1>& lhs, const SVector<T,D2>& rhs) {
0930   SMatrix<T,D1,D2,MatRepStd<T, D1, D2>> tmp;
0931   for (unsigned int i=0; i< D1; ++i)
0932     for (unsigned int j=0; j< D2; ++j) {
0933       tmp(i,j) = lhs[i]*rhs[j];
0934     }
0935
0936   return tmp;
0937 }
0938 //==============================================================================
0939 // TensorProd (VecExpr x SVector)
0940 //==============================================================================
0941 template <class T, unsigned int D1, unsigned int D2, class A>
0942 inline SMatrix<T,D1,D2,MatRepStd<T, D1, D2>>  TensorProd(const VecExpr<A,T,D1>& lhs, const SVector<T,D2>& rhs) {
0943   SMatrix<T,D1,D2,MatRepStd<T, D1, D2>> tmp;
0944   for (unsigned int i=0; i< D1; ++i)
0945     for (unsigned int j=0; j< D2; ++j)
0946       tmp(i,j) = lhs.apply(i) * rhs.apply(j);
0947
0948   return tmp;
0949 }
0950 //==============================================================================
0951 // TensorProd (SVector x VecExpr)
0952 //==============================================================================
0953 template <class T, unsigned int D1, unsigned int D2, class A>
0954 inline SMatrix<T,D1,D2, MatRepStd<T, D1, D2>> TensorProd(const SVector<T,D1>& lhs, const VecExpr<A,T,D2>& rhs) {
0955   SMatrix<T,D1,D2,MatRepStd<T, D1, D2>> tmp;
0956   for (unsigned int i=0; i< D1; ++i)
0957     for (unsigned int j=0; j< D2; ++j)
0958       tmp(i,j) = lhs.apply(i) * rhs.apply(j);
0959
0960   return tmp;
0961 }
0962
0963 //==============================================================================
0964 // TensorProd (VecExpr x VecExpr)
0965 //==============================================================================
0966
0967 template <class T, unsigned int D1, unsigned int D2, class A, class B>
0968 inline SMatrix<T,D1,D2,MatRepStd<T, D1, D2>> TensorProd(const VecExpr<A,T,D1>& lhs, const VecExpr<B,T,D2>& rhs) {
0969   SMatrix<T,D1,D2,MatRepStd<T, D1, D2>> tmp;
0970   for (unsigned int i=0; i< D1; ++i)
0971     for (unsigned int j=0; j< D2; ++j)
0972       tmp(i,j) = lhs.apply(i) * rhs.apply(j);
0973
0974   return tmp;
0975 }
0976
0977
0978 #endif
0979
0980 // solving a positive defined symmetric linear system using Choleski decompositions
0981 // matrix will be decomposed and the returned vector will be overwritten in vec
0982 // If the user wants to pass const objects need to copy the matrices
0983 // It will work only for symmetric matrices
0984 template <class T, unsigned int D>
0985 bool SolveChol( SMatrix<T, D, D, MatRepSym<T, D>  > & mat,  SVector<T, D> & vec ) {
0986    CholeskyDecomp<T, D> decomp(mat);
0987    return decomp.Solve(vec);
0988 }
0989
0990 /// same function as before but not overwriting the matrix and returning a copy of the vector
0991 /// (this is the slow version)
0992 template <class T, unsigned int D>
0993 SVector<T,D> SolveChol( const SMatrix<T, D, D, MatRepSym<T, D>  > & mat,  const SVector<T, D> & vec, int & ifail  ) {
0994    SMatrix<T, D, D, MatRepSym<T, D> > atmp(mat);
0995    SVector<T,D> vret(vec);
0996    bool ok = SolveChol( atmp, vret);
0997    ifail =  (ok) ? 0 : -1;
0998    return vret;
0999 }
1000
1001
1002
1003   }  // namespace Math
1004
1005 }  // namespace ROOT
1006
1007
1008 #endif  /* ROOT_Math_MatrixFunctions */