Skip to content

Commit da1ebd9

Browse files
devajithvsdpiparo
devajithvs
authored and
dpiparo
committed
[SYCL] Add genvectorx tutorials for SYCL
1 parent 1427667 commit da1ebd9

3 files changed

Lines changed: 251 additions & 2 deletions

File tree

‎tutorials/CMakeLists.txt‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -151,7 +151,11 @@ if(NOT davix)
151151
endif()
152152

153153
if(NOT sycl)
154-
list(APPEND sycl_veto heterogeneous/syclintro.C)
154+
list(APPEND sycl_veto heterogeneous/syclintro.C math/mathcoreGenVectorSYCL.C)
155+
endif()
156+
157+
if(NOT experimental_genvectorx)
158+
list(APPEND sycl_veto math/mathcoreGenVectorSYCL.C)
155159
endif()
156160

157161
if(NOT geom)

‎tutorials/math/index.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,7 @@ Examples showing usage of the GenVector library, generic class templates for mod
5959
| mathcoreVectorCollection.C | Example showing how to write and read a std vector of ROOT::Math LorentzVector in a ROOT tree. |
6060
| mathcoreVectorFloatIO.C | Macro illustrating I/O with Lorentz Vectors of floats. |
6161
| mathcoreVectorIO.C | Example of I/O of a GenVector Lorentz Vectors in a Tree and comparison with legacy TLorentzVector.|
62-
62+
| mathcoreGenVectorSYCL.C | Example macro testing available methods and operation of the GenVectorX classes to be used in SYCL for heterogeneous computing. |
6363

6464
\anchor histograms
6565
## Histogramming related features

‎tutorials/math/mathcoreGenVectorSYCL.C‎

Lines changed: 245 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,245 @@
1+
/// \file
2+
/// \ingroup tutorial_math
3+
/// \notebook -nodraw
4+
/// Example macro testing available methods and operation of the GenVectorX classes.
5+
///
6+
/// \note This tutorial requires ROOT to be built with **SYCL support** and **GenVectorX** enabled.
7+
/// Configure CMake with:
8+
/// `-Dexperimental_adaptivecpp=ON -Dexperimental_genvectorx=ON`
9+
///
10+
/// The results are compared and checked
11+
/// The macro is divided in 3 parts:
12+
/// - testVector3D : tests of the 3D Vector classes
13+
/// - testPoint3D : tests of the 3D Point classes
14+
/// - testLorentzVector : tests of the 4D LorentzVector classes
15+
///
16+
/// \macro_code
17+
///
18+
/// \author Devajith Valaparambil Sreeramaswamy (CERN)
19+
20+
#define ROOT_MATH_ARCH MathSYCL
21+
22+
#include "MathX/Vector3D.h"
23+
#include "MathX/Point3D.h"
24+
25+
#include "MathX/Vector2D.h"
26+
#include "MathX/Point2D.h"
27+
28+
#include "MathX/EulerAngles.h"
29+
30+
#include "MathX/Transform3D.h"
31+
#include "MathX/Translation3D.h"
32+
33+
#include "MathX/Rotation3D.h"
34+
#include "MathX/RotationX.h"
35+
#include "MathX/RotationY.h"
36+
#include "MathX/RotationZ.h"
37+
#include "MathX/Quaternion.h"
38+
#include "MathX/AxisAngle.h"
39+
#include "MathX/RotationZYX.h"
40+
41+
#include "MathX/LorentzRotation.h"
42+
#include "MathX/PtEtaPhiM4D.h"
43+
#include "MathX/LorentzVector.h"
44+
45+
#include "MathX/VectorUtil.h"
46+
47+
#include <sycl/sycl.hpp>
48+
49+
#include <vector>
50+
51+
using namespace ROOT::ROOT_MATH_ARCH;
52+
using namespace ROOT::ROOT_MATH_ARCH::VectorUtil;
53+
54+
typedef DisplacementVector3D<Cartesian3D<double>, GlobalCoordinateSystemTag> GlobalXYZVector;
55+
typedef DisplacementVector3D<Cartesian3D<double>, LocalCoordinateSystemTag> LocalXYZVector;
56+
typedef DisplacementVector3D<Polar3D<double>, GlobalCoordinateSystemTag> GlobalPolar3DVector;
57+
58+
typedef PositionVector3D<Cartesian3D<double>, GlobalCoordinateSystemTag> GlobalXYZPoint;
59+
typedef PositionVector3D<Cartesian3D<double>, LocalCoordinateSystemTag> LocalXYZPoint;
60+
typedef PositionVector3D<Polar3D<double>, GlobalCoordinateSystemTag> GlobalPolar3DPoint;
61+
typedef PositionVector3D<Polar3D<double>, LocalCoordinateSystemTag> LocalPolar3DPoint;
62+
63+
int ntest = 0;
64+
int nfail = 0;
65+
int ok = 0;
66+
67+
int compare(double v1, double v2, double scale = 1.0)
68+
{
69+
ntest = ntest + 1;
70+
71+
// numerical double limit for epsilon
72+
double eps = scale * std::numeric_limits<double>::epsilon();
73+
int iret = 0;
74+
double delta = v2 - v1;
75+
double d = 0;
76+
if (delta < 0)
77+
delta = -delta;
78+
if (v1 == 0 || v2 == 0) {
79+
if (delta > eps) {
80+
iret = 1;
81+
}
82+
}
83+
// skip case v1 or v2 is infinity
84+
else {
85+
d = v1;
86+
87+
if (v1 < 0)
88+
d = -d;
89+
// add also case when delta is small by default
90+
if (delta / d > eps && delta > eps)
91+
iret = 1;
92+
}
93+
94+
return iret;
95+
}
96+
97+
template <class Transform>
98+
bool IsEqual(const Transform &t1, const Transform &t2, unsigned int size)
99+
{
100+
// size should be an enum of the Transform class
101+
std::vector<double> x1(size);
102+
std::vector<double> x2(size);
103+
t1.GetComponents(x1.begin(), x1.end());
104+
t2.GetComponents(x2.begin(), x2.end());
105+
bool ret = true;
106+
unsigned int i = 0;
107+
while (ret && i < size) {
108+
// from TMath::AreEqualRel(x1,x2,2*eps)
109+
bool areEqual =
110+
std::abs(x1[i] - x2[i]) < std::numeric_limits<double>::epsilon() * (std::abs(x1[i]) + std::abs(x2[i]));
111+
ret &= areEqual;
112+
i++;
113+
}
114+
return ret;
115+
}
116+
117+
int testVector3D()
118+
{
119+
std::cout << "\n************************************************************************\n " << " Vector 3D Test"
120+
<< "\n************************************************************************\n";
121+
122+
sycl::buffer<int, 1> ok_buf(&ok, sycl::range<1>(1));
123+
sycl::default_selector device_selector;
124+
sycl::queue queue(device_selector);
125+
126+
{
127+
queue.submit([&](sycl::handler &cgh) {
128+
auto ok_device = ok_buf.get_access<sycl::access::mode::read_write>(cgh);
129+
cgh.single_task<class testVector3D>([=]() {
130+
// test the vector tags
131+
132+
GlobalXYZVector vg(1., 2., 3.);
133+
GlobalXYZVector vg2(vg);
134+
GlobalPolar3DVector vpg(vg);
135+
136+
ok_device[0] += compare(vpg.R(), vg2.R());
137+
138+
// std::cout << vg2 << std::endl;
139+
140+
double r = vg.Dot(vpg);
141+
ok_device[0] += compare(r, vg.Mag2());
142+
143+
GlobalXYZVector vcross = vg.Cross(vpg);
144+
ok_device[0] += compare(vcross.R(), 0.0, 10);
145+
146+
// std::cout << vg.Dot(vpg) << std::endl;
147+
// std::cout << vg.Cross(vpg) << std::endl;
148+
149+
GlobalXYZVector vg3 = vg + vpg;
150+
ok_device[0] += compare(vg3.R(), 2 * vg.R());
151+
152+
GlobalXYZVector vg4 = vg - vpg;
153+
ok_device[0] += compare(vg4.R(), 0.0, 10);
154+
});
155+
});
156+
}
157+
158+
if (ok == 0)
159+
std::cout << "\t\t OK " << std::endl;
160+
161+
return ok;
162+
}
163+
164+
int testPoint3D()
165+
{
166+
std::cout << "\n************************************************************************\n " << " Point 3D Tests"
167+
<< "\n************************************************************************\n";
168+
169+
sycl::buffer<int, 1> ok_buf(&ok, sycl::range<1>(1));
170+
sycl::default_selector device_selector;
171+
sycl::queue queue(device_selector);
172+
173+
{
174+
queue.submit([&](sycl::handler &cgh) {
175+
auto ok_device = ok_buf.get_access<sycl::access::mode::read_write>(cgh);
176+
cgh.single_task<class testPoint3D>([=]() {
177+
// test the vector tags
178+
179+
GlobalXYZPoint pg(1., 2., 3.);
180+
GlobalXYZPoint pg2(pg);
181+
182+
GlobalPolar3DPoint ppg(pg);
183+
184+
ok_device[0] += compare(ppg.R(), pg2.R());
185+
// std::cout << pg2 << std::endl;
186+
187+
GlobalXYZVector vg(pg);
188+
189+
double r = pg.Dot(vg);
190+
ok_device[0] += compare(r, pg.Mag2());
191+
192+
GlobalPolar3DVector vpg(pg);
193+
GlobalXYZPoint pcross = pg.Cross(vpg);
194+
ok_device[0] += compare(pcross.R(), 0.0, 10);
195+
196+
GlobalPolar3DPoint pg3 = ppg + vg;
197+
ok_device[0] += compare(pg3.R(), 2 * pg.R());
198+
199+
GlobalXYZVector vg4 = pg - ppg;
200+
ok_device[0] += compare(vg4.R(), 0.0, 10);
201+
202+
// operator -
203+
XYZPoint q1(1., 2., 3.);
204+
XYZPoint q2 = -1. * q1;
205+
XYZVector v2 = -XYZVector(q1);
206+
ok_device[0] += compare(XYZVector(q2) == v2, true);
207+
});
208+
});
209+
}
210+
211+
if (ok == 0)
212+
std::cout << "\t OK " << std::endl;
213+
214+
return ok;
215+
}
216+
217+
int testLorentzVector()
218+
{
219+
std::cout << "\n************************************************************************\n "
220+
<< " Lorentz Vector Tests"
221+
<< "\n************************************************************************\n";
222+
223+
LorentzVector<PtEtaPhiM4D<float>> v1(1, 2, 3, 4);
224+
LorentzVector<PtEtaPhiM4D<float>> v2(5, 6, 7, 8);
225+
ok += compare(v1.DeltaR(v2), 4.60575f);
226+
// Result cross-validated using:
227+
// TLorentzVector t1, t2;
228+
// t1.SetPtEtaPhiE(1,2,3,4); t2.SetPtEtaPhiE(5,6,7,8);
229+
// t1.DeltaR(t2)
230+
231+
if (ok == 0)
232+
std::cout << "\t OK " << std::endl;
233+
234+
return ok;
235+
}
236+
237+
void mathcoreGenVectorSYCL()
238+
{
239+
240+
testVector3D();
241+
testPoint3D();
242+
testLorentzVector();
243+
244+
std::cout << "\n\nNumber of tests " << ntest << " failed = " << nfail << std::endl;
245+
}

0 commit comments

Comments
 (0)