ROOT  6.06/09
Reference Guide
MatrixFunctions.h
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Authors: T. Glebe, L. Moneta 2005
3 
4 #ifndef ROOT_Math_MatrixFunctions
5 #define ROOT_Math_MatrixFunctions
6 // ********************************************************************
7 //
8 // source:
9 //
10 // type: source code
11 //
12 // created: 20. Mar 2001
13 //
14 // author: Thorsten Glebe
15 // HERA-B Collaboration
16 // Max-Planck-Institut fuer Kernphysik
17 // Saupfercheckweg 1
18 // 69117 Heidelberg
19 // Germany
20 // E-mail: T.Glebe@mpi-hd.mpg.de
21 //
22 // Description: Functions/Operators special to Matrix
23 //
24 // changes:
25 // 20 Mar 2001 (TG) creation
26 // 20 Mar 2001 (TG) Added Matrix * Vector multiplication
27 // 21 Mar 2001 (TG) Added transpose, product
28 // 11 Apr 2001 (TG) transpose() speed improvment by removing rows(), cols()
29 // through static members of Matrix and Expr class
30 //
31 // ********************************************************************
32 
33 //doxygen tag
34 /**
35  @defgroup MatrixFunctions Matrix Template Functions
36  @ingroup SMatrixGroup
37 
38  These function apply to matrices (and also Matrix expression) and can return a
39  matrix expression of a particular defined type, like in the matrix multiplication or
40  a vector, like in the matrix-vector product or a scalar like in the Similarity vector-matrix product.
41 */
42 
43 #ifndef ROOT_Math_BinaryOpPolicy
44 #include "Math/BinaryOpPolicy.h"
45 #endif
46 #ifndef ROOT_Math_Expression
47 #include "Math/Expression.h"
48 #endif
49 #ifndef ROOT_Math_HelperOps
50 #include "Math/HelperOps.h"
51 #endif
52 #ifndef ROOT_Math_CholeskyDecomp
53 #include "Math/CholeskyDecomp.h"
54 #endif
55 
56 namespace ROOT {
57 
58  namespace Math {
59 
60  template <class T, unsigned int D>
61  class SVector;
62 
63 #ifdef XXX
64 //==============================================================================
65 // SMatrix * SVector
66 //==============================================================================
67 template <class T, unsigned int D1, unsigned int D2, class R>
68 SVector<T,D1> operator*(const SMatrix<T,D1,D2,R>& rhs, const SVector<T,D2>& lhs)
69 {
70  SVector<T,D1> tmp;
71  for(unsigned int i=0; i<D1; ++i) {
72  const unsigned int rpos = i*D2;
73  for(unsigned int j=0; j<D2; ++j) {
74  tmp[i] += rhs.apply(rpos+j) * lhs.apply(j);
75  }
76  }
77  return tmp;
78 }
79 #endif
80 
81 
82 // matrix-vector product:
83 // use apply(i) function for matrices. Tested (11/05/06) with using (i,j) but
84 // performances are slightly worse (not clear why)
85 
86 //==============================================================================
87 // meta_row_dot
88 //==============================================================================
89 template <unsigned int I>
90 struct meta_row_dot {
91  template <class A, class B>
92  static inline typename A::value_type f(const A& lhs, const B& rhs,
93  const unsigned int offset) {
94  return lhs.apply(offset+I) * rhs.apply(I) + meta_row_dot<I-1>::f(lhs,rhs,offset);
95  }
96 };
97 
98 
99 //==============================================================================
100 // meta_row_dot<0>
101 //==============================================================================
102 template <>
103 struct meta_row_dot<0> {
104  template <class A, class B>
105  static inline typename A::value_type f(const A& lhs, const B& rhs,
106  const unsigned int offset) {
107  return lhs.apply(offset) * rhs.apply(0);
108  }
109 };
110 
111 //==============================================================================
112 // VectorMatrixRowOp
113 //==============================================================================
114 template <class Matrix, class Vector, unsigned int D2>
116 public:
117 
118  typedef typename Vector::value_type T;
119 
120  ///
121  VectorMatrixRowOp(const Matrix& lhs, const Vector& rhs) :
122  lhs_(lhs), rhs_(rhs) {}
123 
124  ///
126 
127  /// calc \f$ \sum_{j} a_{ij} * v_j \f$
128  inline typename Matrix::value_type apply(unsigned int i) const {
129  return meta_row_dot<D2-1>::f(lhs_, rhs_, i*D2);
130  }
131 
132  // check if passed pointer is in use
133  // check only the vector since this is a vector expression
134  inline bool IsInUse (const T * p) const {
135  return rhs_.IsInUse(p);
136  }
137 
138 
139 protected:
140  const Matrix& lhs_;
141  const Vector& rhs_;
142 };
143 
144 
145 //==============================================================================
146 // meta_col_dot
147 //==============================================================================
148 template <unsigned int I>
149 struct meta_col_dot {
150  template <class Matrix, class Vector>
151  static inline typename Matrix::value_type f(const Matrix& lhs, const Vector& rhs,
152  const unsigned int offset) {
153  return lhs.apply(Matrix::kCols*I+offset) * rhs.apply(I) +
154  meta_col_dot<I-1>::f(lhs,rhs,offset);
155  }
156 };
157 
158 
159 //==============================================================================
160 // meta_col_dot<0>
161 //==============================================================================
162 template <>
163 struct meta_col_dot<0> {
164  template <class Matrix, class Vector>
165  static inline typename Matrix::value_type f(const Matrix& lhs, const Vector& rhs,
166  const unsigned int offset) {
167  return lhs.apply(offset) * rhs.apply(0);
168  }
169 };
170 
171 //==============================================================================
172 // VectorMatrixColOp
173 //==============================================================================
174 /**
175  Class for Vector-Matrix multiplication
176 
177  @ingroup Expression
178  */
179 template <class Vector, class Matrix, unsigned int D1>
181 public:
182 
183  typedef typename Vector::value_type T;
184  ///
185  VectorMatrixColOp(const Vector& lhs, const Matrix& rhs) :
186  lhs_(lhs), rhs_(rhs) {}
187 
188  ///
190 
191  /// calc \f$ \sum_{j} a_{ij} * v_j \f$
192  inline typename Matrix::value_type apply(unsigned int i) const {
193  return meta_col_dot<D1-1>::f(rhs_, lhs_, i);
194  }
195 
196  // check if passed pointer is in use
197  // check only the vector since this is a vector expression
198  inline bool IsInUse (const T * p) const {
199  return lhs_.IsInUse(p);
200  }
201 
202 
203 protected:
204  const Vector& lhs_;
205  const Matrix& rhs_;
206 };
207 
208 /**
209  Matrix * Vector multiplication \f$ a(i) = \sum_{j} M(i,j) * b(j) \f$
210  returning a vector expression
211 
212  @ingroup MatrixFunctions
213  */
214 //==============================================================================
215 // operator*: SMatrix * SVector
216 //==============================================================================
217 template <class T, unsigned int D1, unsigned int D2, class R>
218 inline VecExpr<VectorMatrixRowOp<SMatrix<T,D1,D2,R>,SVector<T,D2>, D2>, T, D1>
219  operator*(const SMatrix<T,D1,D2,R>& lhs, const SVector<T,D2>& rhs) {
220  typedef VectorMatrixRowOp<SMatrix<T,D1,D2,R>,SVector<T,D2>, D2> VMOp;
221  return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
222 }
223 
224 //==============================================================================
225 // operator*: SMatrix * Expr<A,T,D2>
226 //==============================================================================
227 template <class A, class T, unsigned int D1, unsigned int D2, class R>
228 inline VecExpr<VectorMatrixRowOp<SMatrix<T,D1,D2,R>, VecExpr<A,T,D2>, D2>, T, D1>
229  operator*(const SMatrix<T,D1,D2,R>& lhs, const VecExpr<A,T,D2>& rhs) {
230  typedef VectorMatrixRowOp<SMatrix<T,D1,D2,R>,VecExpr<A,T,D2>, D2> VMOp;
231  return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
232 }
233 
234 //==============================================================================
235 // operator*: Expr<A,T,D1,D2> * SVector
236 //==============================================================================
237 template <class A, class T, unsigned int D1, unsigned int D2, class R>
238 inline VecExpr<VectorMatrixRowOp<Expr<A,T,D1,D2,R>, SVector<T,D2>, D2>, T, D1>
239  operator*(const Expr<A,T,D1,D2,R>& lhs, const SVector<T,D2>& rhs) {
240  typedef VectorMatrixRowOp<Expr<A,T,D1,D2,R>,SVector<T,D2>, D2> VMOp;
241  return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
242 }
243 
244 //==============================================================================
245 // operator*: Expr<A,T,D1,D2> * VecExpr<B,T,D2>
246 //==============================================================================
247 template <class A, class B, class T, unsigned int D1, unsigned int D2, class R>
248 inline VecExpr<VectorMatrixRowOp<Expr<A,T,D1,D2,R>, VecExpr<B,T,D2>, D2>, T, D1>
249  operator*(const Expr<A,T,D1,D2,R>& lhs, const VecExpr<B,T,D2>& rhs) {
250  typedef VectorMatrixRowOp<Expr<A,T,D1,D2,R>,VecExpr<B,T,D2>, D2> VMOp;
251  return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
252 }
253 
254 //==============================================================================
255 // operator*: SVector * SMatrix
256 //==============================================================================
257 template <class T, unsigned int D1, unsigned int D2, class R>
258 inline VecExpr<VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2,R>, D1>, T, D2>
259  operator*(const SVector<T,D1>& lhs, const SMatrix<T,D1,D2,R>& rhs) {
260  typedef VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2,R>, D1> VMOp;
261  return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
262 }
263 
264 //==============================================================================
265 // operator*: SVector * Expr<A,T,D1,D2>
266 //==============================================================================
267 template <class A, class T, unsigned int D1, unsigned int D2, class R>
268 inline VecExpr<VectorMatrixColOp<SVector<T,D1>, Expr<A,T,D1,D2,R>, D1>, T, D2>
269  operator*(const SVector<T,D1>& lhs, const Expr<A,T,D1,D2,R>& rhs) {
270  typedef VectorMatrixColOp<SVector<T,D1>, Expr<A,T,D1,D2,R>, D1> VMOp;
271  return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
272 }
273 
274 //==============================================================================
275 // operator*: VecExpr<A,T,D1> * SMatrix
276 //==============================================================================
277 template <class A, class T, unsigned int D1, unsigned int D2, class R>
278 inline VecExpr<VectorMatrixColOp<VecExpr<A,T,D1>, SMatrix<T,D1,D2,R>, D1>, T, D2>
279  operator*(const VecExpr<A,T,D1>& lhs, const SMatrix<T,D1,D2,R>& rhs) {
280  typedef VectorMatrixColOp<VecExpr<A,T,D1>, SMatrix<T,D1,D2,R>, D1> VMOp;
281  return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
282 }
283 
284 //==============================================================================
285 // operator*: VecExpr<A,T,D1> * Expr<B,T,D1,D2>
286 //==============================================================================
287 template <class A, class B, class T, unsigned int D1, unsigned int D2, class R>
288 inline VecExpr<VectorMatrixColOp<VecExpr<A,T,D1>, Expr<B,T,D1,D2,R>, D1>, T, D2>
289  operator*(const VecExpr<A,T,D1>& lhs, const Expr<B,T,D1,D2,R>& rhs) {
290  typedef VectorMatrixColOp<VecExpr<A,T,D1>, Expr<B,T,D1,D2,R>, D1> VMOp;
291  return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
292 }
293 
294 //==============================================================================
295 // meta_matrix_dot
296 //==============================================================================
297 template <unsigned int I>
299 
300  template <class MatrixA, class MatrixB>
301  static inline typename MatrixA::value_type f(const MatrixA& lhs,
302  const MatrixB& rhs,
303  const unsigned int offset) {
304  return lhs.apply(offset/MatrixB::kCols*MatrixA::kCols + I) *
305  rhs.apply(MatrixB::kCols*I + offset%MatrixB::kCols) +
306  meta_matrix_dot<I-1>::f(lhs,rhs,offset);
307  }
308 
309  // multiplication using i and j indeces
310  template <class MatrixA, class MatrixB>
311  static inline typename MatrixA::value_type g(const MatrixA& lhs,
312  const MatrixB& rhs,
313  unsigned int i,
314  unsigned int j) {
315  return lhs(i, I) * rhs(I , j) +
316  meta_matrix_dot<I-1>::g(lhs,rhs,i,j);
317  }
318 };
319 
320 
321 //==============================================================================
322 // meta_matrix_dot<0>
323 //==============================================================================
324 template <>
325 struct meta_matrix_dot<0> {
326 
327  template <class MatrixA, class MatrixB>
328  static inline typename MatrixA::value_type f(const MatrixA& lhs,
329  const MatrixB& rhs,
330  const unsigned int offset) {
331  return lhs.apply(offset/MatrixB::kCols*MatrixA::kCols) *
332  rhs.apply(offset%MatrixB::kCols);
333  }
334 
335  // multiplication using i and j
336  template <class MatrixA, class MatrixB>
337  static inline typename MatrixA::value_type g(const MatrixA& lhs,
338  const MatrixB& rhs,
339  unsigned int i, unsigned int j) {
340  return lhs(i,0) * rhs(0,j);
341  }
342 
343 };
344 
345 //==============================================================================
346 // MatrixMulOp
347 //==============================================================================
348 /**
349  Class for Matrix-Matrix multiplication
350 
351  @ingroup Expression
352  */
353 template <class MatrixA, class MatrixB, class T, unsigned int D>
354 class MatrixMulOp {
355 public:
356  ///
357  MatrixMulOp(const MatrixA& lhs, const MatrixB& rhs) :
358  lhs_(lhs), rhs_(rhs) {}
359 
360  ///
362 
363  /// calc \f$\sum_{j} a_{ik} * b_{kj}\f$
364  inline T apply(unsigned int i) const {
365  return meta_matrix_dot<D-1>::f(lhs_, rhs_, i);
366  }
367 
368  inline T operator() (unsigned int i, unsigned j) const {
369  return meta_matrix_dot<D-1>::g(lhs_, rhs_, i, j);
370  }
371 
372  inline bool IsInUse (const T * p) const {
373  return lhs_.IsInUse(p) || rhs_.IsInUse(p);
374  }
375 
376 
377 protected:
378  const MatrixA& lhs_;
379  const MatrixB& rhs_;
380 };
381 
382 
383 /**
384  Matrix * Matrix multiplication , \f$ C(i,j) = \sum_{k} A(i,k) * B(k,j)\f$
385  returning a matrix expression
386 
387  @ingroup MatrixFunctions
388  */
389 //==============================================================================
390 // operator* (SMatrix * SMatrix, binary)
391 //==============================================================================
392 template < class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
394  operator*(const SMatrix<T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
395  typedef MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, T,D> MatMulOp;
396  return Expr<MatMulOp,T,D1,D2,
397  typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
398 }
399 
400 //==============================================================================
401 // operator* (SMatrix * Expr, binary)
402 //==============================================================================
403 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
404 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>
405  operator*(const SMatrix<T,D1,D,R1>& lhs, const Expr<A,T,D,D2,R2>& rhs) {
406  typedef MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>,T,D> MatMulOp;
407  return Expr<MatMulOp,T,D1,D2,
408  typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
409 }
410 
411 //==============================================================================
412 // operator* (Expr * SMatrix, binary)
413 //==============================================================================
414 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
415 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>
416  operator*(const Expr<A,T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
417  typedef MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>,T,D> MatMulOp;
418  return Expr<MatMulOp,T,D1,D2,
419  typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
420 }
421 
422 //==============================================================================
423 // operator* (Expr * Expr, binary)
424 //==============================================================================
425 template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
426 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>
427  operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
428  typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, T,D> MatMulOp;
430 }
431 
432 
433 
434 #ifdef XXX
435 //==============================================================================
436 // MatrixMulOp
437 //==============================================================================
438 template <class MatrixA, class MatrixB, unsigned int D>
439 class MatrixMulOp {
440 public:
441  ///
442  MatrixMulOp(const MatrixA& lhs, const MatrixB& rhs) :
443  lhs_(lhs), rhs_(rhs) {}
444 
445  ///
446  ~MatrixMulOp() {}
447 
448  /// calc $\sum_{j} a_{ik} * b_{kj}$
449  inline typename MatrixA::value_type apply(unsigned int i) const {
450  return meta_matrix_dot<D-1>::f(lhs_, rhs_, i);
451  }
452 
453 protected:
454  const MatrixA& lhs_;
455  const MatrixB& rhs_;
456 };
457 
458 
459 //==============================================================================
460 // operator* (SMatrix * SMatrix, binary)
461 //==============================================================================
462 template < class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
463 inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
464  operator*(const SMatrix<T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
465  typedef MatrixMulOp<SMatrix<T,D1,D,R1>, SMatrix<T,D,D2,R2>, D> MatMulOp;
466  return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
467 }
468 
469 //==============================================================================
470 // operator* (SMatrix * Expr, binary)
471 //==============================================================================
472 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
473 inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
474  operator*(const SMatrix<T,D1,D,R1>& lhs, const Expr<A,T,D,D2,R2>& rhs) {
475  typedef MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>, D> MatMulOp;
476  return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
477 }
478 
479 //==============================================================================
480 // operator* (Expr * SMatrix, binary)
481 //==============================================================================
482 template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
483 inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>, D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
484  operator*(const Expr<A,T,D1,D,R1>& lhs, const SMatrix<T,D,D2,R2>& rhs) {
485  typedef MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>, D> MatMulOp;
486  return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
487 }
488 
489 //=============================================================================
490 // operator* (Expr * Expr, binary)
491 //=============================================================================
492 template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
493 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>
494  operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
495  typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, D> MatMulOp;
496  return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
497 }
498 #endif
499 
500 //==============================================================================
501 // TransposeOp
502 //==============================================================================
503 /**
504  Class for Transpose Operations
505 
506  @ingroup Expression
507  */
508 template <class Matrix, class T, unsigned int D1, unsigned int D2=D1>
509 class TransposeOp {
510 public:
511  ///
512  TransposeOp( const Matrix& rhs) :
513  rhs_(rhs) {}
514 
515  ///
517 
518  ///
519  inline T apply(unsigned int i) const {
520  return rhs_.apply( (i%D1)*D2 + i/D1);
521  }
522  inline T operator() (unsigned int i, unsigned j) const {
523  return rhs_( j, i);
524  }
525 
526  inline bool IsInUse (const T * p) const {
527  return rhs_.IsInUse(p);
528  }
529 
530 protected:
531  const Matrix& rhs_;
532 };
533 
534 
535 /**
536  Matrix Transpose B(i,j) = A(j,i)
537  returning a matrix expression
538 
539  @ingroup MatrixFunctions
540  */
541 //==============================================================================
542 // transpose
543 //==============================================================================
544 template <class T, unsigned int D1, unsigned int D2, class R>
547  typedef TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2> MatTrOp;
548 
550 }
551 
552 //==============================================================================
553 // transpose
554 //==============================================================================
555 template <class A, class T, unsigned int D1, unsigned int D2, class R>
556 inline Expr<TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2>, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
557  Transpose(const Expr<A,T,D1,D2,R>& rhs) {
558  typedef TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2> MatTrOp;
559 
561 }
562 
563 
564 #ifdef ENABLE_TEMPORARIES_TRANSPOSE
565 // sometimes is faster to create a temp, not clear why
566 
567 //==============================================================================
568 // transpose
569 //==============================================================================
570 template <class T, unsigned int D1, unsigned int D2, class R>
571 inline SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
572  Transpose(const SMatrix<T,D1,D2, R>& rhs) {
573  typedef TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2> MatTrOp;
574 
575  return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
576  ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
577 }
578 
579 //==============================================================================
580 // transpose
581 //==============================================================================
582 template <class A, class T, unsigned int D1, unsigned int D2, class R>
583 inline SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
584  Transpose(const Expr<A,T,D1,D2,R>& rhs) {
585  typedef TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2> MatTrOp;
586 
587  return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
588  ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
589 }
590 
591 #endif
592 
593 
594 #ifdef OLD
595 //==============================================================================
596 // product: SMatrix/SVector calculate v^T * A * v
597 //==============================================================================
598 template <class T, unsigned int D, class R>
599 inline T Product(const SMatrix<T,D,D,R>& lhs, const SVector<T,D>& rhs) {
600  return Dot(rhs, lhs * rhs);
601 }
602 
603 //==============================================================================
604 // product: SVector/SMatrix calculate v^T * A * v
605 //==============================================================================
606 template <class T, unsigned int D, class R>
607 inline T Product(const SVector<T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
608  return Dot(lhs, rhs * lhs);
609 }
610 
611 //==============================================================================
612 // product: SMatrix/Expr calculate v^T * A * v
613 //==============================================================================
614 template <class A, class T, unsigned int D, class R>
615 inline T Product(const SMatrix<T,D,D,R>& lhs, const VecExpr<A,T,D>& rhs) {
616  return Dot(rhs, lhs * rhs);
617 }
618 
619 //==============================================================================
620 // product: Expr/SMatrix calculate v^T * A * v
621 //==============================================================================
622 template <class A, class T, unsigned int D, class R>
623 inline T Product(const VecExpr<A,T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
624  return Dot(lhs, rhs * lhs);
625 }
626 
627 //==============================================================================
628 // product: SVector/Expr calculate v^T * A * v
629 //==============================================================================
630 template <class A, class T, unsigned int D, class R>
631 inline T Product(const SVector<T,D>& lhs, const Expr<A,T,D,D,R>& rhs) {
632  return Dot(lhs, rhs * lhs);
633 }
634 
635 //==============================================================================
636 // product: Expr/SVector calculate v^T * A * v
637 //==============================================================================
638 template <class A, class T, unsigned int D, class R>
639 inline T Product(const Expr<A,T,D,D,R>& lhs, const SVector<T,D>& rhs) {
640  return Dot(rhs, lhs * rhs);
641 }
642 
643 //==============================================================================
644 // product: Expr/Expr calculate v^T * A * v
645 //==============================================================================
646 template <class A, class B, class T, unsigned int D, class R>
647 inline T Product(const Expr<A,T,D,D,R>& lhs, const VecExpr<B,T,D>& rhs) {
648  return Dot(rhs, lhs * rhs);
649 }
650 
651 //==============================================================================
652 // product: Expr/Expr calculate v^T * A * v
653 //==============================================================================
654 template <class A, class B, class T, unsigned int D, class R>
655 inline T Product(const VecExpr<A,T,D>& lhs, const Expr<B,T,D,D,R>& rhs) {
656  return Dot(lhs, rhs * lhs);
657 }
658 #endif
659 
660 /**
661  Similarity Vector - Matrix Product: v^T * A * v
662  returning a scalar value of type T \f$ s = \sum_{i,j} v(i) * A(i,j) * v(j)\f$
663 
664  @ingroup MatrixFunctions
665  */
666 
667 //==============================================================================
668 // product: SMatrix/SVector calculate v^T * A * v
669 //==============================================================================
670 template <class T, unsigned int D, class R>
671 inline T Similarity(const SMatrix<T,D,D,R>& lhs, const SVector<T,D>& rhs) {
672  return Dot(rhs, lhs * rhs);
673 }
674 
675 //==============================================================================
676 // product: SVector/SMatrix calculate v^T * A * v
677 //==============================================================================
678 template <class T, unsigned int D, class R>
679 inline T Similarity(const SVector<T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
680  return Dot(lhs, rhs * lhs);
681 }
682 
683 //==============================================================================
684 // product: SMatrix/Expr calculate v^T * A * v
685 //==============================================================================
686 template <class A, class T, unsigned int D, class R>
687 inline T Similarity(const SMatrix<T,D,D,R>& lhs, const VecExpr<A,T,D>& rhs) {
688  return Dot(rhs, lhs * rhs);
689 }
690 
691 //==============================================================================
692 // product: Expr/SMatrix calculate v^T * A * v
693 //==============================================================================
694 template <class A, class T, unsigned int D, class R>
695 inline T Similarity(const VecExpr<A,T,D>& lhs, const SMatrix<T,D,D,R>& rhs) {
696  return Dot(lhs, rhs * lhs);
697 }
698 
699 //==============================================================================
700 // product: SVector/Expr calculate v^T * A * v
701 //==============================================================================
702 template <class A, class T, unsigned int D, class R>
703 inline T Similarity(const SVector<T,D>& lhs, const Expr<A,T,D,D,R>& rhs) {
704  return Dot(lhs, rhs * lhs);
705 }
706 
707 //==============================================================================
708 // product: Expr/SVector calculate v^T * A * v
709 //==============================================================================
710 template <class A, class T, unsigned int D, class R>
711 inline T Similarity(const Expr<A,T,D,D,R>& lhs, const SVector<T,D>& rhs) {
712  return Dot(rhs, lhs * rhs);
713 }
714 
715 //==============================================================================
716 // product: Expr/Expr calculate v^T * A * v
717 //==============================================================================
718 template <class A, class B, class T, unsigned int D, class R>
719 inline T Similarity(const Expr<A,T,D,D,R>& lhs, const VecExpr<B,T,D>& rhs) {
720  return Dot(rhs, lhs * rhs);
721 }
722 
723 //==============================================================================
724 // product: Expr/Expr calculate v^T * A * v
725 //==============================================================================
726 template <class A, class B, class T, unsigned int D, class R>
727 inline T Similarity(const VecExpr<A,T,D>& lhs, const Expr<B,T,D,D,R>& rhs) {
728  return Dot(lhs, rhs * lhs);
729 }
730 
731 
732 /**
733  Similarity Matrix Product : B = U * A * U^T for A symmetric
734  returning a symmetric matrix expression:
735  \f$ B(i,j) = \sum_{k,l} U(i,k) * A(k,l) * U(j,l) \f$
736 
737  @ingroup MatrixFunctions
738  */
739 //==============================================================================
740 // product: SMatrix/SMatrix calculate M * A * M^T where A is a symmetric matrix
741 // return matrix will be nrows M x nrows M
742 //==============================================================================
743 template <class T, unsigned int D1, unsigned int D2, class R>
745  SMatrix<T,D1,D2, MatRepStd<T,D1,D2> > tmp = lhs * rhs;
746  typedef SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym;
747  SMatrixSym mret;
748  AssignSym::Evaluate(mret, tmp * Transpose(lhs) );
749  return mret;
750 }
751 
752 //==============================================================================
753 // product: SMatrix/SMatrix calculate M * A * M^T where A is a symmetric matrix
754 // return matrix will be nrowsM x nrows M
755 // M is a matrix expression
756 //==============================================================================
757 template <class A, class T, unsigned int D1, unsigned int D2, class R>
759  SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = lhs * rhs;
760  typedef SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym;
761  SMatrixSym mret;
762  AssignSym::Evaluate(mret, tmp * Transpose(lhs) );
763  return mret;
764 }
765 
766 #ifdef XXX
767  // not needed (
768 //==============================================================================
769 // product: SMatrix/SMatrix calculate M * A * M where A and M are symmetric matrices
770 // return matrix will be nrows M x nrows M
771 //==============================================================================
772 template <class T, unsigned int D1>
773 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) {
774  SMatrix<T,D1,D1, MatRepStd<T,D1,D1> > tmp = lhs * rhs;
775  typedef SMatrix<T,D1,D1,MatRepSym<T,D1> > SMatrixSym;
776  SMatrixSym mret;
777  AssignSym::Evaluate(mret, tmp * lhs );
778  return mret;
779 }
780 #endif
781 
782 
783 /**
784  Transpose Similarity Matrix Product : B = U^T * A * U for A symmetric
785  returning a symmetric matrix expression: \f$ B(i,j) = \sum_{k,l} U(k,i) * A(k,l) * U(l,j) \f$
786 
787  @ingroup MatrixFunctions
788  */
789 //==============================================================================
790 // product: SMatrix/SMatrix calculate M^T * A * M where A is a symmetric matrix
791 // return matrix will be ncolsM x ncols M
792 //==============================================================================
793 template <class T, unsigned int D1, unsigned int D2, class R>
795  SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = rhs * lhs;
796  typedef SMatrix<T,D2,D2,MatRepSym<T,D2> > SMatrixSym;
797  SMatrixSym mret;
798  AssignSym::Evaluate(mret, Transpose(lhs) * tmp );
799  return mret;
800 }
801 
802 //==============================================================================
803 // product: SMatrix/SMatrix calculate M^T * A * M where A is a symmetric matrix
804 // return matrix will be ncolsM x ncols M
805 // M is a matrix expression
806 //==============================================================================
807 template <class A, class T, unsigned int D1, unsigned int D2, class R>
809  SMatrix<T,D1,D2,MatRepStd<T,D1,D2> > tmp = rhs * lhs;
810  typedef SMatrix<T,D2,D2,MatRepSym<T,D2> > SMatrixSym;
811  SMatrixSym mret;
812  AssignSym::Evaluate(mret, Transpose(lhs) * tmp );
813  return mret;
814 }
815 
816 
817 
818 
819 
820 // //==============================================================================
821 // // Mult * (Expr * Expr, binary) with a symmetric result
822 // // the operation is done only for half
823 // //==============================================================================
824 // template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
825 // 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>
826 // operator*(const Expr<A,T,D1,D,R1>& lhs, const Expr<B,T,D,D2,R2>& rhs) {
827 // typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, T,D> MatMulOp;
828 // return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
829 // }
830 
831 
832 
833 //==============================================================================
834 // TensorMulOp
835 //==============================================================================
836 /**
837  Class for Tensor Multiplication (outer product) of two vectors
838  giving a matrix
839 
840  @ingroup Expression
841  */
842 template <class Vector1, class Vector2>
843 class TensorMulOp {
844 public:
845  ///
846  TensorMulOp( const Vector1 & lhs, const Vector2 & rhs) :
847  lhs_(lhs),
848  rhs_(rhs) {}
849 
850  ///
852 
853  /// Vector2::kSize is the number of columns in the resulting matrix
854  inline typename Vector1::value_type apply(unsigned int i) const {
855  return lhs_.apply( i/ Vector2::kSize) * rhs_.apply( i % Vector2::kSize );
856  }
857  inline typename Vector1::value_type operator() (unsigned int i, unsigned j) const {
858  return lhs_.apply(i) * rhs_.apply(j);
859  }
860 
861  inline bool IsInUse (const typename Vector1::value_type * ) const {
862  return false;
863  }
864 
865 
866 protected:
867 
868  const Vector1 & lhs_;
869  const Vector2 & rhs_;
870 
871 };
872 
873 
874 
875 /**
876  Tensor Vector Product : M(i,j) = v(i) * v(j)
877  returning a matrix expression
878 
879  @ingroup VectFunction
880  */
881 
882 #ifndef _WIN32
883 
884  // Tensor Prod (use default MatRepStd for the returned expression
885  // cannot make a symmetric matrix
886 //==============================================================================
887 // TensorProd (SVector x SVector)
888 //==============================================================================
889 template <class T, unsigned int D1, unsigned int D2>
890 inline Expr<TensorMulOp<SVector<T,D1>, SVector<T,D2> >, T, D1, D2 >
891  TensorProd(const SVector<T,D1>& lhs, const SVector<T,D2>& rhs) {
892  typedef TensorMulOp<SVector<T,D1>, SVector<T,D2> > TVMulOp;
893  return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
894 }
895 
896 //==============================================================================
897 // TensorProd (VecExpr x SVector)
898 //==============================================================================
899  template <class T, unsigned int D1, unsigned int D2, class A>
900 inline Expr<TensorMulOp<VecExpr<A,T,D1>, SVector<T,D2> >, T, D1, D2 >
901  TensorProd(const VecExpr<A,T,D1>& lhs, const SVector<T,D2>& rhs) {
902  typedef TensorMulOp<VecExpr<A,T,D1>, SVector<T,D2> > TVMulOp;
903  return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
904 }
905 
906 //==============================================================================
907 // TensorProd (SVector x VecExpr)
908 //==============================================================================
909  template <class T, unsigned int D1, unsigned int D2, class A>
910 inline Expr<TensorMulOp<SVector<T,D1>, VecExpr<A,T,D2> >, T, D1, D2 >
911  TensorProd(const SVector<T,D1>& lhs, const VecExpr<A,T,D2>& rhs) {
912  typedef TensorMulOp<SVector<T,D1>, VecExpr<A,T,D2> > TVMulOp;
913  return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
914 }
915 
916 
917 //==============================================================================
918 // TensorProd (VecExpr x VecExpr)
919 //==============================================================================
920  template <class T, unsigned int D1, unsigned int D2, class A, class B>
921 inline Expr<TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2> >, T, D1, D2 >
922  TensorProd(const VecExpr<A,T,D1>& lhs, const VecExpr<B,T,D2>& rhs) {
923  typedef TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2> > TVMulOp;
924  return Expr<TVMulOp,T,D1,D2>(TVMulOp(lhs,rhs));
925 }
926 
927 #endif
928 #ifdef _WIN32
929 /// case of WINDOWS - problem using Expression ( C1001: INTERNAL COMPILER ERROR )
930 
931 //==============================================================================
932 // TensorProd (SVector x SVector)
933 //==============================================================================
934 template <class T, unsigned int D1, unsigned int D2>
935 inline SMatrix<T,D1,D2> TensorProd(const SVector<T,D1>& lhs, const SVector<T,D2>& rhs) {
936  SMatrix<T,D1,D2> tmp;
937  for (unsigned int i=0; i< D1; ++i)
938  for (unsigned int j=0; j< D2; ++j) {
939  tmp(i,j) = lhs[i]*rhs[j];
940  }
941 
942  return tmp;
943 }
944 //==============================================================================
945 // TensorProd (VecExpr x SVector)
946 //==============================================================================
947  template <class T, unsigned int D1, unsigned int D2, class A>
948 inline SMatrix<T,D1,D2> TensorProd(const VecExpr<A,T,D1>& lhs, const SVector<T,D2>& rhs) {
949  SMatrix<T,D1,D2> tmp;
950  for (unsigned int i=0; i< D1; ++i)
951  for (unsigned int j=0; j< D2; ++j)
952  tmp(i,j) = lhs.apply(i) * rhs.apply(j);
953 
954  return tmp;
955 }
956 //==============================================================================
957 // TensorProd (SVector x VecExpr)
958 //==============================================================================
959  template <class T, unsigned int D1, unsigned int D2, class A>
960 inline SMatrix<T,D1,D2> TensorProd(const SVector<T,D1>& lhs, const VecExpr<A,T,D2>& rhs) {
961  SMatrix<T,D1,D2> tmp;
962  for (unsigned int i=0; i< D1; ++i)
963  for (unsigned int j=0; j< D2; ++j)
964  tmp(i,j) = lhs.apply(i) * rhs.apply(j);
965 
966  return tmp;
967 }
968 
969 //==============================================================================
970 // TensorProd (VecExpr x VecExpr)
971 //==============================================================================
972 
973  template <class T, unsigned int D1, unsigned int D2, class A, class B>
974 inline SMatrix<T,D1,D2 > TensorProd(const VecExpr<A,T,D1>& lhs, const VecExpr<B,T,D2>& rhs) {
975  SMatrix<T,D1,D2> tmp;
976  for (unsigned int i=0; i< D1; ++i)
977  for (unsigned int j=0; j< D2; ++j)
978  tmp(i,j) = lhs.apply(i) * rhs.apply(j);
979 
980  return tmp;
981 }
982 
983 
984 #endif
985 
986 // solving a positive defined symmetric linear system using Choleski decompositions
987 // matrix will be decomposed and the returned vector will be overwritten in vec
988 // If the user wants to pass const objects need to copy the matrices
989 // It will work only for symmetric matrices
990 template <class T, unsigned int D>
991 bool SolveChol( SMatrix<T, D, D, MatRepSym<T, D> > & mat, SVector<T, D> & vec ) {
992  CholeskyDecomp<T, D> decomp(mat);
993  return decomp.Solve(vec);
994 }
995 
996 /// same function as before but not overwriting the matrix and returning a copy of the vector
997 /// (this is the slow version)
998 template <class T, unsigned int D>
999 SVector<T,D> SolveChol( const SMatrix<T, D, D, MatRepSym<T, D> > & mat, const SVector<T, D> & vec, int & ifail ) {
1001  SVector<T,D> vret(vec);
1002  bool ok = SolveChol( atmp, vret);
1003  ifail = (ok) ? 0 : -1;
1004  return vret;
1005 }
1006 
1007 
1008 
1009  } // namespace Math
1010 
1011 } // namespace ROOT
1012 
1013 
1014 #endif /* ROOT_Math_MatrixFunctions */
T Dot(const SVector< T, D > &lhs, const SVector< T, D > &rhs)
Vector dot product.
Definition: Functions.h:166
Class for Vector-Matrix multiplication.
T Similarity(const SMatrix< T, D, D, R > &lhs, const SVector< T, D > &rhs)
Similarity Vector - Matrix Product: v^T * A * v returning a scalar value of type T ...
static Double_t Product(const Double_t *x, const Float_t *y)
Product.
Definition: TCTUB.cxx:101
static double B[]
Class for Tensor Multiplication (outer product) of two vectors giving a matrix.
static MatrixA::value_type g(const MatrixA &lhs, const MatrixB &rhs, unsigned int i, unsigned int j)
Class for Matrix-Matrix multiplication.
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
VectorMatrixRowOp(const Matrix &lhs, const Vector &rhs)
double T(double x)
Definition: ChebyshevPol.h:34
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
Definition: HelperOps.h:35
bool IsInUse(const typename Vector1::value_type *) const
T apply(unsigned int i) const
Matrix::value_type apply(unsigned int i) const
calc
Class for Transpose Operations.
Expr< TransposeOp< SMatrix< T, D1, D2, R >, T, D1, D2 >, T, D2, D1, typename TranspPolicy< T, D1, D2, R >::RepType > Transpose(const SMatrix< T, D1, D2, R > &rhs)
Matrix Transpose B(i,j) = A(j,i) returning a matrix expression.
bool SolveChol(SMatrix< T, D, D, MatRepSym< T, D > > &mat, SVector< T, D > &vec)
T operator()(unsigned int i, unsigned j) const
bool IsInUse(const T *p) const
TensorMulOp(const Vector1 &lhs, const Vector2 &rhs)
static double A[]
header file containing the templated implementation of matrix inversion routines for use with ROOT's ...
bool IsInUse(const T *p) const
class to compute the Cholesky decomposition of a matrix
SMatrix: a generic fixed size D1 x D2 Matrix class.
TransposeOp(const Matrix &rhs)
Expression wrapper class for Matrix objects.
Definition: Expression.h:134
Vector1::value_type apply(unsigned int i) const
Vector2::kSize is the number of columns in the resulting matrix.
T operator()(unsigned int i, unsigned j) const
VectorMatrixColOp(const Vector &lhs, const Matrix &rhs)
bool IsInUse(const T *p) const
static Matrix::value_type f(const Matrix &lhs, const Vector &rhs, const unsigned int offset)
T apply(unsigned int i) const
calc
bool Solve(V &rhs) const
solves a linear system for the given right hand side
Expr< TensorMulOp< SVector< T, D1 >, SVector< T, D2 > >, T, D1, D2 > TensorProd(const SVector< T, D1 > &lhs, const SVector< T, D2 > &rhs)
Tensor Vector Product : M(i,j) = v(i) * v(j) returning a matrix expression.
static void Evaluate(SMatrix< T, D, D, MatRepSym< T, D > > &lhs, const Expr< A, T, D, D, R > &rhs)
assign a symmetric matrix from an expression
Definition: HelperOps.h:156
bool IsInUse(const T *p) const
MatrixMulOp(const MatrixA &lhs, const MatrixB &rhs)
static A::value_type f(const A &lhs, const B &rhs, const unsigned int offset)
static MatrixA::value_type f(const MatrixA &lhs, const MatrixB &rhs, const unsigned int offset)
static A::value_type f(const A &lhs, const B &rhs, const unsigned int offset)
Namespace for new Math classes and functions.
static MatrixA::value_type g(const MatrixA &lhs, const MatrixB &rhs, unsigned int i, unsigned int j)
Expression wrapper class for Vector objects.
Definition: Expression.h:64
Vector1::value_type operator()(unsigned int i, unsigned j) const
static MatrixA::value_type f(const MatrixA &lhs, const MatrixB &rhs, const unsigned int offset)
SMatrix< T, D2, D2, MatRepSym< T, D2 > > SimilarityT(const SMatrix< T, D1, D2, R > &lhs, const SMatrix< T, D1, D1, MatRepSym< T, D1 > > &rhs)
Transpose Similarity Matrix Product : B = U^T * A * U for A symmetric returning a symmetric matrix ex...
static Matrix::value_type f(const Matrix &lhs, const Vector &rhs, const unsigned int offset)
#define I(x, y, z)
Matrix::value_type apply(unsigned int i) const
calc
AxisAngle operator*(RotationX const &r1, AxisAngle const &r2)
Multiplication of an axial rotation by an AxisAngle.
SVector: a generic fixed size Vector class.