From 462b57ff12004a835fbbfcfcef07505c1a38ac61 Mon Sep 17 00:00:00 2001
From: Lorenzo Moneta <Lorenzo.Moneta@cern.ch>
Date: Tue, 25 Apr 2006 13:54:01 +0000
Subject: [PATCH] add IsInUse function in SMAtrix and Expression to create
automatically a temporary object in operations like A = B * A
git-svn-id: http://root.cern.ch/svn/root/trunk@14834 27541ba8-7e3a-0410-8455-c3a389f83636
---
smatrix/inc/Math/Expression.h | 47 +++++++++---
smatrix/inc/Math/HelperOps.h | 75 +++++++++++++++----
smatrix/inc/Math/MatrixFunctions.h | 43 ++++++++++-
.../inc/Math/MatrixRepresentationsStatic.h | 5 +-
smatrix/inc/Math/SMatrix.h | 33 ++++++--
smatrix/inc/Math/SMatrix.icc | 10 ++-
6 files changed, 178 insertions(+), 35 deletions(-)
diff --git a/smatrix/inc/Math/Expression.h b/smatrix/inc/Math/Expression.h
index 17d391aa8c8..6c4fa176749 100644
--- a/smatrix/inc/Math/Expression.h
+++ b/smatrix/inc/Math/Expression.h
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: Expression.h,v 1.8 2006/03/09 09:24:04 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: Expression.h,v 1.9 2006/03/20 17:11:44 moneta Exp $
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_Expression
@@ -53,11 +53,11 @@ namespace ROOT {
+// template <class T, unsigned int D, unsigned int D2> class MatRepStd;
- // template <class T, unsigned int D, unsigned int D2> class MatRepStd;
+template <class ExprType, class T, unsigned int D >
+class VecExpr {
- template <class ExprType, class T, unsigned int D >
- class VecExpr {
public:
typedef T value_type;
@@ -108,11 +108,11 @@ private:
- template <class T, unsigned int D, unsigned int D2> class MatRepStd;
+template <class T, unsigned int D, unsigned int D2> class MatRepStd;
- template <class ExprType, class T, unsigned int D, unsigned int D2 = 1,
- class R1=MatRepStd<T,D,D2> >
- class Expr {
+template <class ExprType, class T, unsigned int D, unsigned int D2 = 1,
+ class R1=MatRepStd<T,D,D2> >
+class Expr {
public:
typedef T value_type;
@@ -130,6 +130,15 @@ public:
inline T operator() (unsigned int i, unsigned j) const {
return rhs_(i,j);
}
+
+ /**
+ function to determine if any use operand
+ is being used (has same memory adress)
+ */
+ inline bool IsInUse (const T * p) const {
+ return rhs_.IsInUse(p);
+ }
+
#ifdef OLD_IMPL
@@ -214,6 +223,10 @@ public:
return Operator::apply(lhs_(i,j), rhs_(i,j) );
}
+ inline bool IsInUse (const T * p) const {
+ return lhs_.IsInUse(p) || rhs_.IsInUse(p);
+ }
+
protected:
const LHS& lhs_;
@@ -249,6 +262,11 @@ public:
return Operator::apply(lhs_(i,j), rhs_(i,j) );
}
+ inline bool IsInUse (const T * p) const {
+ // no need to check left since we copy it
+ return rhs_.IsInUse(p);
+ }
+
protected:
const LHS lhs_;
@@ -260,7 +278,7 @@ protected:
//==============================================================================
/**
Special case of BinaryOp where for the wight argument a copy is stored instead of a reference
- This is use in the coase for example of constant where we cannot store by reference
+ This is use in the case for example of constant where we cannot store by reference
but need to copy since Constant is a temporary object
*/
//==============================================================================
@@ -282,6 +300,11 @@ public:
return Operator::apply(lhs_(i,j), rhs_(i,j) );
}
+ inline bool IsInUse (const T * p) const {
+ // no need for right since we copied
+ return lhs_.IsInUse(p);
+ }
+
protected:
const LHS& lhs_;
@@ -318,6 +341,10 @@ public:
return Operator::apply(rhs_(i,j));
}
+ inline bool IsInUse (const T * p) const {
+ return rhs_.IsInUse(p);
+ }
+
protected:
const RHS& rhs_;
@@ -349,6 +376,8 @@ public:
inline T operator() (unsigned int /*i */, unsigned int /*j */ ) const { return rhs_; }
+ //inline bool IsInUse (const T * ) const { return false; }
+
protected:
const T rhs_; // no need for reference. It is a fundamental type normally
diff --git a/smatrix/inc/Math/HelperOps.h b/smatrix/inc/Math/HelperOps.h
index 0f57015f931..51e3c1d0410 100644
--- a/smatrix/inc/Math/HelperOps.h
+++ b/smatrix/inc/Math/HelperOps.h
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: HelperOps.h,v 1.6 2006/03/30 10:33:44 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: HelperOps.h,v 1.7 2006/03/30 16:18:05 moneta Exp $
// Authors: J. Palacios 2006
#ifndef ROOT_Math_HelperOps_h
@@ -39,14 +39,29 @@ namespace ROOT {
{
static void Evaluate(SMatrix<T,D1,D2,R1>& lhs, const Expr<A,T,D1,D2,R2>& rhs)
{
- //for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = rhs.apply(i);
- unsigned int l = 0;
- for(unsigned int i=0; i<D1; ++i)
- for(unsigned int j=0; j<D2; ++j) {
- lhs.fRep[l] = rhs(i,j);
- l++;
- }
+ if (! rhs.IsInUse(lhs.begin() ) ) {
+ unsigned int l = 0;
+ for(unsigned int i=0; i<D1; ++i)
+ for(unsigned int j=0; j<D2; ++j) {
+ lhs.fRep[l] = rhs(i,j);
+ l++;
+ }
+ }
+ // lhs is in use in expression, need to create a temporary with the result
+ else {
+ T tmp[D1*D2];
+ unsigned int l = 0;
+ for(unsigned int i=0; i<D1; ++i)
+ for(unsigned int j=0; j<D2; ++j) {
+ tmp[l] = rhs(i,j);
+ l++;
+ }
+ // copy now the temp object
+ for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = tmp[i];
+ }
+
}
+
};
// specialization in case of symmetric expression to symmetric matrices :
template <class T,
@@ -58,13 +73,28 @@ namespace ROOT {
static void Evaluate(SMatrix<T,D1,D2,MatRepSym<T,D1> >& lhs,
const Expr<A,T,D1,D2,MatRepSym<T,D1> >& rhs)
{
- unsigned int l = 0;
- for(unsigned int i=0; i<D1; ++i)
- // storage of symmetric matrix is in lower block
- for(unsigned int j=0; j<=i; ++j) {
- lhs.fRep.Array()[l] = rhs(i,j);
- l++;
- }
+ if (! rhs.IsInUse(lhs.begin() ) ) {
+ unsigned int l = 0;
+ for(unsigned int i=0; i<D1; ++i)
+ // storage of symmetric matrix is in lower block
+ for(unsigned int j=0; j<=i; ++j) {
+ lhs.fRep.Array()[l] = rhs(i,j);
+ l++;
+ }
+ }
+ // create a temporary object to store result
+ else {
+ std::cout << "create temp obj since it is in use" << std::endl;
+ T tmp[MatRepSym<T,D1>::kSize];
+ unsigned int l = 0;
+ for(unsigned int i=0; i<D1; ++i)
+ for(unsigned int j=0; j<=i; ++j) {
+ tmp[l] = rhs(i,j);
+ l++;
+ }
+ // copy now the temp object
+ for(unsigned int i=0; i<MatRepSym<T,D1>::kSize; ++i) lhs.fRep.Array()[i] = tmp[i];
+ }
}
};
@@ -101,6 +131,21 @@ namespace ROOT {
l++;
}
}
+ // assign a symmetric matrix from a general matrix that we assume is symmetric
+ template <class T,
+ unsigned int D,
+ class R>
+ static void Evaluate(SMatrix<T,D,D,MatRepSym<T,D> >& lhs, const SMatrix<T,D,D,R>& rhs)
+ {
+ //for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = rhs.apply(i);
+ unsigned int l = 0;
+ for(unsigned int i=0; i<D; ++i)
+ // storage of symmetric matrix is in lower block
+ for(unsigned int j=0; j<=i; ++j) {
+ lhs.fRep.Array()[l] = rhs(i,j);
+ l++;
+ }
+ }
}; // struct AssignSym
diff --git a/smatrix/inc/Math/MatrixFunctions.h b/smatrix/inc/Math/MatrixFunctions.h
index 60f48d149e6..a190e2d4123 100644
--- a/smatrix/inc/Math/MatrixFunctions.h
+++ b/smatrix/inc/Math/MatrixFunctions.h
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: MatrixFunctions.h,v 1.8 2006/03/09 09:24:04 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: MatrixFunctions.h,v 1.9 2006/03/20 17:11:44 moneta Exp $
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_MatrixFunctions
@@ -308,6 +308,11 @@ public:
return meta_matrix_dot<D-1>::g(lhs_, rhs_, i, j);
}
+ inline bool IsInUse (const T * p) const {
+ return lhs_.IsInUse(p) || rhs_.IsInUse(p);
+ }
+
+
protected:
const MatrixA& lhs_;
const MatrixB& rhs_;
@@ -358,6 +363,7 @@ inline Expr<MatrixMulOp<Expr<A,T,D1,D,R1>, Expr<B,T,D,D2,R2>,T,D>, T, D1, D2, ty
}
+
#ifdef XXX
//==============================================================================
// MatrixMulOp
@@ -445,6 +451,10 @@ public:
return rhs_( j, i);
}
+ inline bool IsInUse (const T * p) const {
+ return rhs_.IsInUse(p);
+ }
+
protected:
const Matrix& rhs_;
};
@@ -472,6 +482,37 @@ inline Expr<TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2>, T, D2, D1, typename TranspPo
return Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs));
}
+
+#ifdef ENABLE_TEMPORARIES_TRANSPOSE
+// sometimes is faster to create a temp, not clear why
+
+//==============================================================================
+// transpose
+//==============================================================================
+template <class T, unsigned int D1, unsigned int D2, class R>
+inline SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
+ Transpose(const SMatrix<T,D1,D2, R>& rhs) {
+ typedef TransposeOp<SMatrix<T,D1,D2,R>,T,D1,D2> MatTrOp;
+
+ return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
+ ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
+}
+
+//==============================================================================
+// transpose
+//==============================================================================
+template <class A, class T, unsigned int D1, unsigned int D2, class R>
+inline SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
+ Transpose(const Expr<A,T,D1,D2,R>& rhs) {
+ typedef TransposeOp<Expr<A,T,D1,D2,R>,T,D1,D2> MatTrOp;
+
+ return SMatrix< T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>
+ ( Expr<MatTrOp, T, D2, D1, typename TranspPolicy<T,D1,D2,R>::RepType>(MatTrOp(rhs)) );
+}
+
+#endif
+
+
#ifdef OLD
//==============================================================================
// product: SMatrix/SVector calculate v^T * A * v
diff --git a/smatrix/inc/Math/MatrixRepresentationsStatic.h b/smatrix/inc/Math/MatrixRepresentationsStatic.h
index 6e40a4764cf..acdf79fe414 100644
--- a/smatrix/inc/Math/MatrixRepresentationsStatic.h
+++ b/smatrix/inc/Math/MatrixRepresentationsStatic.h
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: MatrixRepresentationsStatic.h,v 1.4 2006/03/14 17:11:19 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: MatrixRepresentationsStatic.h,v 1.5 2006/03/17 15:11:35 moneta Exp $
// Authors: L. Moneta, J. Palacios 2006
#ifndef ROOT_Math_MatrixRepresentationsStatic_h
@@ -198,7 +198,8 @@ namespace ROOT {
private:
T fArray[kSize];
- RowOffsets<D> * fOff;
+
+ RowOffsets<D> * fOff; //! transient
};
diff --git a/smatrix/inc/Math/SMatrix.h b/smatrix/inc/Math/SMatrix.h
index 2f9c214e0ca..f28e456da0f 100644
--- a/smatrix/inc/Math/SMatrix.h
+++ b/smatrix/inc/Math/SMatrix.h
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: SMatrix.h,v 1.18 2006/03/30 16:18:05 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: SMatrix.h,v 1.19 2006/04/20 13:13:21 moneta Exp $
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_SMatrix
@@ -119,7 +119,7 @@ public:
/**
construct from an expression.
- In case of symmetric matrices doesnot work if expression is of type general
+ In case of symmetric matrices does not work if expression is of type general
matrices. In case one needs to force the assignment from general to symmetric, one can use the
ROOT::Math::AssignSym::Evaluate function.
*/
@@ -220,7 +220,7 @@ public:
};
#endif
/** @name --- Access functions --- */
- /** access the parse tree. Index starts from zero and follows C convention for accessing
+ /** access the parse tree with the index starting from zero and following the C convention for the order in accessing
the matrix elements.
*/
T apply(unsigned int i) const;
@@ -230,7 +230,19 @@ public:
/// return pointer to internal array
T* Array();
- /** @name --- STL-like interface --- */
+ /** @name --- STL-like interface ---
+ The iterators access the matrix element in the order how they are
+ stored in memory. The C (row-major) convention is used, and in the
+ case of symmetric matrices the iterator spans only the lower diagonal
+ block. For example for a symmetric 3x3 matrices the order of the 6
+ elements \f${a_0,...a_5}\f$ is:
+ \f[
+ M = \left( \begin{array}{ccc}
+ a_0 & a_1 & a_3 \\
+ a1 & a_2 & a_4 \\
+ a3 & a_4 & a_5 \end{array} \right)
+ \f]
+ */
/** STL iterator interface. */
iterator begin();
@@ -280,11 +292,11 @@ public:
bool operator<(const Expr<A,T,D1,D2,R2>& rhs) const;
/**
- read only access to matrix element. Indeces start from 0
+ read only access to matrix element, with indices starting from 0
*/
const T& operator()(unsigned int i, unsigned int j) const;
/**
- read/write access to matrix element. Indeces start from 0
+ read/write access to matrix element with indices starting from 0
*/
T& operator()(unsigned int i, unsigned int j);
@@ -483,7 +495,14 @@ SMatrix<T,D1,D2,R>& Place_in_col(const VecExpr<A,T,D>& rhs,
#endif
-
+ /**
+ check if matrix is sharing same memory location.
+ This functionis used by the expression to avoid the alias problem when
+ evaluating them. In case matrix is in use, a temporary object is automatically
+ created evaluating the expression. Then the correct result is obtained for operations
+ like A = B * A
+ */
+ bool IsInUse(const T* p) const;
// submatrices
diff --git a/smatrix/inc/Math/SMatrix.icc b/smatrix/inc/Math/SMatrix.icc
index 2347c346f79..51298651941 100644
--- a/smatrix/inc/Math/SMatrix.icc
+++ b/smatrix/inc/Math/SMatrix.icc
@@ -1,4 +1,4 @@
-// @(#)root/smatrix:$Name: $:$Id: SMatrix.icc,v 1.18 2006/03/30 10:33:44 moneta Exp $
+// @(#)root/smatrix:$Name: $:$Id: SMatrix.icc,v 1.19 2006/03/30 16:18:05 moneta Exp $
// Authors: T. Glebe, L. Moneta 2005
#ifndef ROOT_Math_SMatrix_icc
@@ -874,6 +874,14 @@ SMatrix<T,D1,D2,R>::SMatrix(const SVector<T, N > & v, bool lower ) {
}
+template <class T, unsigned int D1, unsigned int D2, class R>
+bool SMatrix<T,D1,D2,R>::IsInUse( const T * p) const {
+ return p == fRep.Array();
+}
+
+
+
+
} // namespace Math
} // namespace ROOT
--
GitLab