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