ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
vectorOperation.cxx
Go to the documentation of this file.
1 // test performance of all vectors operations +,- and *
2 
3 
4 
5 #include <cassert>
6 
7 #ifdef USE_VDT
8 #include "vdtMath.h"
9 #endif
10 
11 //#define USE_VC
12 #ifdef USE_VC
13 #include "Vc/Vc"
14 #include <Vc/Allocator>
15 typedef Vc::double_v Double_type;
16 #define ZERO Vc::Zero
17 #else
18 typedef double Double_type;
19 #define ZERO 0
20 #endif
21 
22 
23 #include "TRandom2.h"
24 #include "TStopwatch.h"
25 
26 #include <vector>
27 #include <iostream>
28 #include <iomanip>
29 
30 #ifdef DIM_2
31 #ifdef USE_POINT
32 #include "Math/Point2D.h"
33 typedef ROOT::Math::XYPoint VecType;
34 #elif USE_TVECTOR
35 #include "TVector2.h"
36 typedef TVector2 VecType;
37 #else
38 #include "Math/Vector2D.h"
39 typedef ROOT::Math::XYVector VecType;
40 #endif
41 
42 #ifndef USE_ROOT
43 #define VSUM(v) v.x() + v.y()
44 #else
45 #define VSUM(v) v.X() + v.Y()
46 #endif
47 
48 
49 #elif DIM_3 // 3 Dimensions
50 
51 #ifdef USE_POINT
52 #include "Math/Point3D.h"
53 typedef ROOT::Math::XYZPoint VecType;
54 #elif USE_TVECTOR
55 #include "TVector3.h"
56 typedef TVector3 VecType;
57 #else
58 #include "Math/Vector3D.h"
59 typedef ROOT::Math::XYZVector VecType;
60 #endif
61 
62 #ifndef USE_TVECTOR
63 #define VSUM(v) v.x() + v.y() + v.z()
64 #else
65 #define VSUM(v) v.X() + v.Y() + v.Z()
66 #endif
67 
68 #else // default is 4D
69 
70 #undef USE_POINT
71 #if USE_TVECTOR
72 #include "TLorentzVector.h"
73 typedef TLorentzVector VecType;
74 #else
75 #include "Math/Vector4D.h"
76 typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<Double_type> > VecType;
77 //#ifdef USE_VC
78 //Vc_DECLARE_ALLOCATOR(VecType)
79 //#endif
80 #endif
81 
82 #include "Math/VectorUtil.h"
83 
84 #ifndef USE_TVECTOR
85 #define VSUM(v) v.x() + v.y() + v.z()
86 //#define VSUM(v) v.x()
87 #else
88 #define VSUM(v) v.X() + v.Y() + v.Z()
89 #endif
90 
91 
92 #endif
93 
94 //#define VLISTSIZE 8
95 #define VLISTSIZE 100
96 
97 #ifdef USE_VC
98 const int N = VLISTSIZE/Vc::double_v::Size;
99 #else
100 const int N = VLISTSIZE;
101 #endif
102 
103 const int NLOOP = 5000000;
104 //const int NLOOP = 1;
105 
106 const int NOP = NLOOP*VLISTSIZE;
107 const double TSCALE = 1.E9/double(NOP);
108 
109 template<class Vector>
110 class TestVector {
111 public:
112 
113  TestVector();
114 
115  void Operations();
116 
117  void Add();
118  void Add2();
119  void Sub();
120  void Sub2();
121  void Scale();
122  void Scale2();
123  void Divide();
124  void Divide2();
125 
126  void MathFunction_sin();
127  void MathFunction_exp();
128  void MathFunction_log();
129  void MathFunction_atan();
130  void MathFunction_atan2();
131 
132  void Boost();
133 
134  void Read();
135 
136  void PrintSummary();
137 
138  void PrintResult(Double_type s);
139 
140 private:
141 
142  std::vector<Vector> vlist;
143  std::vector<Vector> vlist2;
144  std::vector<Double_type> scale;
145  std::vector <std::vector <double > > vcoords;
146  double fTime[50]; // timing results
147  int fTest;
148 };
149 
150 //utility function to aggregate vectors in a vc type
151 template <class Vector, class Vector_V, int Size>
152 void MakeVcVector( const Vector * vlist, Vector_V & vret ) {
153  const int dim = 4;
154  typename Vector_V::Scalar vcoord[dim];
155  typename Vector::Scalar coord[dim];
156  for (int i = 0; i < Size; ++i) {
157  vlist[i].GetCoordinates(coord);
158  for (int j = 0; j < dim; ++j)
159  vcoord[j][i] = coord[j];
160  }
161  vret.SetCootdinates(vcoord);
162  return vret;
163 }
164 // template <class Vector, class Vector_V, int Size>
165 // void UnpackVcVector( const Vector * vlist, Vector_V & vret ) {
166 // const int dim = 4;
167 // typename Vector_V::Scalar vcoord[dim];
168 // typename Vector::Scalar coord[dim];
169 // for (int i = 0; i < Size; ++i) {
170 // vlist[i].GetCoordinates(coord);
171 // for (int j = 0; j < dim; ++j)
172 // vcoord[j][i] = coord[j];
173 // }
174 // vret.SetCootdinates(vcoord);
175 // return vret;
176 // }
177 
178 
179 template<class Vector>
180 TestVector<Vector>::TestVector() :
181  vlist(N),
182  vlist2(N),
183  scale(N),
184 #ifdef USE_VC
185  vcoords(N*Vc::double_v::Size),
186 #else
187  vcoords(N),
188 #endif
189  fTest(0)
190 {
191  // create list of vectors and fill them
192 
193  TRandom2 r(111);
194 
195  double coord[4];
196  for (int i = 0; i< N; ++i) {
197 #ifndef USE_VC
198  Double_type x = r.Uniform(-1,1);
199  Double_type y = r.Uniform(-1,1);
200  Double_type z = r.Uniform(-1,1);
201  Double_type t = r.Uniform(2,10);
202  Double_type s = r.Uniform(0,1);
203  coord[0] = x; coord[1] = y; coord[2] = z; coord[3] = t;
204  vcoords[i] = std::vector<double>(coord,coord+4);
205 #else
206  Double_type x = 0.;
207  Double_type y = 0.;
208  Double_type z = 0.;
209  Double_type t = 0.;
210  Double_type s = 0.;
211  for (int j = 0; j< Vc::double_v::Size; ++j) {
212  x[j] = r.Uniform(-1,1);
213  y[j] = r.Uniform(-1,1);
214  z[j] = r.Uniform(-1,1);
215  t[j] = r.Uniform(2,10);
216  s[j] = r.Uniform(0,1);
217  coord[0] = x[j]; coord[1] = y[j]; coord[2] = z[j]; coord[3] = t[j];
218  vcoords[i*Vc::double_v::Size+j] = std::vector<double>(coord,coord+4);
219  }
220 #endif
221 
222 #ifdef DIM_2
223  vlist[i] = Vector( x, y );
224 #elif DIM_3
225  vlist[i] = Vector( x, y, z);
226 #else // 4D
227  vlist[i] = Vector( x, y, z, t);
228 #endif
229  scale[i] = s;
230  }
231 
232  std::cout << "test using " << typeid(vlist[0]).name() << std::endl;
233  std::cout << "Vector used " << vlist[0] << std::endl;
234  std::cout << "Vector used " << vlist[1] << std::endl;
235 
236  // create second list of vectors which is same vector shifted by 1
237  for (int i = 0; i< N; ++i) {
238 #ifndef USE_VC
239  vlist2[i] = (i < N-1) ? vlist[i+1] : vlist[0];
240 #else
241  Double_type x1 = vlist[i].X();
242  Double_type y1 = vlist[i].Y();
243  Double_type z1 = vlist[i].Z();
244  Double_type t1 = vlist[i].E();
245  Double_type x2 = (i< N-1) ? vlist[i+1].X() : vlist[0].X();
246  Double_type y2 = (i< N-1) ? vlist[i+1].Y() : vlist[0].Y();
247  Double_type z2 = (i< N-1) ? vlist[i+1].Z() : vlist[0].Z();
248  Double_type t2 = (i< N-1) ? vlist[i+1].E() : vlist[0].E();
249  Double_type x;
250  Double_type y;
251  Double_type z;
252  Double_type t;
253  int j = 0;
254  for (j = 0; j< Vc::double_v::Size-1; ++j) {
255  x[j] = x1[j+1];
256  y[j] = y1[j+1];
257  z[j] = z1[j+1];
258  t[j] = t1[j+1];
259  }
260  j = Vc::double_v::Size-1;
261  x[j] = x2[0];
262  y[j] = y2[0];
263  z[j] = z2[0];
264  t[j] = t2[0];
265  vlist2[i] = Vector( x, y, z, t);
266 #endif
267  }
268 
269 }
270 
271 
272 template<class Vector>
273 void TestVector<Vector>::PrintResult(Double_type s)
274  // print result
275 {
276 #ifndef USE_VC
277  std::cout << "value " << s << std::endl << std::endl;
278 #else
279  Double_t s2 = 0;
280  for (int i = 0; i < Vc::double_v::Size; ++i)
281  s2 += s[i];
282 // std::cout << "s = " << s << " sum ";
283  std::cout << "value " << s2 << std::endl << std::endl;
284 #endif
285 }
286 
287 template<class Vector>
288 void TestVector<Vector>::Read()
289  // just read vector
290 {
291  TStopwatch w;
292  w.Start();
293  Double_type s(0.0);
294  for (int l = 0; l<NLOOP; ++l) {
295  //Vector v0 = vlist[l%N];
296  //s += v0.X();
297  //s = 0;
298  for (int i = 0; i< N; ++i) {
299  Vector v3 = vlist[i];
300  s += VSUM(v3);
301  // if (l == 0) {
302  // std::cout << v3 << " sum " << VSUM(v3) << " total sum " << s << std::endl;
303  // }
304  }
305  }
306 
307  std::cout << "Time for read v3 :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
308  PrintResult(s);
309  fTime[fTest++] = w.CpuTime();
310 }
311 
312 
313 template<class Vector>
314 void TestVector<Vector>::Add()
315  // normal addition
316 {
317  TStopwatch w;
318  w.Start();
319  Double_type s(0.0);
320  for (int l = 0; l<NLOOP; ++l) {
321  for (int i = 0; i< N; ++i) {
322  Vector v3 = vlist[i] + vlist2[i];
323  s += VSUM(v3);
324 // std::cout << vlist[i] << " " << vlist2[i] << std::endl;
325  //s += v3.X();
326 // PrintResult(s);
327  }
328  }
329 
330  std::cout << "Time for v3 = v1 + v2 :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
331  PrintResult(s);
332  fTime[fTest++] = w.CpuTime()*TSCALE;
333 }
334 
335 template<class Vector>
336 void TestVector<Vector>::Add2()
337 {
338  // self addition
339  TStopwatch w;
340  Vector v3;
341  w.Start();
342  Double_type s(0.0);
343  for (int l = 0; l<NLOOP; ++l) {
344  v3 = Vector();
345  for (int i = 0; i< N; ++i) {
346  v3 += vlist[i];
347  }
348  // if I put s inside inner loop break autivec
349  // but also Vc cannot work (will compute something different)
350  //s += v3.X(); (if I use this compiler can auto-vectorize)
351  s += VSUM(v3);
352  }
353 
354  std::cout << "Time for v2 += v1 :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
355  PrintResult(s);
356  fTime[fTest++] = w.CpuTime()*TSCALE;
357 }
358 
359 template<class Vector>
360 void TestVector<Vector>::Sub()
361 {
362  // normal sub
363  TStopwatch w;
364  w.Start();
365  Double_type s(0.0);
366  for (int l = 0; l<NLOOP; ++l) {
367  for (int i = 0; i< N; ++i) {
368  Vector v3 = vlist[i] - vlist2[i];
369  s += VSUM(v3);
370  }
371  }
372 
373  std::cout << "Time for v3 = v1 - v2 :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
374  PrintResult(s);
375  fTime[fTest++] = w.CpuTime()*TSCALE;
376 }
377 
378 template<class Vector>
379 void TestVector<Vector>::Sub2()
380 {
381  // self subtruction
382  TStopwatch w;
383  Vector v3;
384  w.Start();
385  Double_type s(0.0);
386  for (int l = 0; l<NLOOP; ++l) {
387  for (int i = 0; i< N; ++i) {
388  v3 -= vlist[i];
389  s += VSUM(v3);
390  }
391  }
392 
393  std::cout << "Time for v2 -= v1 :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
394  PrintResult(s);
395  fTime[fTest++] = w.CpuTime()*TSCALE;
396 }
397 
398 
399 template<class Vector>
400 void TestVector<Vector>::Scale()
401 {
402 // normal multiply
403  TStopwatch w;
404  w.Start();
405  Double_type s(0.0);
406  for (int l = 0; l<NLOOP; ++l) {
407  for (int i = 0; i< N; ++i) {
408  Vector v3 = scale[i]*vlist[i];
409  s += VSUM(v3);
410  }
411  }
412 
413  std::cout << "Time for v2 = A * v1 :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
414  PrintResult(s);
415  fTime[fTest++] = w.CpuTime()*TSCALE;
416 }
417 
418 template<class Vector>
419 void TestVector<Vector>::Scale2()
420 {
421  // self scale
422  TStopwatch w;
423  Vector v3;
424  w.Start();
425  Double_type s(0.0);
426  for (int l = 0; l<NLOOP; ++l) {
427  for (int i = 0; i< N; ++i) {
428  v3 = vlist[i];
429  v3 *= scale[i];
430  s += VSUM(v3);
431  }
432  }
433 
434  std::cout << "Time for v *= a :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
435  PrintResult(s);
436  fTime[fTest++] = w.CpuTime()*TSCALE;
437 }
438 
439 template<class Vector>
440 void TestVector<Vector>::Divide()
441 {
442 // normal divide
443  TStopwatch w;
444  w.Start();
445  Double_type s(0.0);
446  for (int l = 0; l<NLOOP; ++l) {
447  for (int i = 0; i< N; ++i) {
448  Vector v3 = vlist[i]/scale[i];
449  s += VSUM(v3);
450  }
451  }
452 
453  std::cout << "Time for v2 = v1 / a :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
454  PrintResult(s);
455  fTime[fTest++] = w.CpuTime()*TSCALE;
456 }
457 
458 template<class Vector>
459 void TestVector<Vector>::Divide2()
460 {
461  // self divide
462  TStopwatch w;
463  Vector v3;
464  w.Start();
465  Double_type s(0.0);
466  for (int l = 0; l<NLOOP; ++l) {
467  for (int i = 0; i< N; ++i) {
468  v3 = vlist[i];
469  v3 /= scale[i];
470  s += VSUM(v3);
471  }
472  }
473 
474  std::cout << "Time for v /= a :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
475  PrintResult(s);
476  fTime[fTest++] = w.CpuTime()*TSCALE;
477 }
478 
479 #ifdef CASE1
480 template<class Vector>
481 void TestVector<Vector>::Operations()
482 {
483  // test several operations
484  TStopwatch w;
485  Vector v1;
486  Vector v2;
487  Vector v3;
488  Vector v4;
489  w.Start();
490  Double_type s(0.0);
491  const int n = sqrt(N)+0.5;
492  for (int l = 0; l<NLOOP; ++l) {
493 #ifndef USE_VC
494  for (int i = 0; i< n; ++i) {
495  for (int j = 0; j< i; ++j) {
496  s += ROOT::Math::VectorUtil::InvariantMass(vlist[i], vlist[j] );
497  //std::cout << "inv mass of " << vlist[i] << " " << vlist[j] << " is " << s << std::endl;
498  }
499  //v2 = vlist[i]*scale[i];
500  // v1 *= scale[i-1];
501  // v2 /= scale[i];
502  // v3 = v1 + v2;
503  // v4 = v1 - v2;
504  //s += ROOT::Math::VectorUtil::InvariantMass2(v3,v4);
505  //s += ROOT::Math::VectorUtil::InvariantMass2(v1,v2);
506  //s += std::sin(v1.E() );
507  }
508 #else
509  const int nn = sqrt(vcoords.size()) + 0.5;
510  int ncomb = nn*(nn-1)/2;
511  std::vector<double> vc1x(ncomb+4);
512  std::vector<double> vc2x(ncomb+4);
513  std::vector<double> vc1y(ncomb+4);
514  std::vector<double> vc2y(ncomb+4);
515  std::vector<double> vc1z(ncomb+4);
516  std::vector<double> vc2z(ncomb+4);
517  std::vector<double> vc1t(ncomb+4);
518  std::vector<double> vc2t(ncomb+4);
519  double c1[4];
520  double c2[4];
521  int k = 0;
522  for (int i = 0; i< nn; ++i) {
523  std::copy(vcoords[i].begin(), vcoords[i].end(),c1);
524  //std::cout << "vcoord " << vcoords[i][0] << " " << c1[0] << std::endl;
525  for (int j = 0; j< i; ++j) {
526  std::copy(vcoords[j].begin(), vcoords[j].end(),c2);
527  vc1x[k] = c1[0];
528  vc2x[k] = c2[0];
529  vc1y[k] = c1[1];
530  vc2y[k] = c2[1];
531  vc1z[k] = c1[2];
532  vc2z[k] = c2[2];
533  vc1t[k] = c1[3];
534  vc2t[k] = c2[3];
535  k++;
536  }
537  }
538  int ncomb2 = ncomb/Vc::double_v::Size;
539  if (ncomb%Vc::Double_v::Size != 0) ncomb2 += 1;
540  Vector v1; Vector v2;
541  for (int i = 0; i< ncomb2; ++i) {
542 
543  typename Vector::Scalar cv[4];
544  cv[0].load( &vc1x[i*Vc::double_v::Size], Vc::Unaligned );
545  cv[1].load( &vc1y[i*Vc::double_v::Size], Vc::Unaligned );
546  cv[2].load( &vc1z[i*Vc::double_v::Size], Vc::Unaligned );
547  cv[3].load( &vc1t[i*Vc::double_v::Size], Vc::Unaligned );
548  //std::cout << cv[0] << " " << vc1x[i*Vc::double_v::Size] << std::endl;
549  v1.SetCoordinates( cv);
550 
551  typename Vector::Scalar cv2[4];
552  cv2[0].load( &vc2x[i*Vc::double_v::Size], Vc::Unaligned );
553  cv2[1].load( &vc2y[i*Vc::double_v::Size], Vc::Unaligned );
554  cv2[2].load( &vc2z[i*Vc::double_v::Size], Vc::Unaligned );
555  cv2[3].load( &vc2t[i*Vc::double_v::Size], Vc::Unaligned );
556  //std::cout << cv[0] << " " << vc2x[i*Vc::double_v::Size] << std::endl;
557  v2.SetCoordinates( cv2);
558 
559  s+= ROOT::Math::VectorUtil::InvariantMass(v1, v2);
560  //std::cout << "inv mass of " << v1 << " " << v2 << " is " << s << std::endl;
561  }
562 
563 #endif
564  }
565 
566  std::cout << "Time for Operation :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
567  PrintResult(s);
568  fTime[fTest++] = w.CpuTime()*TSCALE;
569 }
570 
571 #else
572 
573 template<class Vector>
574 void TestVector<Vector>::Operations()
575 {
576  // test several operations
577  TStopwatch w;
578  Vector v1;
579  Vector v2;
580  Vector v3;
581  Vector v4;
582  w.Start();
583  Double_type s(0.0);
584 
585 #ifndef USE_VC
586  const int n = sqrt(N)+0.5;
587 #endif
588 
589  const int nn = sqrt(vcoords.size()) + 0.5;
590  int ncomb = nn*(nn-1)/2;
591  std::vector<double> vc1x(ncomb+4);
592  std::vector<double> vc2x(ncomb+4);
593  std::vector<double> vc1y(ncomb+4);
594  std::vector<double> vc2y(ncomb+4);
595  std::vector<double> vc1z(ncomb+4);
596  std::vector<double> vc2z(ncomb+4);
597  std::vector<double> vc1t(ncomb+4);
598  std::vector<double> vc2t(ncomb+4);
599  double c1[4];
600  double c2[4];
601  int k = 0;
602  for (int i = 0; i< nn; ++i) {
603  std::copy(vcoords[i].begin(), vcoords[i].end(),c1);
604  //std::cout << "vcoord " << vcoords[i][0] << " " << c1[0] << std::endl;
605  for (int j = 0; j< i; ++j) {
606  std::copy(vcoords[j].begin(), vcoords[j].end(),c2);
607  vc1x[k] = c1[0];
608  vc2x[k] = c2[0];
609  vc1y[k] = c1[1];
610  vc2y[k] = c2[1];
611  vc1z[k] = c1[2];
612  vc2z[k] = c2[2];
613  vc1t[k] = c1[3];
614  vc2t[k] = c2[3];
615  k++;
616  }
617  }
618 
619 
620 
621  for (int l = 0; l<NLOOP; ++l) {
622 #ifndef USE_VC
623  for (int i = 0; i< n; ++i) {
624  for (int j = 0; j< i; ++j) {
625  s += ROOT::Math::VectorUtil::InvariantMass(vlist[i], vlist[j] );
626  //std::cout << "inv mass of " << vlist[i] << " " << vlist[j] << " is " << s << std::endl;
627  }
628  //v2 = vlist[i]*scale[i];
629  // v1 *= scale[i-1];
630  // v2 /= scale[i];
631  // v3 = v1 + v2;
632  // v4 = v1 - v2;
633  //s += ROOT::Math::VectorUtil::InvariantMass2(v3,v4);
634  //s += ROOT::Math::VectorUtil::InvariantMass2(v1,v2);
635  //s += std::sin(v1.E() );
636  }
637 #else
638  int ncomb2 = ncomb/Vc::double_v::Size;
639  if (ncomb%Vc::double_v::Size != 0) ncomb2 += 1;
640  Vector v1; Vector v2;
641  for (int i = 0; i< ncomb2; ++i) {
642 
643  typename Vector::Scalar cv[4];
644  cv[0].load( &vc1x[i*Vc::double_v::Size], Vc::Unaligned );
645  cv[1].load( &vc1y[i*Vc::double_v::Size], Vc::Unaligned );
646  cv[2].load( &vc1z[i*Vc::double_v::Size], Vc::Unaligned );
647  cv[3].load( &vc1t[i*Vc::double_v::Size], Vc::Unaligned );
648  //std::cout << cv[0] << " " << vc1x[i*Vc::double_v::Size] << std::endl;
649  v1.SetCoordinates( cv);
650 
651  typename Vector::Scalar cv2[4];
652  cv2[0].load( &vc2x[i*Vc::double_v::Size], Vc::Unaligned );
653  cv2[1].load( &vc2y[i*Vc::double_v::Size], Vc::Unaligned );
654  cv2[2].load( &vc2z[i*Vc::double_v::Size], Vc::Unaligned );
655  cv2[3].load( &vc2t[i*Vc::double_v::Size], Vc::Unaligned );
656  //std::cout << cv[0] << " " << vc2x[i*Vc::double_v::Size] << std::endl;
657  v2.SetCoordinates( cv2);
658 
659  s+= ROOT::Math::VectorUtil::InvariantMass(v1, v2);
660  // std::cout << s << std::endl;
661  // PrintResult(s);
662  //std::cout << "inv mass of " << v1 << " " << v2 << " is " << s << std::endl;
663  }
664 
665  // std::cout << " \n" << std::endl;
666  // std::cout << "End trial " << l << std::endl;
667  // PrintResult(s);
668  // std::cout << " \n" << std::endl;
669 #endif
670  }
671 
672  std::cout << "Time for Operation :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
673  PrintResult(s);
674  fTime[fTest++] = w.CpuTime()*TSCALE;
675 }
676 #endif
677 
678 
679 template<class Vector>
680 void TestVector<Vector>::Boost()
681 {
682  // test several operations
683  TStopwatch w;
684  Vector v1;
685  Vector v2;
686  w.Start();
687  Double_type s(0.0);
688  for (int l = 0; l<NLOOP; ++l) {
689  for (int i = 0; i< N; ++i) {
690  v1 = vlist[i];
691  v2 = ROOT::Math::VectorUtil::boostX(v1,scale[i]);
692  // s += Vc::atan2(v1.Y(),v1.X()) - Vc::atan2(v2.Y(),v2.X());
693  //s += VSUM(v2);
694  }
695  }
696 
697  std::cout << "Time for Boost :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
698  PrintResult(s);
699  fTime[fTest++] = w.CpuTime()*TSCALE;
700 }
701 
702 
703 
704 template<class Vector>
705 void TestVector<Vector>::MathFunction_exp()
706 {
707  // test math function
708  TStopwatch w;
709  Vector v1;
710  w.Start();
711  Double_type s(0.0);
712  for (int l = 0; l<NLOOP/10; ++l) {
713  for (int i = 0; i< N; ++i) {
714  v1 = vlist[i];
715  s += std::exp( v1.X() - v1.Y() );
716  }
717  }
718 
719  std::cout << "Time for Exp Function :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
720  PrintResult(s);
721  fTime[fTest++] = w.CpuTime()*TSCALE*10;
722 }
723 
724 
725 
726 template<class Vector>
727 void TestVector<Vector>::MathFunction_log()
728 {
729  // test several operations
730  TStopwatch w;
731  Vector v1;
732  w.Start();
733  Double_type s(0.0);
734  for (int l = 0; l<NLOOP/10; ++l) {
735  for (int i = 0; i< N; ++i) {
736  v1 = vlist[i];
737  s += std::log( std::abs(v1.X() ) );
738  }
739  }
740 
741  std::cout << "Time for Log Function :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
742  PrintResult(s);
743  fTime[fTest++] = w.CpuTime()*TSCALE*10;
744 }
745 
746 template<class Vector>
747 void TestVector<Vector>::MathFunction_sin()
748 {
749  // test math function
750  TStopwatch w;
751  Vector v1;
752  w.Start();
753  Double_type s(0.0);
754  for (int l = 0; l<NLOOP/100; ++l) {
755  for (int i = 0; i< N; ++i) {
756  v1 = vlist[i];
757  s += std::sin( v1.X() );
758  }
759  }
760 
761  std::cout << "Time for Sin Function :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
762  PrintResult(s);
763  fTime[fTest++] = w.CpuTime()*TSCALE*100;
764 }
765 
766 template<class Vector>
767 void TestVector<Vector>::MathFunction_atan()
768 {
769  // test several operations
770  TStopwatch w;
771  Vector v1;
772  w.Start();
773  Double_type s(0.0);
774  for (int l = 0; l<NLOOP/100; ++l) {
775  for (int i = 0; i< N; ++i) {
776  v1 = vlist[i]; ///scale[i];
777  s += std::atan( v1.Y()/v1.X() );
778  }
779  }
780 
781  std::cout << "Time for Atan Function :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
782  PrintResult(s);
783  fTime[fTest++] = w.CpuTime()*TSCALE*100;
784 }
785 
786 template<class Vector>
787 void TestVector<Vector>::MathFunction_atan2()
788 {
789  // test several operations
790  TStopwatch w;
791  Vector v1;
792  w.Start();
793  Double_type s(0.0);
794  for (int l = 0; l<NLOOP/100; ++l) {
795  for (int i = 0; i< N; ++i) {
796  v1 = vlist[i]; ///scale[i];
797  s += std::atan2( v1.Y(),v1.X() );
798  }
799  }
800 
801  std::cout << "Time for Atan2 Function :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
802  PrintResult(s);
803  fTime[fTest++] = w.CpuTime()*TSCALE*100;
804 }
805 
806 #ifdef HAS_VDT // use VDT
807 template<class Vector>
808 void TestVector<Vector>::MathFunction()
809 {
810  // test several operations
811  TStopwatch w;
812  Vector v;
813  double x[N];
814  double r[N];
815  w.Start();
816  Double_type s(0.0);
817  for (int l = 0; l<NLOOP; ++l) {
818  for (int i = 0; i< N; ++i) {
819  Vector v = vlist[i];
820  //x[i] = v.X()/v.Pt() ;
821  x[i] = v.X()/v.Y();
822  }
823  vdt::fast_atanv(N,x,r);
824  for (int i = 0; i< N; ++i) {
825  s += x[i];
826  }
827 
828  //s += sin(v1.E() );
829  }
830 
831  std::cout << "Time for MathFUnction(VDT) :\t" << w.RealTime() << "\t" << w.CpuTime() << std::endl;
832  PrintResult(s);
833  fTime[fTest++] = w.CpuTime()*TSCALE;
834 }
835 
836 #endif
837 
838 template<class Vector>
839 void TestVector<Vector>::PrintSummary()
840 {
841  std::cout << "\nResults for " << typeid(vlist[0]).name() << std::endl;
842  std::cout << " v3 = v1+v2"
843  << " v2 += v1 "
844  << " v3 = v1-v2"
845  << " v2 -= v1 "
846  << " v2 = a*v1 "
847  << " v1 *= a "
848  << " v2 = v1/a "
849  << " v1 /= a "
850  << " log "
851  << " exp "
852  << " sin "
853  << " atan "
854  << " atan2 "
855  << std::endl;
856 
857  // start from 3
858  for (int i = 3; i < fTest; ++i) {
859  std::cout << std::setw(8) << fTime[i] << " ";
860  }
861  std::cout << std::endl << std::endl;
862 }
863 
864 int main() {
865  TestVector<VecType> t;
866 
867 #ifdef USE_VC
868  std::cout << "testing using Vc: Vc size is " << Vc::double_v::Size << " Looping on " << N << " vectors" << std::endl;
869  std::cout << "Implementation type: " << VC_IMPL << std::endl;
870 #else
871  std::cout << "testing using standard double's. Looping on " << N << " vectors" << std::endl;
872 #endif
873 
874 
875  t.Read();
876 
877  t.Operations();
878  t.Boost();
879 
880 
881 #ifndef USE_POINT
882  t.Add();
883  t.Add2();
884  t.Sub();
885  t.Sub2();
886 #endif
887  t.Scale();
888  t.Scale2();
889 #ifndef USE_TVECTOR
890  t.Divide();
891  t.Divide2();
892 #endif
893 
894 
895  t.MathFunction_log();
896  t.MathFunction_exp();
897  t.MathFunction_sin();
898  t.MathFunction_atan();
899  t.MathFunction_atan2();
900 
901 
902 
903 
904 
905  // summurize test
906  t.PrintSummary();
907 }
ROOT::Math::LorentzVector
Class describing a generic LorentzVector in the 4D space-time, using the specified coordinate system ...
Definition: LorentzVector.h:54
VC_IMPL
#define VC_IMPL
Definition: global.h:429
TStopwatch::RealTime
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:108
TStopwatch::Start
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
v1
const Double_t * v1
Definition: TArcBall.cxx:33
Add
TMatrixT< Element > & Add(TMatrixT< Element > &target, Element scalar, const TMatrixT< Element > &source)
Modify addition: target += scalar * source.
Definition: TMatrixT.cxx:2925
xmlio::Size
const char * Size
Definition: TXMLSetup.cxx:56
TRandom2
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition: TRandom2.h:29
ROOT::Math::VectorUtil::boostX
LVector boostX(const LVector &v, T beta)
Boost a generic Lorentz Vector class along the X direction with a factor beta The only requirement on...
Definition: VectorUtil.h:358
VSUM
#define VSUM(v)
Definition: vectorOperation.cxx:85
ROOT::Math::PositionVector3D
Class describing a generic position vector (point) in 3 dimensions.
Definition: PositionVector3D.h:63
TStopwatch::CpuTime
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:123
main
int main()
Definition: vectorOperation.cxx:330
NLOOP
const int NLOOP
Definition: vectorOperation.cxx:103
sqrt
double sqrt(double)
x
Double_t x[n]
Definition: legend1.C:17
ROOT::Vc::AVX::abs
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
Definition: vector.h:450
ROOT::Math::VectorUtil::InvariantMass
Vector1::Scalar InvariantMass(const Vector1 &v1, const Vector2 &v2)
return the invariant mass of two LorentzVector The only requirement on the LorentzVector is that they...
Definition: VectorUtil.h:217
N
const int N
Definition: vectorOperation.cxx:99
sin
double sin(double)
z
Float_t z[5]
Definition: Ifit.C:16
NOP
const int NOP
Definition: vectorOperation.cxx:106
t
TThread * t[5]
Definition: threadsh1.C:13
MakeVcVector
void MakeVcVector(const Vector *vlist, Vector_V &vret)
Definition: vectorOperation.cxx:152
r
ROOT::R::TRInterface & r
Definition: Object.C:4
ROOT::Vc::AVX::double_v
Vector< double > double_v
Definition: vector.h:416
VLISTSIZE
#define VLISTSIZE
Definition: vectorOperation.cxx:95
v
SVector< double, 2 > v
Definition: Dict.h:5
ROOT::Vc::double_v
VECTOR_NAMESPACE::double_v double_v
Definition: vector.h:80
Scalar
#define Scalar
Definition: global.h:83
TSCALE
const double TSCALE
Definition: vectorOperation.cxx:107
TMath::E
Double_t E()
Definition: TMath.h:54
qtexample.w
tuple w
Definition: qtexample.py:51
l
TLine * l
Definition: textangle.C:4
ROOT::Math::PositionVector2D
Class describing a generic position vector (point) in 2 dimensions.
Definition: PositionVector2D.h:56
Double_type
double Double_type
Definition: vectorOperation.cxx:18
Double_t
double Double_t
Definition: RtypesCore.h:55
atan2
double atan2(double, double)
y
Double_t y[n]
Definition: legend1.C:17
atan
double atan(double)
VecType
ROOT::Math::XYZTVector VecType
Definition: vectorOperation.cxx:87
name
#define name(a, b)
Definition: linkTestLib0.cpp:5
else
else
Definition: TBase64.cxx:55
Z
std::complex< float_v > Z
Definition: main.cpp:120
exp
double exp(double)
n
const Int_t n
Definition: legend1.C:16
log
double log(double)
ROOT::Math::DisplacementVector2D
Class describing a generic displacement vector in 2 dimensions.
Definition: DisplacementVector2D.h:66
TStopwatch
Stopwatch class.
Definition: TStopwatch.h:30