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