Logo ROOT   6.18/05
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
48namespace 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//==============================================================================
59template <class T, unsigned int D1, unsigned int D2, class R>
60SVector<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//==============================================================================
81template <unsigned int I>
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//==============================================================================
94template <>
95struct 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//==============================================================================
106template <class Matrix, class Vector, unsigned int D2>
108public:
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
131protected:
132 const Matrix& lhs_;
133 const Vector& rhs_;
134};
135
136
137//==============================================================================
138// meta_col_dot
139//==============================================================================
140template <unsigned int I>
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//==============================================================================
154template <>
155struct 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 */
171template <class Vector, class Matrix, unsigned int D1>
173public:
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
195protected:
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//==============================================================================
209template <class T, unsigned int D1, unsigned int D2, class R>
213 return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
214}
215
216//==============================================================================
217// operator*: SMatrix * Expr<A,T,D2>
218//==============================================================================
219template <class A, class T, unsigned int D1, unsigned int D2, class R>
220inline VecExpr<VectorMatrixRowOp<SMatrix<T,D1,D2,R>, VecExpr<A,T,D2>, D2>, T, D1>
223 return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
224}
225
226//==============================================================================
227// operator*: Expr<A,T,D1,D2> * SVector
228//==============================================================================
229template <class A, class T, unsigned int D1, unsigned int D2, class R>
230inline 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) {
233 return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
234}
235
236//==============================================================================
237// operator*: Expr<A,T,D1,D2> * VecExpr<B,T,D2>
238//==============================================================================
239template <class A, class B, class T, unsigned int D1, unsigned int D2, class R>
240inline VecExpr<VectorMatrixRowOp<Expr<A,T,D1,D2,R>, VecExpr<B,T,D2>, D2>, T, D1>
243 return VecExpr<VMOp, T, D1>(VMOp(lhs,rhs));
244}
245
246//==============================================================================
247// operator*: SVector * SMatrix
248//==============================================================================
249template <class T, unsigned int D1, unsigned int D2, class R>
250inline VecExpr<VectorMatrixColOp<SVector<T,D1>, SMatrix<T,D1,D2,R>, D1>, T, D2>
253 return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
254}
255
256//==============================================================================
257// operator*: SVector * Expr<A,T,D1,D2>
258//==============================================================================
259template <class A, class T, unsigned int D1, unsigned int D2, class R>
260inline 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) {
263 return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
264}
265
266//==============================================================================
267// operator*: VecExpr<A,T,D1> * SMatrix
268//==============================================================================
269template <class A, class T, unsigned int D1, unsigned int D2, class R>
270inline VecExpr<VectorMatrixColOp<VecExpr<A,T,D1>, SMatrix<T,D1,D2,R>, D1>, T, D2>
273 return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
274}
275
276//==============================================================================
277// operator*: VecExpr<A,T,D1> * Expr<B,T,D1,D2>
278//==============================================================================
279template <class A, class B, class T, unsigned int D1, unsigned int D2, class R>
280inline VecExpr<VectorMatrixColOp<VecExpr<A,T,D1>, Expr<B,T,D1,D2,R>, D1>, T, D2>
283 return VecExpr<VMOp, T, D2>(VMOp(lhs,rhs));
284}
285
286//==============================================================================
287// meta_matrix_dot
288//==============================================================================
289template <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//==============================================================================
316template <>
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 */
345template <class MatrixA, class MatrixB, class T, unsigned int D>
347public:
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 {
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
369protected:
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//==============================================================================
384template < class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
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//==============================================================================
395template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
396inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
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//==============================================================================
406template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
407inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>,T,D>, T, D1, D2, typename MultPolicy<T,R1,R2>::RepType>
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//==============================================================================
417template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
418inline 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>
422}
423
424
425
426#ifdef XXX
427//==============================================================================
428// MatrixMulOp
429//==============================================================================
430template <class MatrixA, class MatrixB, unsigned int D>
431class MatrixMulOp {
432public:
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 {
443 }
444
445protected:
446 const MatrixA& lhs_;
447 const MatrixB& rhs_;
448};
449
450
451//==============================================================================
452// operator* (SMatrix * SMatrix, binary)
453//==============================================================================
454template < class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
455inline 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;
458 return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
459}
460
461//==============================================================================
462// operator* (SMatrix * Expr, binary)
463//==============================================================================
464template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
465inline 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;
468 return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
469}
470
471//==============================================================================
472// operator* (Expr * SMatrix, binary)
473//==============================================================================
474template <class A, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
475inline 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;
478 return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
479}
480
481//=============================================================================
482// operator* (Expr * Expr, binary)
483//=============================================================================
484template <class A, class B, class T, unsigned int D1, unsigned int D, unsigned int D2, class R1, class R2>
485inline 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;
488 return Expr<MatMulOp,T,D1,D2,typename MultPolicy<T,R1,R2>::RepType>(MatMulOp(lhs,rhs));
489}
490#endif
491
492//==============================================================================
493// TransposeOp
494//==============================================================================
495/**
496 Class for Transpose Operations
497
498 @ingroup Expression
499 */
500template <class Matrix, class T, unsigned int D1, unsigned int D2=D1>
502public:
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
522protected:
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//==============================================================================
536template <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//==============================================================================
547template <class A, class T, unsigned int D1, unsigned int D2, class R>
548inline Expr<TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2>, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
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//==============================================================================
562template <class T, unsigned int D1, unsigned int D2, class R>
563inline SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
564 Transpose(const SMatrix<T,D1,D2, R>& rhs) {
565 typedef TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2> MatTrOp;
566
567 return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
568 ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
569}
570
571//==============================================================================
572// transpose
573//==============================================================================
574template <class A, class T, unsigned int D1, unsigned int D2, class R>
575inline SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
576 Transpose(const Expr<A,T,D1,D2,R>& rhs) {
577 typedef TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2> MatTrOp;
578
579 return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
580 ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
581}
582
583#endif
584
585
586#ifdef OLD
587//==============================================================================
588// product: SMatrix/SVector calculate v^T * A * v
589//==============================================================================
590template <class T, unsigned int D, class R>
591inline 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//==============================================================================
598template <class T, unsigned int D, class R>
599inline 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//==============================================================================
606template <class A, class T, unsigned int D, class R>
607inline 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//==============================================================================
614template <class A, class T, unsigned int D, class R>
615inline 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//==============================================================================
622template <class A, class T, unsigned int D, class R>
623inline 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//==============================================================================
630template <class A, class T, unsigned int D, class R>
631inline 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//==============================================================================
638template <class A, class B, class T, unsigned int D, class R>
639inline 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//==============================================================================
646template <class A, class B, class T, unsigned int D, class R>
647inline 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//==============================================================================
662template <class T, unsigned int D, class R>
663inline 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//==============================================================================
670template <class T, unsigned int D, class R>
671inline 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//==============================================================================
678template <class A, class T, unsigned int D, class R>
679inline 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//==============================================================================
686template <class A, class T, unsigned int D, class R>
687inline 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//==============================================================================
694template <class A, class T, unsigned int D, class R>
695inline 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//==============================================================================
702template <class A, class T, unsigned int D, class R>
703inline 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//==============================================================================
710template <class A, class B, class T, unsigned int D, class R>
711inline 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//==============================================================================
718template <class A, class B, class T, unsigned int D, class R>
719inline 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//==============================================================================
735template <class T, unsigned int D1, unsigned int D2, class R>
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//==============================================================================
749template <class A, class T, unsigned int D1, unsigned int D2, class R>
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//==============================================================================
764template <class T, unsigned int D1>
765inline 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//==============================================================================
785template <class T, unsigned int D1, unsigned int D2, class R>
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//==============================================================================
799template <class A, class T, unsigned int D1, unsigned int D2, class R>
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 */
834template <class Vector1, class Vector2>
836public:
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
858protected:
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//==============================================================================
881template <class T, unsigned int D1, unsigned int 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>
892inline 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) {
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>
902inline 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) {
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>
913inline Expr<TensorMulOp<VecExpr<A,T,D1>, VecExpr<B,T,D2> >, T, D1, D2 >
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//==============================================================================
926template <class T, unsigned int D1, unsigned int D2>
927inline SMatrix<T,D1,D2,MatRepStd<T, D1, D2>> TensorProd(const SVector<T,D1>& lhs, const SVector<T,D2>& rhs) {
928 SMatrix<T,D1,D2,MatRepStd<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//==============================================================================
939template <class T, unsigned int D1, unsigned int D2, class A>
940inline SMatrix<T,D1,D2,MatRepStd<T, D1, D2>> TensorProd(const VecExpr<A,T,D1>& lhs, const SVector<T,D2>& rhs) {
941 SMatrix<T,D1,D2,MatRepStd<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//==============================================================================
951template <class T, unsigned int D1, unsigned int D2, class A>
952inline SMatrix<T,D1,D2, MatRepStd<T, D1, D2>> TensorProd(const SVector<T,D1>& lhs, const VecExpr<A,T,D2>& rhs) {
953 SMatrix<T,D1,D2,MatRepStd<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
965template <class T, unsigned int D1, unsigned int D2, class A, class B>
966inline SMatrix<T,D1,D2,MatRepStd<T, D1, D2>> TensorProd(const VecExpr<A,T,D1>& lhs, const VecExpr<B,T,D2>& rhs) {
967 SMatrix<T,D1,D2,MatRepStd<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
982template <class T, unsigned int D>
983bool 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)
990template <class T, unsigned int D>
991SVector<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 */
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:102
@ 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
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:164
#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.
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
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)