Skip to content

Extend templation in GenVector to better support vectorised Vc types. #394

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
vgvassilev merged 54 commits into from
Mar 15, 2017
Merged

Extend templation in GenVector to better support vectorised Vc types. #394

Show file tree
Hide file tree
Changes from 26 commits
Commits
Show all changes
54 commits
Select commit Hold shift + click to select a range
76a8f9d
Template the GenVector Plane3D, Transform3D and Translation3D classes…
cjones051073 Feb 25, 2017
7e9c2b7
Use the 'using namespace std' trick with Cartesian3D and Plane3D to f…
cjones051073 Feb 25, 2017
b15d62a
Extend the 'using namespace std' trick to a few more classes
cjones051073 Feb 25, 2017
2b2fc17
Change class to typename
cjones051073 Feb 25, 2017
2e74f0d
Explicitly cast some constants to the scalar type (needed for Vc types)
cjones051073 Feb 25, 2017
c3175c6
Remove empty implementation files
cjones051073 Feb 25, 2017
9bdfb8e
Template the Boost classes
cjones051073 Feb 25, 2017
c97245c
Remove empty implementation files for Boost classes
cjones051073 Feb 25, 2017
f12d368
Clean up some indentations
cjones051073 Feb 25, 2017
bbbbbd0
Add some missing template types to the Boost classes
cjones051073 Feb 25, 2017
9a0ee77
Remove empty Boost implementation file
cjones051073 Feb 25, 2017
371d50d
Small cleanups
cjones051073 Feb 25, 2017
dcb0b30
Clean up some initialisations
cjones051073 Feb 25, 2017
29e7e30
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Feb 27, 2017
14e73b4
Add Scalar/Vector Plane3D::Normalise() methods using SFINAE
cjones051073 Feb 27, 2017
27bbb57
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Feb 27, 2017
196a00b
Roll back some changes for now, as they are causing problems in the t…
cjones051073 Feb 27, 2017
0ad5828
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Feb 27, 2017
1928c75
Add specific vectorised implementations of a few methods
cjones051073 Feb 28, 2017
5ccdaea
Extend ostream operator for Displacement and Position vectors to supp…
cjones051073 Feb 28, 2017
ef1aea9
Fix the vector abs method (not fabs) and Vector Unit() for vectorised…
cjones051073 Feb 28, 2017
6cac16a
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Feb 28, 2017
68ba28a
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 1, 2017
939b7cb
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 4, 2017
d9b810b
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 6, 2017
776ce1a
Add override qualifier
cjones051073 Mar 6, 2017
fbe74d2
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 10, 2017
12b43b0
Fix initialisation bug
cjones051073 Mar 10, 2017
8367bdf
remove space
cjones051073 Mar 10, 2017
0f404ba
Some small changes following review comments
cjones051073 Mar 10, 2017
65a701d
remove spaces
cjones051073 Mar 10, 2017
4ba8a1b
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 10, 2017
ff94a5a
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 11, 2017
b8ca822
Fix formatting using clang-format
cjones051073 Mar 11, 2017
b32cb6a
some more 'using namespace std' updates
cjones051073 Mar 11, 2017
febe783
Fix more clang format issues
cjones051073 Mar 11, 2017
a628a2a
fix (hopefully) one more clang format issue
cjones051073 Mar 11, 2017
89b2d26
clang format is picky...
cjones051073 Mar 11, 2017
6bf0368
Merge branch 'master' of github.com:root-mirror/root into GenVector-A…
cjones051073 Mar 11, 2017
bd96bd0
revert changes in VectorUtils.h
cjones051073 Mar 13, 2017
dea8d05
Add new test for GenVector with Vc types
cjones051073 Mar 13, 2017
a18221a
extend test to 96 photons
cjones051073 Mar 13, 2017
eb0b8ed
reorder headers
cjones051073 Mar 13, 2017
7da1945
remove comment that number of test photons must be multiple of 8, as …
cjones051073 Mar 13, 2017
a271eab
Add <cmath> header includes
cjones051073 Mar 13, 2017
decbe11
revert back to using std:: as it seems to cause issues with gcc build…
cjones051073 Mar 14, 2017
58c6315
Include Vc before ROOT in test application. Needed for std:: support...
cjones051073 Mar 14, 2017
17f42a3
Add some more missing std::
cjones051073 Mar 14, 2017
434396d
More math updates
cjones051073 Mar 14, 2017
b145899
Change enable_if implementation in GenVector
cjones051073 Mar 14, 2017
ad8d3a5
update comment
cjones051073 Mar 14, 2017
b018a03
Extend the GenVector Vc test with a timing test that asserts that the…
cjones051073 Mar 14, 2017
bb1e386
add fabs on timing diff and extend Vc genvector expected diff to 25%
cjones051073 Mar 14, 2017
1786383
do not return a failure from the Vc speed test, as test could be run …
cjones051073 Mar 14, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion math/genvector/inc/Math/GenVector/BitReproducible.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ class BitReproducibleException : public std::exception
public:
BitReproducibleException(const std::string & w) throw() : fMsg(w) {}
~BitReproducibleException() throw() {}
const char* what() const throw() { return fMsg.c_str(); }
const char* what() const throw() override { return fMsg.c_str(); }
private:
std::string fMsg;
}; // DoubConvException
Expand Down
24 changes: 14 additions & 10 deletions math/genvector/inc/Math/GenVector/Cartesian2D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -89,8 +89,9 @@ public :
Scalar X() const { return fX;}
Scalar Y() const { return fY;}
Scalar Mag2() const { return fX*fX + fY*fY; }
Scalar R() const { return std::sqrt( Mag2());}
Scalar Phi() const { return (fX==0 && fY==0) ? 0.0 : atan2(fY,fX);}
Scalar R() const { using namespace std; return sqrt( Mag2()); }
Copy link
Member

@vgvassilev vgvassilev Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm, what is the advantage of the using namespace std; pattern, instead of using return std::sqrt...?

Copy link
Contributor Author

@cjones051073 cjones051073 Mar 10, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is required to support Vc types. For Vc types you need to use the versions of sqrt etc. shipped as part of that library. Using std::sqrt does not allow this. The namespace trick means you use the std:: versions when appropriate, but also allows other implementations when not.

Copy link
Member

@vgvassilev vgvassilev Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Huh, thanks for explaining! Neat and subtle at the same time. Btw, there is an inherited extra space after the open paren.

Copy link
Member

@vgvassilev vgvassilev Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Btw, it's @lmoneta's call but I'd prefer to have a short comment saying what the pattern does, or maybe a #define R__enable_vc_types using namespace std; This would make this easier to read by an unarmed eye.

Copy link
Contributor Author

@cjones051073 cjones051073 Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I dislike preprocessor directives unless absolutely required, and I don't think here adding one really helps. To my eye the pattern as is is clear, but as you say its up to @lmoneta . I would be OK with adding some sort of comment somewhere ?

Copy link
Member

@pcanal pcanal Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think a comment would be as clear and more concise than a ifdef

Copy link
Member

@lmoneta lmoneta Mar 13, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hi,

I don't understand this. Are you sure it is needed for Vc ? In the past I have seen Vc replacing automatically its vectorised function implementations in the std namespace

Copy link
Contributor Author

@cjones051073 cjones051073 Mar 13, 2017 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am sure I ran into some issue that lead me to doing this... Also note the templation is in principle there for other types than Vc. This change in principle supports types that do not extend the math functions under std::, as well as those that do.

Scalar Phi() const { using namespace std;
return (fX==Scalar(0) && fY==Scalar(0)) ? Scalar(0) : atan2(fY,fX);}

/**
set the x coordinate value keeping y constant
Expand Down Expand Up @@ -124,8 +125,9 @@ public :
rotate by an angle
*/
void Rotate(Scalar angle) {
Scalar s = std::sin(angle);
Scalar c = std::cos(angle);
using namespace std;
const Scalar s = sin(angle);
const Scalar c = cos(angle);
SetCoordinates( c*fX - s*fY, s*fX + c * fY );
}

Expand Down Expand Up @@ -161,10 +163,11 @@ public :
template <class T2>
explicit Cartesian2D( const Polar2D<T2> & v )
{
Scalar r = v.R(); // re-using this instead of calling v.X() and v.Y()
using namespace std;
const Scalar r = v.R(); // re-using this instead of calling v.X() and v.Y()
// is the speed improvement
fX = r * std::cos(v.Phi());
fY = r * std::sin(v.Phi());
fX = r * cos(v.Phi());
fY = r * sin(v.Phi());
}
// Technical note: This works even though only Polar2Dfwd.h is
// included (and in fact, including Polar2D.h would cause circularity
Expand All @@ -174,9 +177,10 @@ public :
template <class T2>
Cartesian2D & operator = (const Polar2D<T2> & v)
{
Scalar r = v.R();
fX = r * std::cos(v.Phi());
fY = r * std::sin(v.Phi());
using namespace std;
const Scalar r = v.R();
fX = r * cos(v.Phi());
fY = r * sin(v.Phi());
return *this;
}

Expand Down
33 changes: 18 additions & 15 deletions math/genvector/inc/Math/GenVector/Cartesian3D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,6 @@ public :
return *this;
}


/**
Set internal data based on an array of 3 Scalar numbers
*/
Expand All @@ -111,15 +110,17 @@ public :
Scalar Z() const { return fZ;}
Scalar Mag2() const { return fX*fX + fY*fY + fZ*fZ;}
Scalar Perp2() const { return fX*fX + fY*fY ;}
Scalar Rho() const { return std::sqrt( Perp2());}
Scalar R() const { return std::sqrt( Mag2());}
Scalar Theta() const { return (fX==0 && fY==0 && fZ==0) ?
0 : atan2(Rho(),Z());}
Scalar Phi() const { return (fX==0 && fY==0) ? 0 : atan2(fY,fX);}
Scalar Rho() const { using namespace std; return sqrt( Perp2() ); }
Scalar R() const { using namespace std; return sqrt( Mag2() ); }
Scalar Theta() const { using namespace std;
return (fX==Scalar(0) && fY==Scalar(0) && fZ==Scalar(0)) ?
Scalar(0) : atan2(Rho(),Z());}
Scalar Phi() const { using namespace std;
return (fX==Scalar(0) && fY==Scalar(0)) ? Scalar(0) : atan2(fY,fX);}

// pseudorapidity
Scalar Eta() const {
return Impl::Eta_FromRhoZ ( Rho(),fZ);
return Impl::Eta_FromRhoZ( Rho(), fZ );
}

/**
Expand Down Expand Up @@ -149,7 +150,7 @@ public :
/**
scale the vector by a scalar quantity a
*/
void Scale(Scalar a) { fX *= a; fY *= a; fZ *= a; }
void Scale(Scalar a) { fX *= a; fY *= a; fZ *= a; }

/**
negate the vector
Expand Down Expand Up @@ -190,10 +191,11 @@ public :
template <class T2>
explicit Cartesian3D( const Polar3D<T2> & v ) : fZ (v.Z())
{
T rho = v.Rho(); // re-using this instead of calling v.X() and v.Y()
// is the speed improvement
fX = rho * std::cos(v.Phi());
fY = rho * std::sin(v.Phi());
using namespace std;
const T rho = v.Rho(); // re-using this instead of calling v.X() and v.Y()
// is the speed improvement
Copy link
Member

@vgvassilev vgvassilev Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could you not break the comments here? Maybe add it as a single comment a line before?

fX = rho * cos(v.Phi());
Copy link
Member

@vgvassilev vgvassilev Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If the caching is giving the speedup, wouldn't it be better to cache v.Phi(), too?

Copy link
Contributor Author

@cjones051073 cjones051073 Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe... This is the original code and my aim was not to change anything I did not have to.

fY = rho * sin(v.Phi());
}
// Technical note: This works even though only Polar3Dfwd.h is
// included (and in fact, including Polar3D.h would cause circularity
Expand All @@ -203,9 +205,10 @@ public :
template <class T2>
Cartesian3D & operator = (const Polar3D<T2> & v)
{
T rho = v.Rho();
fX = rho * std::cos(v.Phi());
fY = rho * std::sin(v.Phi());
using namespace std;
const T rho = v.Rho();
fX = rho * cos(v.Phi());
fY = rho * sin(v.Phi());
fZ = v.Z();
return *this;
}
Expand Down
19 changes: 12 additions & 7 deletions math/genvector/inc/Math/GenVector/Cylindrical3D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -106,10 +106,12 @@ public :
{rho=fRho; zz=fZ; phi=fPhi;}

private:
inline static Scalar pi() { return M_PI; }
inline void Restrict() {
inline static Scalar pi() { return Scalar(M_PI); }
inline void Restrict()
{
using namespace std;
if ( fPhi <= -pi() || fPhi > pi() )
fPhi = fPhi - std::floor( fPhi/(2*pi()) +.5 ) * 2*pi();
fPhi = fPhi - floor( fPhi/(2*pi()) +.5 ) * 2*pi();
return;
}
public:
Expand All @@ -120,13 +122,16 @@ public :
Scalar Z() const { return fZ; }
Scalar Phi() const { return fPhi; }

Scalar X() const { return fRho*std::cos(fPhi); }
Scalar Y() const { return fRho*std::sin(fPhi); }
Scalar X() const { using namespace std; return fRho*cos(fPhi); }
Scalar Y() const { using namespace std; return fRho*sin(fPhi); }

Scalar Mag2() const { return fRho*fRho + fZ*fZ; }
Scalar R() const { return std::sqrt( Mag2()); }
Scalar R() const { using namespace std; return sqrt( Mag2()); }
Scalar Perp2() const { return fRho*fRho; }
Scalar Theta() const { return (fRho==0 && fZ==0 ) ? 0 : atan2(fRho,fZ); }
Scalar Theta() const {
using namespace std;
return (fRho==Scalar(0) && fZ==Scalar(0) ) ? Scalar(0) : atan2(fRho,fZ);
}

// pseudorapidity - use same implementation as in Cartesian3D
Scalar Eta() const {
Expand Down
94 changes: 62 additions & 32 deletions math/genvector/inc/Math/GenVector/DisplacementVector3D.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,11 +249,11 @@ namespace ROOT {
*/
template <class IT>
void GetCoordinates( IT begin) const {
Scalar a,b,c = 0;
GetCoordinates (a,b,c);
*begin++ = a;
*begin++ = b;
*begin = c;
Scalar a,b,c = Scalar(0);
Copy link
Member

@vgvassilev vgvassilev Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't that leave a and b uninitialized?

Copy link
Contributor Author

@cjones051073 cjones051073 Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Quite right.... Will fix.

Copy link
Contributor Author

@cjones051073 cjones051073 Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fix pushed.

GetCoordinates (a,b,c);
*begin++ = a;
*begin++ = b;
*begin = c;
}

/**
Expand All @@ -262,8 +262,8 @@ namespace ROOT {
then (x, y, z) are converted to that form)
*/
DisplacementVector3D<CoordSystem, Tag>& SetXYZ (Scalar a, Scalar b, Scalar c) {
fCoordinates.SetXYZ(a,b,c);
return *this;
fCoordinates.SetXYZ(a,b,c);
return *this;
}

// ------------------- Equality -----------------
Expand Down Expand Up @@ -333,13 +333,26 @@ namespace ROOT {
Scalar Perp2() const { return fCoordinates.Perp2();}

/**
return unit vector parallel to this
return unit vector parallel to this (scalar)
*/
DisplacementVector3D Unit() const {
Scalar tot = R();
template< typename SCALAR = Scalar >
Copy link
Member

@vgvassilev vgvassilev Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I personally find it clearer enabling it as a template template default arg.
template <typename SCALAR = Scalar, typename = std::enable_if< std::is_arithmetic<SCALAR>::value> and have a clear return type of the function.

Copy link
Contributor Author

@cjones051073 cjones051073 Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I prefer this one ;) So lets get a third vote.

Copy link
Member

@pcanal pcanal Mar 10, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So if I understand correctly the choice is between:

template< typename SCALAR = Scalar >
typename std::enable_if< std::is_arithmetic<SCALAR>::value, DisplacementVector3D >::type
Unit() const {

and

template <typename SCALAR = Scalar, typename = std::enable_if< std::is_arithmetic<SCALAR>::value >
DisplacementVector3 
Unit() const {

If that is correct, I prefer the 2nd one as the first one give the wrong impression that the return type is dependent on the template parameter (unless you actually read the template detail :) ) while the first clearly indicate that there is a limitation on the SCALAR type/template-parameter.

typename std::enable_if< std::is_arithmetic<SCALAR>::value, DisplacementVector3D >::type
Unit() const {
const auto tot = R();
return tot == 0 ? *this : DisplacementVector3D(*this) / tot;
}

/**
return unit vector parallel to this (scalar)
*/
template< typename SCALAR = Scalar >
typename std::enable_if< !std::is_arithmetic<SCALAR>::value, DisplacementVector3D >::type
Unit() const {
SCALAR tot = R();
tot( tot == SCALAR(0) ) = SCALAR(1);
return DisplacementVector3D(*this) / tot;
}

// ------ Setting of individual elements present in coordinate system ------

/**
Expand Down Expand Up @@ -620,35 +633,52 @@ namespace ROOT {
// ------------- I/O to/from streams -------------

template< class char_t, class traits_t, class T, class U >
inline
std::basic_ostream<char_t,traits_t> &
operator << ( std::basic_ostream<char_t,traits_t> & os
inline
typename std::enable_if< std::is_arithmetic<typename DisplacementVector3D<T,U>::Scalar>::value,
std::basic_ostream<char_t,traits_t> & >::type
operator << ( std::basic_ostream<char_t,traits_t> & os
, DisplacementVector3D<T,U> const & v
)
{
if( !os ) return os;

typename T::Scalar a, b, c;
v.GetCoordinates(a, b, c);

if( detail::get_manip( os, detail::bitforbit ) ) {
detail::set_manip( os, detail::bitforbit, '\00' );
typedef GenVector_detail::BitReproducible BR;
BR::Output(os, a);
BR::Output(os, b);
BR::Output(os, c);
if ( os ) {

typename T::Scalar a, b, c;
v.GetCoordinates(a, b, c);

if( detail::get_manip( os, detail::bitforbit ) ) {
detail::set_manip( os, detail::bitforbit, '\00' );
typedef GenVector_detail::BitReproducible BR;
BR::Output(os, a);
BR::Output(os, b);
BR::Output(os, c);
}
else {
os << detail::get_manip( os, detail::open ) << a
<< detail::get_manip( os, detail::sep ) << b
<< detail::get_manip( os, detail::sep ) << c
<< detail::get_manip( os, detail::close );
}
}
else {
os << detail::get_manip( os, detail::open ) << a
<< detail::get_manip( os, detail::sep ) << b
<< detail::get_manip( os, detail::sep ) << c
<< detail::get_manip( os, detail::close );
}

return os;

} // op<< <>()

template< class char_t, class traits_t, class T, class U >
inline
typename std::enable_if< !std::is_arithmetic<typename DisplacementVector3D<T,U>::Scalar>::value,
std::basic_ostream<char_t,traits_t> & >::type
operator << ( std::basic_ostream<char_t,traits_t> & os
, DisplacementVector3D<T,U> const & v
)
{
if ( os ) {
os << "{ ";
for ( std::size_t i = 0; i < PositionVector3D<T,U>::Scalar::Size; ++i ) {
os << "(" << v.x()[i] << "," << v.y()[i] << "," << v.z()[i] << ") ";
}
os << "}";
}
return os;
} // op<< <>()

template< class char_t, class traits_t, class T, class U >
inline
Expand Down
Loading