ROOT  6.06/09
Reference Guide
testOperations.cxx
Go to the documentation of this file.
1 #define ENABLE_TEMPORARIES
2 
3 
4 #ifdef USE_VC
5 #include "Vc/Vc"
6 #endif
7 
8 #include "Math/SVector.h"
9 #include "Math/SMatrix.h"
10 
11 
12 #include "TMatrixD.h"
13 #include "TVectorD.h"
14 
15 #include "TRandom3.h"
16 #include "TH1D.h"
17 #include "TProfile.h"
18 #include "TFile.h"
19 
20 //#define HAVE_CLHEP
21 #define TEST_SYM
22 
23 #ifdef TEST_ALL_MATRIX_SIZES
24 #define REPORT_TIME
25 #endif
26 #ifndef NITER
27 #define NITER 1 // number of iterations
28 #endif
29 #ifndef NLOOP_MIN
30 #define NLOOP_MIN 100;
31 #endif
32 
33 
34 #ifdef HAVE_CLHEP
35 #include "CLHEP/Matrix/SymMatrix.h"
36 #include "CLHEP/Matrix/Matrix.h"
37 #include "CLHEP/Matrix/Vector.h"
38 #endif
39 
40 //#define DEBUG
41 
42 #include <iostream>
43 
44 // #ifndef NDIM1
45 // #define NDIM1 5
46 // #endif
47 // #ifndef NDIM2
48 // #define NDIM2 5
49 // #endif
50 
51 
52 int NLOOP;
53 
54 
55 //#define NLOOP 1
56 
57 //#define DEBUG
58 
59 #define NLISTSIZE 64 // size of matrix/vector lists
60 
61 
62 #ifdef USE_VC
63 typedef Vc::double_v Stype;
64 const int NLIST = NLISTSIZE / Vc::double_v::Size;
65 #else
66 typedef double Stype;
67 const int NLIST = NLISTSIZE;
68 #endif
69 
70 
71 Stype gResultSum;
72 
73 
74 #include "matrix_op_vec.h"
75 #include "matrix_util.h"
76 #include <map>
77 
78 
79 
80 template<unsigned int ND1, unsigned int ND2 >
81 int test_smatrix_op() {
82 
83  // need to write explicitly the dimensions
84 
85 
86  typedef SMatrix<Stype, ND1, ND1> MnMatrixNN;
87  typedef SMatrix<Stype, ND2, ND2> MnMatrixMM;
88  typedef SMatrix<Stype, ND1, ND2> MnMatrixNM;
89  typedef SMatrix<Stype, ND2 , ND1> MnMatrixMN;
90  typedef SVector<Stype, ND1> MnVectorN;
91  typedef SVector<Stype, ND2> MnVectorM;
92 
93 
94 
95  int first = ND1; //Can change the size of the matrices
96  int second = ND2;
97 
98 
99  std::cout << "************************************************\n";
100  std::cout << " SMatrix operations test " << first << " x " << second << std::endl;
101  std::cout << "************************************************\n";
102 
103 
104  double t_veq, t_meq, t_vad, t_mad, t_dot, t_mv, t_gmv, t_mm, t_prd, t_inv, t_vsc, t_msc, t_ama, t_tra;
105  double totTime1, totTime2;
106 
107 
108 
109  int npass = NITER;
110  TRandom3 r(111);
111 
112  Stype r1[NLIST];
113  Stype r2[NLIST];
114 
115  gResultSum = 0;
116 
117  for (int k = 0; k < npass; k++) {
118 
119 
120  MnMatrixNM A[NLIST];
121  MnMatrixMN B[NLIST];
122  MnMatrixNN C[NLIST];
123  MnMatrixMM D[NLIST];
124  MnVectorN v[NLIST];
125  MnVectorM p[NLIST];
126 
127  TStopwatch w;
128 
129  // fill matrices with random data
130  for (int j = 0; j < NLIST; j++) fillRandomMat(r,A[j],first,second);
131  for (int j = 0; j < NLIST; j++) fillRandomMat(r,B[j],second,first);
132  for (int j = 0; j < NLIST; j++) fillRandomMat(r,C[j],first,first);
133  for (int j = 0; j < NLIST; j++) fillRandomMat(r,D[j],second,second);
134  for (int j = 0; j < NLIST; j++) fillRandomVec(r,v[j],first);
135  for (int j = 0; j < NLIST; j++) fillRandomVec(r,p[j],second);
136 
137 
138 
139 #ifdef DEBUG
140  std::cout << "pass " << k << std::endl;
141  if (k == 0) {
142  std::cout << " A = " << A[0] << std::endl;
143  std::cout << " B = " << B[0] << std::endl;
144  std::cout << " C = " << C[0] << std::endl;
145  std::cout << " D = " << D[0] << std::endl;
146  std::cout << " v = " << v[0] << std::endl;
147  std::cout << " p = " << p[0] << std::endl;
148  }
149 #endif
150 
151  w.Start();
152 
153 
154 
155 
156  MnVectorN v1[NLIST]; testMV(A,v,t_mv,v1);
157  //if (k == 0) v1.Print(std::cout);
158  MnVectorN v2[NLIST]; testGMV(A,v,v1,t_gmv,v2);
159  //if (k == 0) v2.Print(std::cout);
160  MnMatrixNN C0[NLIST]; testMM(A,B,C,t_mm,C0);
161  //if (k == 0) C0.Print(std::cout);
162  MnMatrixNN C1[NLIST]; testATBA_S(B,C0,t_ama,C1);
163  //if (k == 0) C1.Print(std::cout);
164 
165  MnMatrixNN C2[NLIST]; testInv_S(C1,t_inv,C2);
166  MnMatrixNN C3[NLIST]; testInvFast_S(C2,t_inv,C3);
167  //MnMatrixNN C2 = C1;
168 
169  MnVectorN v3[NLIST]; testVeq(v,t_veq,v3);
170  MnVectorN v4[NLIST]; testVad(v2,v3,t_vad,v4);
171  MnVectorN v5[NLIST]; testVscale(v4,2.0,t_vsc,v5);
172  MnMatrixNN C4[NLIST]; testMeq(C,t_meq,C4);
173  MnMatrixNN C5[NLIST]; testMad(C3,C4,t_mad,C5);
174  MnMatrixNN C6[NLIST]; testMscale(C5,0.5,t_msc,C6);
175  MnMatrixNN C7[NLIST]; testMT_S(C6,t_tra,C7);
176 
177 #ifdef DEBUG
178  if (k == 0) {
179  std::cout << " C6 = " << C5 << std::endl;
180  std::cout << " v5 = " << v5 << std::endl;
181  }
182 #endif
183 
184  testDot_S(v3,v5,r1,t_dot);
185 
186  testInnerProd_S(C7,v5,r2,t_prd);
187 
188 
189  w.Stop();
190  totTime1 = w.RealTime();
191  totTime2 = w.CpuTime();
192 
193  }
194  //tr.dump();
195 
196  //double totTime = t_veq + t_meq + t_vad + t_mad + t_dot + t_mv + t_gmv + t_mm + t_prd + t_inv + t_vsc + t_msc + t_ama + t_tra;
197  std::cout << "Total Time = " << totTime1 << " (s) " << " cpu " << totTime2 << " (s) " << std::endl;
198  std::cerr << "SMatrix: r1[0] = " << r1[0] << " r2[0] = " << r2[0] << std::endl;
199  std::cerr << "SMatrix: r1[N] = " << r1[NLIST-1] << " r2[N] = " << r2[NLIST-1] << std::endl;
200  std::cerr << "sum of results = " << gResultSum << std::endl;
201 
202  return 0;
203 }
204 
205 
206 
207 #ifdef TEST_SYM
208 template<unsigned int ND1, unsigned int ND2>
209 int test_smatrix_sym_op() {
210 
211  // need to write explicitly the dimensions
212 
213 
214  typedef SMatrix<Stype, ND1, ND1, MatRepSym<Stype,ND1> > MnSymMatrixNN;
215  typedef SMatrix<Stype, ND1, ND1 > MnMatrixNN;
216  typedef SVector<Stype, ND1> MnVectorN;
217 
218 
219 
220  int first = ND1; //Can change the size of the matrices
221 
222 
223  std::cout << "************************************************\n";
224  std::cout << " SMatrixSym operations test " << first << " x " << first << std::endl;
225  std::cout << "************************************************\n";
226 
227 
228  double t_meq, t_mad, t_mv, t_gmv, t_mm, t_prd, t_inv, t_msc, t_ama = 0;
229  double totTime1, totTime2;
230 
231 
232 
233  Stype r1[NLIST];
234  int npass = NITER;
235  TRandom3 r(111);
236  gResultSum = 0;
237 
238  for (int k = 0; k < npass; k++) {
239 
240 
241  MnSymMatrixNN A[NLIST];
242  MnSymMatrixNN B[NLIST];
243  MnMatrixNN C[NLIST];
244  MnVectorN v[NLIST];
245 
246 
247  TStopwatch w;
248 
249  // fill matrices with random data
250  for (int j = 0; j < NLIST; j++) fillRandomSym(r,A[j],first);
251  for (int j = 0; j < NLIST; j++) fillRandomSym(r,B[j],first);
252  for (int j = 0; j < NLIST; j++) fillRandomMat(r,C[j],first,first);
253  for (int j = 0; j < NLIST; j++) fillRandomVec(r,v[j],first);
254 
255 
256 
257 #ifdef DEBUG
258  std::cout << "pass " << k << std::endl;
259  if (k == 0) {
260  std::cout << " A = " << A[0] << std::endl;
261  std::cout << " B = " << B[0] << std::endl;
262  std::cout << " C = " << C[0] << std::endl;
263  std::cout << " v = " << v[0] << std::endl;
264  }
265 #endif
266 
267  w.Start();
268 
269  MnVectorN v1[NLIST]; testMV(A,v,t_mv,v1);
270  MnVectorN v2[NLIST]; testGMV(A,v,v1,t_gmv,v2);
271  MnMatrixNN C0[NLIST]; testMM(A,B,C,t_mm,C0);
272  MnSymMatrixNN C1[NLIST]; testATBA_S2(C0,B,t_ama,C1);
273  MnSymMatrixNN C2[NLIST]; testInv_S(A,t_inv,C2);
274  MnSymMatrixNN C3[NLIST]; testInvFast_S(C2,t_inv,C3);
275  MnSymMatrixNN C4[NLIST]; testInvChol_S(C3,t_inv,C4);
276  //C2 = C1;
277  MnSymMatrixNN C5[NLIST]; testMeq(C4,t_meq,C5);
278  MnSymMatrixNN C6[NLIST]; testMad(A,C5,t_mad,C6);
279  MnSymMatrixNN C7[NLIST]; testMscale(C6,0.5,t_msc,C7);
280 
281  testInnerProd_S(C7,v2,r1,t_prd);
282 
283 
284 #ifdef DEBUG
285  std::cout << "output matrices" << std::endl;
286  if (k == 0) {
287  std::cout << " C1 = " << C1 << std::endl;
288  std::cout << " C3 = " << C3 << std::endl;
289  std::cout << " C4 = " << C4 << std::endl;
290  std::cout << " C5 = " << C5 << std::endl;
291  }
292 #endif
293 
294 
295  w.Stop();
296  totTime1 = w.RealTime();
297  totTime2 = w.CpuTime();
298 
299 
300  }
301  //tr.dump();
302 
303  //double totTime = t_meq + t_mv + t_gmv + t_mm + t_prd + t_inv + t_mad + t_msc + t_ama;
304  std::cout << "Total Time = " << totTime1 << " (s) - cpu " << totTime2 << " (s) " << std::endl;
305  std::cerr << "SMatrixSym: r1[0] = " << r1[0] << std::endl;
306  std::cerr << "SMatrixSym: r1[N] = " << r1[NLIST-1] << std::endl;
307  std::cerr << "sum of results = " << gResultSum << std::endl;
308 
309  return 0;
310 }
311 #endif
312 
313 
314 // ROOT test
315 
316 
317 template<unsigned int ND1, unsigned int ND2>
318 int test_tmatrix_op() {
319 
320 #ifdef USE_TMATRIX
321 
322 
323 
324  typedef TMatrixD MnMatrix;
325  typedef TVectorD MnVector;
326 
327 // typedef boost::numeric::ublas::matrix<double> MnMatrix;
328  //typedef HepSymMatrix MnSymMatrixHep;
329 
330 
331  int first = ND1; //Can change the size of the matrices
332  int second = ND2;
333 
334 
335  std::cout << "************************************************\n";
336  std::cout << " TMatrix operations test " << first << " x " << second << std::endl;
337  std::cout << "************************************************\n";
338 
339  double t_veq, t_meq, t_vad, t_mad, t_dot, t_mv, t_gmv, t_mm, t_prd, t_inv, t_vsc, t_msc, t_ama, t_tra = 0;
340  double totTime1, totTime2;
341 
342  double r1,r2;
343  int npass = NITER;
344  TRandom3 r(111);
345  gMatrixCheck = 0;
346 
347  for (int k = 0; k < npass; k++) {
348 
349 
350  MnMatrix A(ND1,ND2);
351  MnMatrix B(ND2,ND1);
352  MnMatrix C(ND1,ND1);
353  MnMatrix D(ND2,ND2);
354  MnVector v(ND1);
355  MnVector p(ND2);
356 
357 
358  TStopwatch w;
359  {
360  // fill matrices with random data
361  fillRandomMat(r,A,first,second);
362  fillRandomMat(r,B,second,first);
363  fillRandomMat(r,C,first,first);
364  fillRandomMat(r,D,second,second);
365 
366  fillRandomVec(r,v,first);
367  fillRandomVec(r,p,second);
368 
369  }
370 
371 #ifdef DEBUG
372  std::cout << "pass " << k << std::endl;
373  if (k == 0) {
374  A.Print(); B.Print(); C.Print(); D.Print(); v.Print(); p.Print();
375  }
376 #endif
377  w.Start();
378 
379 
380  MnVector v1(ND1); testMV_T(A,v,t_mv,v1);
381  //if (k == 0) v1.Print();
382  MnVector v2(ND1); testGMV_T(A,v,v1,t_gmv,v2);
383  //if (k == 0) v2.Print();
384  MnMatrix C0(ND1,ND1); testMM_T(A,B,C,t_mm,C0);
385  //if (k == 0) C0.Print();
386  MnMatrix C1(ND1,ND1); testATBA_T(B,C0,t_ama,C1);
387  //if (k == 0) C1.Print();
388  MnMatrix C2(ND1,ND1); testInv_T(C1,t_inv,C2);
389  //if (k == 0) C2.Print();
390  MnVector v3(ND1); testVeq(v,t_veq,v3);
391  MnVector v4(ND1); testVad_T(v2,v3,t_vad,v4);
392  MnVector v5(ND1); testVscale_T(v4,2.0,t_vsc,v5);
393  MnMatrix C3(ND1,ND1); testMeq(C,t_meq,C3);
394  MnMatrix C4(ND1,ND1); testMad_T(C2,C3,t_mad,C4);
395  //if (k == 0) C4.Print();
396  MnMatrix C5(ND1,ND1); testMscale_T(C4,0.5,t_msc,C5);
397  //if (k == 0) C5.Print();
398  MnMatrix C6(ND1,ND1); testMT_T(C5,t_tra,C6);
399 
400 #ifdef DEBUG
401  if (k == 0) {
402  C6.Print();
403  v5.Print();
404  }
405 #endif
406 
407  r1 = testDot_T(v3,v5,t_dot);
408 
409  r2 = testInnerProd_T(C6,v5,t_prd);
410 
411  //MnMatrix C2b(ND1,ND1); testInv_T2(C1,t_inv2,C2b);
412 
413 
414  w.Stop();
415  totTime1 = w.RealTime();
416  totTime2 = w.CpuTime();
417  }
418  // tr.dump();
419 
420  //double totTime = t_veq + t_meq + t_vad + t_mad + t_dot + t_mv + t_gmv + t_mm + t_prd + t_inv + t_inv2 + t_vsc + t_msc + t_ama + t_tra;
421  std::cout << "Total Time = " << totTime1 << " (s) - cpu " << totTime2 << " (s) " << std::endl;
422  std::cerr << "TMatrix: r1 = " << r1 << " r2 = " << r2 << std::endl;
423 
424 #endif
425 
426  return 0;
427 
428 }
429 
430 
431 
432 #ifdef TEST_SYM
433 template<unsigned int ND1, unsigned int ND2>
434 int test_tmatrix_sym_op() {
435 
436 #ifdef USE_TMATRIX
437 
438  // need to write explicitly the dimensions
439 
440 
441  typedef TMatrixDSym MnSymMatrix;
442  typedef TMatrixD MnMatrix;
443  typedef TVectorD MnVector;
444 
445 
446 
447  int first = ND1; //Can change the size of the matrices
448 
449 
450  std::cout << "************************************************\n";
451  std::cout << " TMatrixSym operations test " << first << " x " << first << std::endl;
452  std::cout << "************************************************\n";
453 
454 
455  double t_meq, t_mad, t_mv, t_gmv, t_mm, t_prd, t_inv, t_msc, t_ama = 0;
456  double totTime1, totTime2;
457 
458 
459 
460  double r1;
461  int npass = NITER;
462  TRandom3 r(111);
463  for (int k = 0; k < npass; k++) {
464 
465 
466  MnSymMatrix A(ND1);
467  MnSymMatrix B(ND1);
468  MnMatrix C(ND1,ND1);
469  MnVector v(ND1);
470 #define N ND1
471 
472  TStopwatch w;
473 
474  {
475  // fill matrices with random data
476  fillRandomSym(r,A,first);
477  fillRandomSym(r,B,first);
478  fillRandomMat(r,C,first,first);
479 
480  fillRandomVec(r,v,first);
481 
482  }
483 
484 
485 #ifdef DEBUG
486  std::cout << "pass " << k << std::endl;
487  if (k == 0) {
488  A.Print(); B.Print(); C.Print(); v.Print();
489  }
490 #endif
491 
492  w.Start();
493 
494  MnVector v1(N); testMV_T(A,v,t_mv,v1);
495  MnVector v2(N); testGMV_T(A,v,v1,t_gmv,v2);
496  MnMatrix C0(N,N); testMM_T(A,B,C,t_mm,C0);
497  MnSymMatrix C1(N); testATBA_T2(C0,B,t_ama,C1);
498  MnSymMatrix C2(N); testInv_T(A,t_inv,C2);
499  MnSymMatrix C3(N); testMeq(C2,t_meq,C3);
500  MnSymMatrix C4(N); testMad_T(A,C3,t_mad,C4);
501  MnSymMatrix C5(N); testMscale_T(C4,0.5,t_msc,C5);
502 
503  r1 = testInnerProd_T(C5,v2,t_prd);
504 
505 #ifdef DEBUG
506  std::cout << "output matrices" << std::endl;
507  if (k == 0) {
508  C1.Print(); C3.Print(); C4.Print(); C5.Print();
509  }
510 #endif
511 
512  w.Stop();
513  totTime1 = w.RealTime();
514  totTime2 = w.CpuTime();
515 
516  }
517  //tr.dump();
518 
519  //double totTime = t_meq + t_mv + t_gmv + t_mm + t_prd + t_inv + t_mad + t_msc + t_ama;
520  std::cout << "Total Time = " << totTime1 << " (s) - cpu " << totTime2 << " (s) " << std::endl;
521  std::cerr << "TMatrixSym: r1 = " << r1 << std::endl;
522 
523 #endif
524 
525  return 0;
526 }
527 #endif // end TEST_SYM
528 
529 #ifdef HAVE_CLHEP
530 
531 template<unsigned int ND1, unsigned int ND2>
532 int test_hepmatrix_op() {
533 
534 
535 
536 
537  typedef HepMatrix MnMatrix;
538  typedef HepVector MnVector;
539 
540 
541 
542  int first = ND1; //Can change the size of the matrices
543  int second = ND2;
544 
545 
546  std::cout << "************************************************\n";
547  std::cout << " HepMatrix operations test " << first << " x " << second << std::endl;
548  std::cout << "************************************************\n";
549 
550  double t_veq, t_meq, t_vad, t_mad, t_dot, t_mv, t_gmv, t_mm, t_prd, t_inv, t_vsc, t_msc, t_ama, t_tra = 0;
551 
552 
553  double totTime1, totTime2;
554 
555  double r1,r2;
556  int npass = NITER;
557  TRandom3 r(111);
558 
559  for (int k = 0; k < npass; k++) {
560 
561 
562  MnMatrix A(ND1,ND2);
563  MnMatrix B(ND2,ND1);
564  MnMatrix C(ND1,ND1);
565  MnMatrix D(ND2,ND2);
566  MnVector v(ND1);
567  MnVector p(ND2);
568 
569  TStopwatch w;
570 
571  {
572  // fill matrices with random data
573  fillRandomMat(r,A,first,second,1);
574  fillRandomMat(r,B,second,first,1);
575  fillRandomMat(r,C,first,first,1);
576  fillRandomMat(r,D,second,second,1);
577 
578  fillRandomVec(r,v,first);
579  fillRandomVec(r,p,second);
580  }
581 
582 #ifdef DEBUG
583  std::cout << "pass " << k << std::endl;
584  if (k == 0) {
585  std::cout << " A = " << A << std::endl;
586  std::cout << " B = " << B << std::endl;
587  std::cout << " C = " << C << std::endl;
588  std::cout << " D = " << D << std::endl;
589  std::cout << " v = " << v << std::endl;
590  std::cout << " p = " << p << std::endl;
591  }
592 #endif
593 
594  w.Start();
595 
596  MnVector v1(ND1); testMV(A,v,t_mv,v1);
597  MnVector v2(ND1); testGMV(A,v,v1,t_gmv,v2);
598  MnMatrix C0(ND1,ND1); testMM_C(A,B,C,t_mm,C0);
599  MnMatrix C1(ND1,ND1); testATBA_C(B,C0,t_ama,C1);
600  //std::cout << " C1 = " << C1 << std::endl;
601  MnMatrix C2(ND1,ND1); testInv_C(C1,t_inv,C2);
602  //std::cout << " C2 = " << C2 << std::endl;
603  MnVector v3(ND1); testVeq(v,t_veq,v3);
604  MnVector v4(ND1); testVad(v2,v3,t_vad,v4);
605  MnVector v5(ND1); testVscale(v4,2.0,t_vsc,v5);
606  MnMatrix C3(ND1,ND1); testMeq_C(C,t_meq,C3);
607  MnMatrix C4(ND1,ND1); testMad_C(C2,C3,t_mad,C4);
608  //std::cout << " C4 = " << C4 << std::endl;
609  MnMatrix C5(ND1,ND1); testMscale_C(C4,0.5,t_msc,C5);
610  //std::cout << " C5 = " << C5 << std::endl;
611  MnMatrix C6(ND1,ND1); testMT_C(C5,t_tra,C6);
612 
613 
614  r1 = testDot_C(v3,v5,t_dot);
615  r2 = testInnerProd_C(C6,v5,t_prd);
616 
617 #ifdef DEBUG
618  if (k == 0) {
619  std::cout << " C6 = " << C6 << std::endl;
620  std::cout << " v5 = " << v5 << std::endl;
621  }
622 #endif
623 
624  // MnMatrix C2b(ND1,ND1); testInv_T2(C1,t_inv2,C2b);
625 
626  w.Stop();
627  totTime1 = w.RealTime();
628  totTime2 = w.CpuTime();
629 
630  }
631  // tr.dump();
632 
633  std::cout << "Total Time = " << totTime1 << " (s) - cpu " << totTime2 << " (s) " << std::endl;
634  std::cerr << "HepMatrix: r1 = " << r1 << " r2 = " << r2 << std::endl;
635 
636  return 0;
637 }
638 
639 
640 #ifdef TEST_SYM
641 template<unsigned int ND1, unsigned int ND2>
642 int test_hepmatrix_sym_op() {
643 
644  // need to write explicitly the dimensions
645 
646 
647  typedef HepSymMatrix MnSymMatrix;
648  typedef HepMatrix MnMatrix;
649  typedef HepVector MnVector;
650 
651 
652 
653  int first = ND1; //Can change the size of the matrices
654 
655 
656  std::cout << "************************************************\n";
657  std::cout << " HepMatrixSym operations test " << first << " x " << first << std::endl;
658  std::cout << "************************************************\n";
659 
660 
661  double t_meq, t_mad, t_mv, t_gmv, t_mm, t_prd, t_inv, t_msc, t_ama = 0;
662 
663  double totTime1, totTime2;
664 
665 
666  double r1;
667  int npass = NITER;
668  TRandom3 r(111);
669  for (int k = 0; k < npass; k++) {
670 
671 
672  MnSymMatrix A(ND1);
673  MnSymMatrix B(ND1);
674  MnMatrix C(ND1,ND1);
675  MnVector v(ND1);
676 #define N ND1
677 
678  TStopwatch w;
679 
680  {
681  // fill matrices with random data
682  fillRandomSym(r,A,first,1);
683  fillRandomSym(r,B,first,1);
684  fillRandomMat(r,C,first,first,1);
685  fillRandomVec(r,v,first);
686 
687  }
688 
689 
690 #ifdef DEBUG
691  std::cout << "pass " << k << std::endl;
692  if (k == 0) {
693  }
694 #endif
695 
696  w.Start();
697 
698  MnVector v1(N); testMV(A,v,t_mv,v1);
699  MnVector v2(N); testGMV(A,v,v1,t_gmv,v2);
700  MnMatrix C0(N,N); testMM_C(A,B,C,t_mm,C0);
701  MnSymMatrix C1(N); testATBA_C2(C0,B,t_ama,C1);
702  MnSymMatrix C2(N); testInv_C(A,t_inv,C2);
703  MnSymMatrix C3(N); testMeq_C(C2,t_meq,C3);
704  MnSymMatrix C4(N); testMad_C(A,C3,t_mad,C4);
705  MnSymMatrix C5(N); testMscale_C(C4,0.5,t_msc,C5);
706 
707  r1 = testInnerProd_C(C5,v2,t_prd);
708 
709 #ifdef DEBUG
710  std::cout << "output matrices" << std::endl;
711  if (k == 0) {
712  }
713 #endif
714 
715  w.Stop();
716  totTime1 = w.RealTime();
717  totTime2 = w.CpuTime();
718 
719  }
720  //tr.dump();
721 
722  std::cout << "Total Time = " << totTime1 << " (s) - cpu " << totTime2 << " (s) " << std::endl;
723  std::cerr << "HepMatrixSym: r1 = " << r1 << std::endl;
724 
725  return 0;
726 }
727 
728 #endif // TEST_SYM
729 #endif // HAVE_CLHEP
730 
731 
732 #if defined(HAVE_CLHEP) && defined (TEST_SYM)
733 #define NTYPES 6
734 #define TEST(N) \
735  MATRIX_SIZE=N; \
736  TEST_TYPE=0; test_smatrix_op<N,N>(); \
737  TEST_TYPE=1; test_tmatrix_op<N,N>(); \
738  TEST_TYPE=2; test_hepmatrix_op<N,N>(); \
739  TEST_TYPE=3; test_smatrix_sym_op<N,N>(); \
740  TEST_TYPE=4; test_tmatrix_sym_op<N,N>(); \
741  TEST_TYPE=5; test_hepmatrix_sym_op<N,N>();
742 #elif !defined(HAVE_CLHEP) && !defined(USE_TMATRIX) && defined (TEST_SYM)
743 #define NTYPES 2
744 #define TEST(N) \
745  MATRIX_SIZE=N; \
746  TEST_TYPE=0; test_smatrix_op<N,N>(); \
747  TEST_TYPE=1; test_smatrix_sym_op<N,N>();
748 #elif !defined(HAVE_CLHEP) && defined(USE_TMATRIX) && defined (TEST_SYM)
749 #define NTYPES 4
750 #define TEST(N) \
751  MATRIX_SIZE=N; \
752  TEST_TYPE=0; test_smatrix_op<N,N>(); \
753  TEST_TYPE=1; test_tmatrix_op<N,N>(); \
754  TEST_TYPE=2; test_smatrix_sym_op<N,N>(); \
755  TEST_TYPE=3; test_tmatrix_sym_op<N,N>();
756 #elif defined(HAVE_CLHEP) && !defined (TEST_SYM)
757 #define NTYPES 3
758 #define TEST(N) \
759  MATRIX_SIZE=N; \
760  TEST_TYPE=0; test_smatrix_op<N,N>(); \
761  TEST_TYPE=1; test_tmatrix_op<N,N>(); \
762  TEST_TYPE=2; test_hepmatrix_op<N,N>();
763 #else
764 #define NTYPES 2
765 #define TEST(N) \
766  TEST_TYPE=0; test_smatrix_op<N,N>(); \
767  TEST_TYPE=1; test_tmatrix_op<N,N>();
768 #endif
769 
770 
771 
772 int TEST_TYPE;
773 int MATRIX_SIZE;
774 #ifdef REPORT_TIME
775 std::vector< std::map<std::string, TH1D *> > testTimeResults(NTYPES);
776 std::vector< std::string > typeNames(NTYPES);
777 
778 void ROOT::Math::test::reportTime(std::string s, double time) {
779  assert( TEST_TYPE >= 0 && TEST_TYPE < NTYPES );
780  std::map<std::string, TH1D * > & result = testTimeResults[TEST_TYPE];
781 
782  std::map<std::string, TH1D * >::iterator pos = result.find(s);
783  TH1D * h = 0;
784  if ( pos != result.end() ) {
785  h = pos->second;
786  }
787  else {
788  // add new elements in map
789  //std::cerr << "insert element in map" << s << typeNames[TEST_TYPE] << std::endl;
790  std::string name = typeNames[TEST_TYPE] + "_" + s;
791  h = new TProfile(name.c_str(), name.c_str(),100,0.5,100.5);
792  //result.insert(std::map<std::string, TH1D * >::value_type(s,h) );
793  result[s] = h;
794  }
795  double scale=1;
796  if (s.find("dot") != std::string::npos ||
797  s.find("V=V") != std::string::npos ||
798  s.find("V+V") != std::string::npos ) scale = 10;
799  h->Fill(double(MATRIX_SIZE),time/double(NLOOP*NITER*scale) );
800 }
801 #endif
802 
803 int testOperations() {
804 
805  NLOOP = 1000*NLOOP_MIN;
806 
807 #ifdef USE_VC
808  std::cout << "Using VC library - size = " << Vc::double_v::Size << " VC_IMPL = " << VC_IMPL << std::endl;
809 #endif
810 
811  std::cout << " making vector/matrix lists of size = " << NLIST << std::endl;
812 
813 
814  initValues();
815 
816 #ifndef NDIM1
817 #define NDIM1 5
818 #endif
819  TEST(NDIM1)
820 
821  return 0;
822 }
823 
824 
825 int main(int argc , char *argv[] ) {
826 
827 
828  std::string fname = "testOperations";
829  if (argc > 1) {
830  std::string platf(argv[1]);
831  fname = fname + "_" + platf;
832  }
833  fname = fname + ".root";
834 
835 
836 #ifdef REPORT_TIME
837  TFile * f = new TFile(fname.c_str(),"recreate");
838 
839  typeNames[0] = "SMatrix";
840  typeNames[1] = "TMatrix";
841 #if !defined(HAVE_CLHEP) && defined (TEST_SYM)
842  typeNames[2] = "SMatrix_sym";
843  typeNames[3] = "TMatrix_sym";
844 #elif defined(HAVE_CLHEP) && defined (TEST_SYM)
845  typeNames[2] = "HepMatrix";
846  typeNames[3] = "SMatrix_sym";
847  typeNames[4] = "TMatrix_sym";
848  typeNames[5] = "HepMatrix_sym";
849 #elif defined(HAVE_CLHEP) && !defined (TEST_SYM)
850  typeNames[2] = "HepMatrix";
851 #endif
852 
853 #endif
854 
855 #ifndef TEST_ALL_MATRIX_SIZES
856 // NLOOP = 1000*NLOOP_MIN
857 // initValues();
858 
859 // TEST(5)
860 // NLOOP = 50*NLOOP_MIN;
861 // TEST(30);
862 
863  return testOperations();
864 
865 #else
866  NLOOP = 5000*NLOOP_MIN;
867  initValues();
868 
869 
870 
871  TEST(2);
872  TEST(3);
873  TEST(4);
874  NLOOP = 1000*NLOOP_MIN
875  TEST(5);
876  TEST(6);
877  TEST(7);
878  TEST(10);
879  NLOOP = 100*NLOOP_MIN;
880  TEST(15);
881  TEST(20);
882  NLOOP = 50*NLOOP_MIN;
883  TEST(30);
884 // NLOOP = NLOOP_MIN;
885 // TEST(50);
886 // TEST(75);
887 // TEST(100);
888 #endif
889 
890 #ifdef REPORT_TIME
891  f->Write();
892  f->Close();
893 #endif
894 
895 }
896 
testATBA_C2
void testATBA_C2(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:550
test_tmatrix_sym_op
int test_tmatrix_sym_op()
Definition: testOperations.cxx:403
npass
const int npass
Definition: testPermute.cxx:21
NLIST
const int NLIST
Definition: testOperations.cxx:67
testVeq
void testVeq(const V &v, double &time, V &result)
Definition: matrix_op.h:33
ROOT::Math::Cephes::B
static double B[]
Definition: SpecFuncCephes.cxx:178
TH1::Fill
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3159
testATBA_T
void testATBA_T(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:311
TRandom3
Random number generator class based on M.
Definition: TRandom3.h:29
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
testMM
void testMM(const A &a, const B &b, const C &c, double &time, C &result)
Definition: matrix_op.h:137
xmlio::Size
const char * Size
Definition: TXMLSetup.cxx:56
assert
#define assert(cond)
Definition: unittest.h:542
testMV_T
void testMV_T(const M &mat, const V &v, double &time, V &result)
Definition: matrix_op.h:262
h
TH1 * h
Definition: legend2.C:5
TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:45
testVscale
void testVscale(const V &v1, double a, double &time, V &result)
Definition: matrix_op.h:85
testATBA_T2
void testATBA_T2(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:401
fillRandomVec
void fillRandomVec(TRandom &r, V &v, unsigned int len, unsigned int start=0, double offset=1)
Definition: matrix_util.h:7
fillRandomSym
void fillRandomSym(TRandom &r, M &m, unsigned int first, unsigned int start=0, double offset=1)
Definition: matrix_util.h:20
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
testMscale
void testMscale(const M &m1, double a, double &time, M &result)
Definition: matrix_op.h:98
TProfile
Profile Historam.
Definition: TProfile.h:34
gResultSum
Stype gResultSum
Definition: testOperations.cxx:71
ROOT::Math::Cephes::A
static double A[]
Definition: SpecFuncCephes.cxx:170
testDot_S
double testDot_S(const V &v1, const V &v2, double &time)
Definition: matrix_op.h:154
TStopwatch::Stop
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
ROOT::Math::SMatrix
SMatrix: a generic fixed size D1 x D2 Matrix class.
Definition: BinaryOperators.h:32
NDIM1
#define NDIM1
testMM_T
void testMM_T(const A &a, const B &b, const C &c, double &time, C &result)
Definition: matrix_op.h:287
testInv_T
void testInv_T(const M &m, double &time, M &result)
Definition: matrix_op.h:350
test_tmatrix_op
int test_tmatrix_op()
Definition: testOperations.cxx:290
testMeq
void testMeq(const M &m, double &time, M &result)
Definition: matrix_op.h:45
initValues
void initValues()
Definition: matrix_op.h:18
testMad
void testMad(const M &m1, const M &m2, double &time, M &result)
Definition: matrix_op.h:71
testOperations
int testOperations()
Definition: testOperations.cxx:762
testInnerProd_C
double testInnerProd_C(const M &a, const V &v, double &time)
Definition: matrix_op.h:456
testInv_S
void testInv_S(const M &m, double &time, M &result)
Definition: matrix_op.h:194
NLOOP_MIN
#define NLOOP_MIN
Definition: testOperations.cxx:30
ROOT::Math::Cephes::C
static double C[]
Definition: SpecFuncCephes.cxx:187
testMscale_T
void testMscale_T(const M &m1, double a, double &time, M &result)
Definition: matrix_op.h:414
r
ROOT::R::TRInterface & r
Definition: Object.C:4
v
SVector< double, 2 > v
Definition: Dict.h:5
TFile::Write
virtual Int_t Write(const char *name=0, Int_t opt=0, Int_t bufsiz=0)
Write memory objects to this file.
Definition: TFile.cxx:2248
ROOT::Vc::double_v
VECTOR_NAMESPACE::double_v double_v
Definition: vector.h:80
r1
unsigned int r1[N_CITIES]
Definition: simanTSP.cxx:321
testMT_S
void testMT_S(const A &a, double &time, C &result)
Definition: matrix_op.h:243
testMM_C
void testMM_C(const A &a, const B &b, const C &c, double &time, C &result)
Definition: matrix_op.h:510
TEST_TYPE
int TEST_TYPE
Definition: testOperations.cxx:731
NLOOP
int NLOOP
Definition: testOperations.cxx:46
testInnerProd_S
double testInnerProd_S(const M &a, const V &v, double &time)
Definition: matrix_op.h:180
gMatrixCheck
R__EXTERN Int_t gMatrixCheck
Definition: TMatrixTBase.h:91
TH1D
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:613
NITER
#define NITER
Definition: testOperations.cxx:27
testDot_T
double testDot_T(const V &v1, const V &v2, double &time)
Definition: matrix_op.h:324
testMeq_C
void testMeq_C(const M &m, double &time, M &result)
Definition: matrix_op.h:472
matrix_util.h
testATBA_C
void testATBA_C(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:537
fillRandomMat
void fillRandomMat(TRandom &r, M &m, unsigned int first, unsigned int second, unsigned int start=0, double offset=1)
Definition: matrix_util.h:13
test_smatrix_sym_op
int test_smatrix_sym_op()
Definition: testOperations.cxx:186
testInv_C
void testInv_C(const M &a, double &time, M &result)
Definition: matrix_op.h:523
f
double f(double x)
Definition: testIntegration.cxx:12
test_smatrix_op
int test_smatrix_op()
Definition: testOperations.cxx:61
testMT_C
void testMT_C(const A &a, double &time, C &result)
Definition: matrix_op.h:562
Stype
double Stype
Definition: testOperations.cxx:66
testInvChol_S
void testInvChol_S(const M *m, double &time, M *result)
Definition: matrix_op_vec.h:258
testGMV_T
void testGMV_T(const M &mat, const V &v1, const V &v2, double &time, V &result)
Definition: matrix_op.h:274
testATBA_S
void testATBA_S(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:210
testGMV
void testGMV(const M &mat, const V &v1, const V &v2, double &time, V &result)
Definition: matrix_op.h:124
testVad
void testVad(const V &v1, const V &v2, double &time, V &result)
Definition: matrix_op.h:59
testInvFast_S
void testInvFast_S(const M *m, double &time, M *result)
Definition: matrix_op_vec.h:239
testVscale_T
void testVscale_T(const V &v1, double a, double &time, V &result)
Definition: matrix_op.h:387
NTYPES
#define NTYPES
Definition: testOperations.cxx:743
matrix_op_vec.h
N
#define N
testInnerProd_T
double testInnerProd_T(const M &a, const V &v, double &time)
Definition: matrix_op.h:337
NLISTSIZE
#define NLISTSIZE
Definition: testOperations.cxx:59
testMad_C
void testMad_C(const M &m1, const M &m2, double &time, M &result)
Definition: matrix_op.h:484
testMT_T
void testMT_T(const A &a, double &time, C &result)
Definition: matrix_op.h:427
testMad_T
void testMad_T(const M &m1, const M &m2, double &time, M &result)
Definition: matrix_op.h:300
MATRIX_SIZE
int MATRIX_SIZE
Definition: testOperations.cxx:732
testDot_C
double testDot_C(const V &v1, const V &v2, double &time)
Definition: matrix_op.h:443
testMV
void testMV(const M &mat, const V &v, double &time, V &result)
Definition: matrix_op.h:112
testVad_T
void testVad_T(const V &v1, const V &v2, double &time, V &result)
Definition: matrix_op.h:375
testMscale_C
void testMscale_C(const M &m1, double a, double &time, M &result)
Definition: matrix_op.h:497
testATBA_S2
void testATBA_S2(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:227
r2
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322
TEST
#define TEST(N)
Definition: testOperations.cxx:744
TFile::Close
virtual void Close(Option_t *option="")
Close a file.
Definition: TFile.cxx:898
ROOT::Math::SVector
SVector: a generic fixed size Vector class.
Definition: BinaryOperators.h:31
TStopwatch
Stopwatch class.
Definition: TStopwatch.h:30
main
int main(int argc, char *argv[])
Definition: testOperations.cxx:773