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