|
| 1 | +#include "MathX/Vector3D.h" |
| 2 | +#include "MathX/Point3D.h" |
| 3 | + |
| 4 | +#include "MathX/Vector2D.h" |
| 5 | +#include "MathX/Point2D.h" |
| 6 | + |
| 7 | +#include "MathX/EulerAngles.h" |
| 8 | + |
| 9 | +#include "MathX/Transform3D.h" |
| 10 | +#include "MathX/Translation3D.h" |
| 11 | + |
| 12 | +#include "MathX/Rotation3D.h" |
| 13 | +#include "MathX/RotationX.h" |
| 14 | +#include "MathX/RotationY.h" |
| 15 | +#include "MathX/RotationZ.h" |
| 16 | +#include "MathX/Quaternion.h" |
| 17 | +#include "MathX/AxisAngle.h" |
| 18 | +#include "MathX/RotationZYX.h" |
| 19 | + |
| 20 | +#include "MathX/LorentzRotation.h" |
| 21 | +#include "MathX/PtEtaPhiM4D.h" |
| 22 | +#include "MathX/LorentzVector.h" |
| 23 | + |
| 24 | +#include "MathX/VectorUtil.h" |
| 25 | + |
| 26 | +#include "MathX/GenVectorX/AccHeaders.h" |
| 27 | + |
| 28 | +#include <sycl/sycl.hpp> |
| 29 | + |
| 30 | +#include <vector> |
| 31 | + |
| 32 | +using namespace ROOT::ROOT_MATH_ARCH; |
| 33 | +using namespace ROOT::ROOT_MATH_ARCH::VectorUtil; |
| 34 | + |
| 35 | +typedef LorentzVector<PtEtaPhiM4D< double>> vec4d; |
| 36 | + |
| 37 | +int compare(double v1, double v2, double scale = 1.0) |
| 38 | +{ |
| 39 | + // ntest = ntest + 1; |
| 40 | + |
| 41 | + // numerical double limit for epsilon |
| 42 | + double eps = scale * std::numeric_limits<double>::epsilon(); |
| 43 | + int iret = 0; |
| 44 | + double delta = v2 - v1; |
| 45 | + double d = 0; |
| 46 | + if (delta < 0) |
| 47 | + delta = -delta; |
| 48 | + if (v1 == 0 || v2 == 0) { |
| 49 | + if (delta > eps) { |
| 50 | + iret = 1; |
| 51 | + } |
| 52 | + } |
| 53 | + // skip case v1 or v2 is infinity |
| 54 | + else { |
| 55 | + d = v1; |
| 56 | + |
| 57 | + if (v1 < 0) |
| 58 | + d = -d; |
| 59 | + // add also case when delta is small by default |
| 60 | + if (delta / d > eps && delta > eps) |
| 61 | + iret = 1; |
| 62 | + } |
| 63 | + |
| 64 | + // if (iret == 0) |
| 65 | + // std::cout << "."; |
| 66 | + // else { |
| 67 | + // int pr = std::cout.precision(18); |
| 68 | + // std::cout << "\nDiscrepancy in " << name << "() : " << v1 << " != " << v2 << " discr = " << int(delta / d / |
| 69 | + // eps) |
| 70 | + // << " (Allowed discrepancy is " << eps << ")\n"; |
| 71 | + // std::cout.precision(pr); |
| 72 | + // // nfail = nfail + 1; |
| 73 | + // } |
| 74 | + return iret; |
| 75 | +} |
| 76 | + |
| 77 | +template <class Transform> |
| 78 | +bool IsEqual(const Transform &t1, const Transform &t2, unsigned int size) |
| 79 | +{ |
| 80 | + // size should be an enum of the Transform class |
| 81 | + std::vector<double> x1(size); |
| 82 | + std::vector<double> x2(size); |
| 83 | + t1.GetComponents(x1.begin(), x1.end()); |
| 84 | + t2.GetComponents(x2.begin(), x2.end()); |
| 85 | + bool ret = true; |
| 86 | + unsigned int i = 0; |
| 87 | + while (ret && i < size) { |
| 88 | + // from TMath::AreEqualRel(x1,x2,2*eps) |
| 89 | + bool areEqual = |
| 90 | + std::abs(x1[i] - x2[i]) < std::numeric_limits<double>::epsilon() * (std::abs(x1[i]) + std::abs(x2[i])); |
| 91 | + ret &= areEqual; |
| 92 | + i++; |
| 93 | + } |
| 94 | + return ret; |
| 95 | +} |
| 96 | + |
| 97 | +vec4d *GenVectors(int n) |
| 98 | +{ |
| 99 | + vec4d *vectors = new vec4d[n]; |
| 100 | + |
| 101 | + // generate n -4 momentum quantities |
| 102 | + for (int i = 0; i < n; ++i) { |
| 103 | + // fill vectors |
| 104 | + vectors[i] = {1., 1., 1., 1.}; |
| 105 | + } |
| 106 | + |
| 107 | + return vectors; |
| 108 | +} |
| 109 | + |
| 110 | +int testInvariantMasses(int N) |
| 111 | +{ |
| 112 | + int iret_host = 0; |
| 113 | + std::cout << "testing Invariant Masses Computation \t:\n"; |
| 114 | + |
| 115 | + sycl::default_selector device_selector; |
| 116 | + sycl::queue queue(device_selector); |
| 117 | + |
| 118 | + auto v1 = GenVectors(N); |
| 119 | + auto v2 = GenVectors(N); |
| 120 | + double *invMasses = new double[N]; |
| 121 | + |
| 122 | + std::cout << "sycl::queue check - selected device:\n" |
| 123 | + << queue.get_device().get_info<sycl::info::device::name>() << std::endl; |
| 124 | + |
| 125 | + { |
| 126 | + |
| 127 | + vec4d *d_v1 = sycl::malloc_device<vec4d>(N, queue); |
| 128 | + vec4d *d_v2 = sycl::malloc_device<vec4d>(N, queue); |
| 129 | + double *d_invMasses = sycl::malloc_device<double>(N, queue); |
| 130 | + |
| 131 | + queue.memcpy(d_v1, v1, N * sizeof(vec4d)); |
| 132 | + queue.memcpy(d_v2, v2, N * sizeof(vec4d)); |
| 133 | + queue.wait(); |
| 134 | + |
| 135 | + queue.submit([&](sycl::handler &cgh) { |
| 136 | + cgh.parallel_for(sycl::range<1>(N), [=](sycl::id<1> indx) { |
| 137 | + vec4d v = d_v2[indx] + d_v2[indx]; |
| 138 | + d_invMasses[indx] = v.M(); |
| 139 | + }); |
| 140 | + }); |
| 141 | + |
| 142 | + queue.wait(); |
| 143 | + queue.memcpy(invMasses, d_invMasses, N * sizeof(double)); |
| 144 | + queue.wait(); |
| 145 | + } |
| 146 | + |
| 147 | + for (int i = 0; i < N; i++) { |
| 148 | + iret_host += (std::abs(invMasses[i] - 2.) > 1e-5); |
| 149 | + } |
| 150 | + |
| 151 | + if (iret_host == 0) |
| 152 | + std::cout << "\tOK\n"; |
| 153 | + else |
| 154 | + std::cout << "\t FAILED\n"; |
| 155 | + |
| 156 | + return iret_host; |
| 157 | +} |
| 158 | + |
| 159 | +int main() |
| 160 | +{ |
| 161 | + int n = 128; |
| 162 | + int ret = testInvariantMasses(n); |
| 163 | + if (ret) |
| 164 | + std::cerr << "test FAILED !!! " << std::endl; |
| 165 | + else |
| 166 | + std::cout << "test OK " << std::endl; |
| 167 | + return ret; |
| 168 | +} |
0 commit comments