4 #ifndef ROOT_Math_MatrixFunctions 5 #define ROOT_Math_MatrixFunctions 52 template <
class T,
unsigned int D>
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)
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);
81 template <
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) {
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);
106 template <
class Matrix,
class Vector,
unsigned int D2>
110 typedef typename Vector::value_type
T;
114 lhs_(lhs), rhs_(rhs) {}
120 inline typename Matrix::value_type
apply(
unsigned int i)
const {
127 return rhs_.IsInUse(p);
140 template <
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) +
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);
171 template <
class Vector,
class Matrix,
unsigned int D1>
175 typedef typename Vector::value_type
T;
178 lhs_(lhs), rhs_(rhs) {}
184 inline typename Matrix::value_type
apply(
unsigned int i)
const {
191 return lhs_.IsInUse(p);
209 template <
class T,
unsigned int D1,
unsigned int D2,
class R>
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>
229 template <
class A,
class T,
unsigned int D1,
unsigned int D2,
class R>
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>
249 template <
class T,
unsigned int D1,
unsigned int D2,
class R>
259 template <
class A,
class T,
unsigned int D1,
unsigned int D2,
class R>
269 template <
class A,
class T,
unsigned int D1,
unsigned int D2,
class R>
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>
289 template <
unsigned int I>
292 template <
class MatrixA,
class MatrixB>
293 static inline typename MatrixA::value_type
f(
const MatrixA& lhs,
295 const unsigned int offset) {
296 return lhs.apply(offset/MatrixB::kCols*MatrixA::kCols +
I) *
297 rhs.apply(MatrixB::kCols*
I + offset%MatrixB::kCols) +
302 template <
class MatrixA,
class MatrixB>
303 static inline typename MatrixA::value_type
g(
const MatrixA& lhs,
307 return lhs(i,
I) * rhs(
I , j) +
319 template <
class MatrixA,
class MatrixB>
320 static inline typename MatrixA::value_type
f(
const MatrixA& lhs,
322 const unsigned int offset) {
323 return lhs.apply(offset/MatrixB::kCols*MatrixA::kCols) *
324 rhs.apply(offset%MatrixB::kCols);
328 template <
class MatrixA,
class MatrixB>
329 static inline typename MatrixA::value_type
g(
const MatrixA& lhs,
331 unsigned int i,
unsigned int j) {
332 return lhs(i,0) * rhs(0,j);
345 template <
class MatrixA,
class MatrixB,
class T,
unsigned int D>
350 lhs_(lhs), rhs_(rhs) {}
356 inline T
apply(
unsigned int i)
const {
365 return lhs_.IsInUse(p) || rhs_.IsInUse(p);
384 template <
class T,
unsigned int D1,
unsigned int D,
unsigned int D2,
class R1,
class R2>
385 inline Expr<MatrixMulOp<SMatrix<T,D1,D,R1>,
SMatrix<T,D,D2,R2>,
T,D>,
T, D1, D2,
typename MultPolicy<T,R1,R2>::RepType>
388 return Expr<MatMulOp,
T,D1,D2,
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>
398 typedef MatrixMulOp<SMatrix<T,D1,D,R1>, Expr<A,T,D,D2,R2>,
T,D> MatMulOp;
399 return Expr<MatMulOp,
T,D1,D2,
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>
409 typedef MatrixMulOp<Expr<A,T,D1,D,R1>, SMatrix<T,D,D2,R2>,
T,D> MatMulOp;
410 return Expr<MatMulOp,
T,D1,D2,
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>
420 typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>,
T,D> MatMulOp;
430 template <
class MatrixA,
class MatrixB,
unsigned int D>
434 MatrixMulOp(
const MatrixA& lhs,
const MatrixB& rhs) :
435 lhs_(lhs), rhs_(rhs) {}
441 inline typename MatrixA::value_type
apply(
unsigned int i)
const {
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>
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>
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>
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>
487 typedef MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>, D> MatMulOp;
500 template <
class Matrix,
class T,
unsigned int D1,
unsigned int D2=D1>
511 inline T
apply(
unsigned int i)
const {
512 return rhs_.apply( (i%D1)*D2 + i/D1);
519 return rhs_.IsInUse(p);
536 template <
class T,
unsigned int D1,
unsigned int D2,
class R>
537 inline Expr<TransposeOp<SMatrix<T,D1,D2,R>,
T,D1,D2>,
T, D2, D1,
typename TranspPolicy<T,D1,D2,R>::RepType>
547 template <
class A,
class T,
unsigned int D1,
unsigned int D2,
class R>
548 inline Expr<TransposeOp<Expr<A,T,D1,D2,R>,
T,D1,D2>,
T, D2, D1,
typename TranspPolicy<T,D1,D2,R>::RepType>
556 #ifdef ENABLE_TEMPORARIES_TRANSPOSE 562 template <
class T,
unsigned int D1,
unsigned int D2,
class R>
574 template <
class A,
class T,
unsigned int D1,
unsigned int D2,
class R>
576 Transpose(
const Expr<A,T,D1,D2,R>& rhs) {
590 template <
class T,
unsigned int D,
class R>
592 return Dot(rhs, lhs * rhs);
598 template <
class T,
unsigned int D,
class R>
600 return Dot(lhs, rhs * lhs);
606 template <
class A,
class T,
unsigned int D,
class R>
608 return Dot(rhs, lhs * rhs);
614 template <
class A,
class T,
unsigned int D,
class R>
616 return Dot(lhs, rhs * lhs);
622 template <
class A,
class T,
unsigned int D,
class R>
624 return Dot(lhs, rhs * lhs);
630 template <
class A,
class T,
unsigned int D,
class R>
632 return Dot(rhs, lhs * rhs);
638 template <
class A,
class B,
class T,
unsigned int D,
class R>
640 return Dot(rhs, lhs * rhs);
646 template <
class A,
class B,
class T,
unsigned int D,
class R>
648 return Dot(lhs, rhs * lhs);
662 template <
class T,
unsigned int D,
class R>
664 return Dot(rhs, lhs * rhs);
670 template <
class T,
unsigned int D,
class R>
672 return Dot(lhs, rhs * lhs);
678 template <
class A,
class T,
unsigned int D,
class R>
680 return Dot(rhs, lhs * rhs);
686 template <
class A,
class T,
unsigned int D,
class R>
688 return Dot(lhs, rhs * lhs);
694 template <
class A,
class T,
unsigned int D,
class R>
696 return Dot(lhs, rhs * lhs);
702 template <
class A,
class T,
unsigned int D,
class R>
704 return Dot(rhs, lhs * rhs);
710 template <
class A,
class B,
class T,
unsigned int D,
class R>
712 return Dot(rhs, lhs * rhs);
718 template <
class A,
class B,
class T,
unsigned int D,
class R>
720 return Dot(lhs, rhs * lhs);
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) {
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) {
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) {
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) {
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) {
834 template <
class Vector1,
class Vector2>
846 inline typename Vector1::value_type
apply(
unsigned int i)
const {
849 inline typename Vector1::value_type
operator() (
unsigned int i,
unsigned j)
const {
850 return lhs_.apply(i) * rhs_.apply(j);
853 inline bool IsInUse (
const typename Vector1::value_type * )
const {
881 template <
class T,
unsigned int D1,
unsigned int D2>
891 template <
class T,
unsigned int D1,
unsigned int D2,
class A>
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 >
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 >
926 template <
class T,
unsigned int D1,
unsigned int D2>
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];
939 template <
class T,
unsigned int D1,
unsigned int D2,
class A>
942 for (
unsigned int i=0; i< D1; ++i)
943 for (
unsigned int j=0; j< D2; ++j)
951 template <
class T,
unsigned int D1,
unsigned int D2,
class A>
954 for (
unsigned int i=0; i< D1; ++i)
955 for (
unsigned int j=0; j< D2; ++j)
965 template <
class T,
unsigned int D1,
unsigned int D2,
class A,
class B>
968 for (
unsigned int i=0; i< D1; ++i)
969 for (
unsigned int j=0; j< D2; ++j)
982 template <
class T,
unsigned int D>
985 return decomp.
Solve(vec);
990 template <
class T,
unsigned int D>
995 ifail = (ok) ? 0 : -1;
T Dot(const SVector< T, D > &lhs, const SVector< T, D > &rhs)
Vector dot product.
Class for Vector-Matrix multiplication.
T Similarity(const SMatrix< T, D, D, R > &lhs, const SVector< T, D > &rhs)
Similarity Vector - Matrix Product: v^T * A * v returning a scalar value of type T ...
static Double_t Product(const Double_t *x, const Float_t *y)
Product.
Class for Tensor Multiplication (outer product) of two vectors giving a matrix.
Class for Matrix-Matrix multiplication.
Namespace for new ROOT classes and functions.
Matrix::value_type apply(unsigned int i) const
calc
VectorMatrixRowOp(const Matrix &lhs, const Vector &rhs)
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
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)
TRObject operator()(const T1 &t1) const
header file containing the templated implementation of matrix inversion routines for use with ROOT'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.
VectorMatrixColOp(const Vector &lhs, const Matrix &rhs)
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
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
Matrix::value_type apply(unsigned int i) const
calc
MatrixMulOp(const MatrixA &lhs, const MatrixB &rhs)
Namespace for new Math classes and functions.
T apply(unsigned int i) const
calc
Expression wrapper class for Vector objects.
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
Vector1::value_type apply(unsigned int i) const
Vector2::kSize is the number of columns in the resulting matrix.
decltype(auto) constexpr apply(F &&f, Tuple &&t)
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