ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
matrix_op.h
Go to the documentation of this file.
1 #ifndef MATRIX_OP_H
2 #define MATRIX_OP_H
3 
4 #include "TestTimer.h"
5 
6 // define functions for matrix operations
7 
8 //#define DEBUG
9 //#ifndef NLOOP
10 //#define NLOOP 1000000
11 //#endif
12 #include <vector>
13 
14 using namespace ROOT::Math;
15 
16 std::vector<float> gV;
17 
18 void initValues() {
19  gV.reserve(10*NLOOP);
20  TRandom3 r;
21  std::cout << "init smearing vector ";
22  for (int l = 0; l < 10*NLOOP; l++)
23  {
24  gV.push_back( r.Rndm() );
25  }
26  std::cout << " with size " << gV.size() << std::endl;
27 
28 }
29 
30 
31 // vector assignment
32 template<class V>
33 void testVeq(const V & v, double & time, V & result) {
34  V vtmp = v;
35  test::Timer t(time,"V=V ");
36  for (int l = 0; l < 10*NLOOP; l++)
37  {
38  vtmp[0] = gV[l];
39  result = vtmp;
40  }
41 }
42 
43 // matrix assignmnent
44 template<class M>
45 void testMeq(const M & m, double & time, M & result) {
46  M mtmp = m;
47  test::Timer t(time,"M=M ");
48  for (int l = 0; l < NLOOP; l++)
49  {
50  mtmp(0,0) = gV[l];
51  result = mtmp;
52  }
53 }
54 
55 
56 
57 // vector sum
58 template<class V>
59 void testVad(const V & v1, const V & v2, double & time, V & result) {
60  V vtmp = v2;
61  test::Timer t(time,"V+V ");
62  for (int l = 0; l < 10*NLOOP; l++)
63  {
64  vtmp[0] = gV[l];
65  result = v1 + vtmp;
66  }
67 }
68 
69 // matrix sum
70 template<class M>
71 void testMad(const M & m1, const M & m2, double & time, M & result) {
72  M mtmp = m2;
73  test::Timer t(time,"M+M ");;
74  for (int l = 0; l < NLOOP; l++)
75  {
76  mtmp(0,0) = gV[l];
77  result = m1; result += mtmp;
78  //M tmp = m1 + mtmp;
79  //result = tmp;
80  }
81 }
82 
83 // vector * constant
84 template<class V>
85 void testVscale(const V & v1, double a, double & time, V & result) {
86  V vtmp = v1;
87  test::Timer t(time,"a*V ");;
88  for (int l = 0; l < NLOOP; l++)
89  {
90  vtmp[0] = gV[l];
91  result = a * vtmp; // v1 * a does not exist in ROOT
92  }
93 }
94 
95 
96 // matrix * constant
97 template<class M>
98 void testMscale(const M & m1, double a, double & time, M & result) {
99  M mtmp = m1;
100  test::Timer t(time,"a*M ");;
101  for (int l = 0; l < NLOOP; l++)
102  {
103  mtmp(0,0) = gV[l];
104  //result = mtmp * a;
105  result = mtmp; result *= a;
106  }
107 }
108 
109 
110 // simple Matrix vector op
111 template<class M, class V>
112 void testMV(const M & mat, const V & v, double & time, V & result) {
113  V vtmp = v;
114  test::Timer t(time,"M*V ");
115  for (int l = 0; l < NLOOP; l++)
116  {
117  vtmp[0] = gV[l];
118  result = mat * vtmp;
119  }
120 }
121 
122 // general matrix vector op
123 template<class M, class V>
124 void testGMV(const M & mat, const V & v1, const V & v2, double & time, V & result) {
125  V vtmp = v1;
126  test::Timer t(time,"M*V+");
127  for (int l = 0; l < NLOOP; l++)
128  {
129  vtmp[0] = gV[l];
130  result = mat * vtmp + v2;
131  }
132 }
133 
134 
135 // general matrix matrix op
136 template<class A, class B, class C>
137 void testMM(const A & a, const B & b, const C & c, double & time, C & result) {
138  B btmp = b;
139  test::Timer t(time,"M*M ");
140  for (int l = 0; l < NLOOP; l++)
141  {
142  btmp(0,0) = gV[l];
143  result = a * btmp + c;
144  }
145 }
146 
147 
148 
149 
150 // specialized functions (depending on the package)
151 
152 //smatrix
153 template<class V>
154 double testDot_S(const V & v1, const V & v2, double & time) {
155  V vtmp = v2;
156  test::Timer t(time,"dot ");
157  double result=0;
158  for (int l = 0; l < 10*NLOOP; l++)
159  {
160  vtmp[0] = gV[l];
161  result = Dot(v1,vtmp);
162  }
163  return result;
164 }
165 
166 
167 // double testDot_S(const std::vector<V*> & w1, const std::vector<V*> & w2, double & time) {
168 // test::Timer t(time,"dot ");
169 // double result=0;
170 // for (int l = 0; l < NLOOP; l++)
171 // {
172 // V & v1 = *w1[l];
173 // V & v2 = *w2[l];
174 // result = Dot(v1,v2);
175 // }
176 // return result;
177 // }
178 
179 template<class M, class V>
180 double testInnerProd_S(const M & a, const V & v, double & time) {
181  V vtmp = v;
182  test::Timer t(time,"prod");
183  double result=0;
184  for (int l = 0; l < NLOOP; l++)
185  {
186  vtmp[0] = gV[l];
187  result = Similarity(vtmp,a);
188  }
189  return result;
190 }
191 
192 //inversion
193 template<class M>
194 void testInv_S( const M & m, double & time, M& result){
195  M mtmp = m;
196  test::Timer t(time,"inv ");
197  int ifail = 0;
198  for (int l = 0; l < NLOOP; l++)
199  {
200  mtmp(0,0) = gV[l];
201  //result = mtmp.InverseFast(ifail);
202  result = mtmp.Inverse(ifail);
203  // assert(ifail == 0);
204  }
205 }
206 
207 
208 // general matrix matrix op
209 template<class A, class B, class C>
210 void testATBA_S(const A & a, const B & b, double & time, C & result) {
211  B btmp = b;
212  test::Timer t(time,"At*M*A");
213  for (int l = 0; l < NLOOP; l++)
214  {
215  //result = Transpose(a) * b * a;
216  //result = a * b * Transpose(a);
217  //result = a * b * a;
218  btmp(0,0) = gV[l];
219  // A at = Transpose(a);
220  C tmp = btmp * Transpose(a);
221  result = a * tmp;
222  }
223 }
224 
225 // general matrix matrix op
226 template<class A, class B, class C>
227 void testATBA_S2(const A & a, const B & b, double & time, C & result) {
228  B btmp = b;
229 
230  test::Timer t(time,"At*M*A");
231  for (int l = 0; l < NLOOP; l++)
232  {
233  //result = Transpose(a) * b * a;
234  //result = a * b * Transpose(a);
235  //result = a * b * a;
236  btmp(0,0) = gV[l];
237  result = SimilarityT(a,b);
238  //result = a * result;
239  }
240 }
241 
242 template<class A, class C>
243 void testMT_S(const A & a, double & time, C & result) {
244  A atmp = a;
245  test::Timer t(time,"Transp");
246  for (int l = 0; l < NLOOP; l++)
247  {
248  //result = Transpose(a) * b * a;
249  //result = a * b * Transpose(a);
250  //result = a * b * a;
251  atmp(0,0) = gV[l];
252  result = Transpose(atmp);
253  }
254 }
255 
256 /////////////////////////////////////
257 // for root
258 //////////////////////////////////
259 
260 // simple Matrix vector op
261 template<class M, class V>
262 void testMV_T(const M & mat, const V & v, double & time, V & result) {
263  V vtmp = v;
264  test::Timer t(time,"M*V ");
265  for (int l = 0; l < NLOOP; l++)
266  {
267  vtmp[0] = gV[l];
268  Add(result,0.0,mat,vtmp);
269  }
270 }
271 
272 // general matrix vector op
273 template<class M, class V>
274 void testGMV_T(const M & mat, const V & v1, const V & v2, double & time, V & result) {
275  V vtmp = v1;
276  test::Timer t(time,"M*V+");
277  for (int l = 0; l < NLOOP; l++)
278  {
279  vtmp[0] = gV[l];
280  memcpy(result.GetMatrixArray(),v2.GetMatrixArray(),v2.GetNoElements()*sizeof(Double_t));
281  Add(result,1.0,mat,vtmp);
282  }
283 }
284 
285 // general matrix matrix op
286 template<class A, class B, class C>
287 void testMM_T(const A & a, const B & b, const C & c, double & time, C & result) {
288  B btmp = b;
289  test::Timer t(time,"M*M ");
290  for (int l = 0; l < NLOOP; l++)
291  {
292  btmp(0,0) = gV[l];
293  result.Mult(a,btmp);
294  result += c;
295  }
296 }
297 
298 // matrix sum
299 template<class M>
300 void testMad_T(const M & m1, const M & m2, double & time, M & result) {
301  M mtmp = m2;
302  test::Timer t(time,"M+M ");
303  for (int l = 0; l < NLOOP; l++)
304  {
305  mtmp(0,0) = gV[l];
306  result.Plus(m1,mtmp);
307  }
308 }
309 
310 template<class A, class B, class C>
311 void testATBA_T(const A & a, const B & b, double & time, C & result) {
312  B btmp = b;
313  test::Timer t(time,"At*M*A");
314  C tmp = a;
315  for (int l = 0; l < NLOOP; l++)
316  {
317  btmp(0,0) = gV[l];
318  tmp.Mult(a,btmp);
319  result.MultT(tmp,a);
320  }
321 }
322 
323 template<class V>
324 double testDot_T(const V & v1, const V & v2, double & time) {
325  V vtmp = v2;
326  test::Timer t(time,"dot ");
327  double result=0;
328  for (int l = 0; l < 10*NLOOP; l++)
329  {
330  vtmp[0] = gV[l];
331  result = Dot(v1,vtmp);
332  }
333  return result;
334 }
335 
336 template<class M, class V>
337 double testInnerProd_T(const M & a, const V & v, double & time) {
338  V vtmp = v;
339  test::Timer t(time,"prod");
340  double result=0;
341  for (int l = 0; l < NLOOP; l++) {
342  vtmp[0] = gV[l];
343  result = a.Similarity(vtmp);
344  }
345  return result;
346 }
347 
348 //inversion
349 template<class M>
350 void testInv_T(const M & m, double & time, M& result){
351  M mtmp = m;
352  test::Timer t(time,"inv ");
353  for (int l = 0; l < NLOOP; l++)
354  {
355  mtmp(0,0) = gV[l];
356  memcpy(result.GetMatrixArray(),mtmp.GetMatrixArray(),mtmp.GetNoElements()*sizeof(Double_t));
357  result.InvertFast();
358  }
359 }
360 
361 template<class M>
362 void testInv_T2(const M & m, double & time, M& result){
363  M mtmp = m;
364  test::Timer t(time,"inv2");
365  for (int l = 0; l < NLOOP; l++)
366  {
367  memcpy(result.GetMatrixArray(),mtmp.GetMatrixArray(),mtmp.GetNoElements()*sizeof(Double_t));
368  result.InvertFast();
369  }
370 }
371 
372 
373 // vector sum
374 template<class V>
375 void testVad_T(const V & v1, const V & v2, double & time, V & result) {
376  V vtmp = v2;
377  test::Timer t(time,"V+V ");;
378  for (int l = 0; l < 10*NLOOP; l++)
379  {
380  vtmp[0] = gV[l];
381  result.Add(v1,vtmp);
382  }
383 }
384 
385 // vector * constant
386 template<class V>
387 void testVscale_T(const V & v1, double a, double & time, V & result) {
388  V vtmp = v1;
389  test::Timer t(time,"a*V ");;
390  for (int l = 0; l < NLOOP; l++)
391  {
392  // result = a * v1;
393  result.Zero();
394  vtmp[0] = gV[l];
395  Add(result,a,vtmp);
396  }
397 }
398 
399 // general matrix matrix op
400 template<class A, class B, class C>
401 void testATBA_T2(const A & a, const B & b, double & time, C & result) {
402  B btmp = b;
403  test::Timer t(time,"At*M*A");
404  for (int l = 0; l < NLOOP; l++)
405  {
406  btmp(0,0) = gV[l];
407  memcpy(result.GetMatrixArray(),btmp.GetMatrixArray(),btmp.GetNoElements()*sizeof(Double_t));
408  result.Similarity(a);
409  }
410 }
411 
412 // matrix * constant
413 template<class M>
414 void testMscale_T(const M & m1, double a, double & time, M & result) {
415  M mtmp = m1;
416  test::Timer t(time,"a*M ");;
417  for (int l = 0; l < NLOOP; l++)
418  {
419  //result = a * m1;
420  result.Zero();
421  mtmp(0,0) = gV[l];
422  Add(result,a,mtmp);
423  }
424 }
425 
426 template<class A, class C>
427 void testMT_T(const A & a, double & time, C & result) {
428  A atmp = a;
429  test::Timer t(time,"Transp");
430  for (int l = 0; l < NLOOP; l++)
431  {
432  atmp(0,0) = gV[l];
433  result.Transpose(atmp);
434  }
435 }
436 
437 ////////////////////////////////////////////
438 // for clhep
439 ////////////////////////////////////////////
440 
441 //smatrix
442 template<class V>
443 double testDot_C(const V & v1, const V & v2, double & time) {
444  V vtmp = v2;
445  test::Timer t(time,"dot ");
446  double result=0;
447  for (int l = 0; l < 10*NLOOP; l++)
448  {
449  vtmp[0] = gV[l];
450  result = dot(v1,vtmp);
451  }
452  return result;
453 }
454 
455 template<class M, class V>
456 double testInnerProd_C(const M & a, const V & v, double & time) {
457  V vtmp = v;
458  test::Timer t(time,"prod");
459  double result=0;
460  for (int l = 0; l < NLOOP; l++)
461  {
462  vtmp[0] = gV[l];
463  V tmp = a*vtmp;
464  result = dot(vtmp,tmp);
465  }
466  return result;
467 }
468 
469 
470 // matrix assignmnent(index starts from 1)
471 template<class M>
472 void testMeq_C(const M & m, double & time, M & result) {
473  M mtmp = m;
474  test::Timer t(time,"M=M ");
475  for (int l = 0; l < NLOOP; l++)
476  {
477  mtmp(1,1) = gV[l];
478  result = mtmp;
479  }
480 }
481 
482 // matrix sum
483 template<class M>
484 void testMad_C(const M & m1, const M & m2, double & time, M & result) {
485  M mtmp = m2;
486  test::Timer t(time,"M+M ");;
487  for (int l = 0; l < NLOOP; l++)
488  {
489  mtmp(1,1) = gV[l];
490  result = m1; result += mtmp;
491  }
492 }
493 
494 
495 // matrix * constant
496 template<class M>
497 void testMscale_C(const M & m1, double a, double & time, M & result) {
498  M mtmp = m1;
499  test::Timer t(time,"a*M ");;
500  for (int l = 0; l < NLOOP; l++)
501  {
502  mtmp(1,1) = gV[l];
503  result = mtmp * a;
504  }
505 }
506 
507 
508 // clhep matrix matrix op (index starts from 1)
509 template<class A, class B, class C>
510 void testMM_C(const A & a, const B & b, const C & c, double & time, C & result) {
511  B btmp = b;
512  test::Timer t(time,"M*M ");
513  for (int l = 0; l < NLOOP; l++)
514  {
515  btmp(1,1) = gV[l];
516  result = a * btmp + c;
517  }
518 }
519 
520 
521 //inversion
522 template<class M>
523 void testInv_C( const M & a, double & time, M& result){
524  M mtmp = a;
525  test::Timer t(time,"inv ");
526  int ifail = 0;
527  for (int l = 0; l < NLOOP; l++)
528  {
529  mtmp(1,1) = gV[l];
530  result = mtmp.inverse(ifail);
531  if (ifail) {std::cout <<"error inverting" << mtmp << std::endl; return; }
532  }
533 }
534 
535 // general matrix matrix op
536 template<class A, class B, class C>
537 void testATBA_C(const A & a, const B & b, double & time, C & result) {
538  B btmp = b;
539  test::Timer t(time,"At*M*A");
540  for (int l = 0; l < NLOOP; l++)
541  {
542  btmp(1,1) = gV[l];
543  //result = a.T() * b * a;
544  result = a * btmp * a.T();
545  }
546 }
547 
548 
549 template<class A, class B, class C>
550 void testATBA_C2(const A & a, const B & b, double & time, C & result) {
551  B btmp = b;
552  test::Timer t(time,"At*M*A");
553  for (int l = 0; l < NLOOP; l++)
554  {
555  btmp(1,1) = gV[l];
556  result = btmp.similarity(a);
557  }
558 }
559 
560 
561 template<class A, class C>
562 void testMT_C(const A & a, double & time, C & result) {
563  A atmp = a;
564  test::Timer t(time,"Transp");
565  for (int l = 0; l < NLOOP; l++)
566  {
567  atmp(1,1) = gV[l];
568  result = atmp.T();
569  }
570 }
571 
572 
573 #endif
T Dot(const SVector< T, D > &lhs, const SVector< T, D > &rhs)
Vector dot product.
Definition: Functions.h:166
void testATBA_C2(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:550
T Similarity(const SMatrix< T, D, D, R > &lhs, const SVector< T, D > &rhs)
Similarity Vector - Matrix Product: v^T * A * v returning a scalar value of type T ...
void testVeq(const V &v, double &time, V &result)
Definition: matrix_op.h:33
static double B[]
void testInv_T2(const M &m, double &time, M &result)
Definition: matrix_op.h:362
void testATBA_T(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:311
Random number generator class based on M.
Definition: TRandom3.h:29
const Double_t * v1
Definition: TArcBall.cxx:33
void testMM(const A &a, const B &b, const C &c, double &time, C &result)
Definition: matrix_op.h:137
TMatrixT< Element > & Add(TMatrixT< Element > &target, Element scalar, const TMatrixT< Element > &source)
Modify addition: target += scalar * source.
Definition: TMatrixT.cxx:2925
return c
void testMV_T(const M &mat, const V &v, double &time, V &result)
Definition: matrix_op.h:262
void testVscale(const V &v1, double a, double &time, V &result)
Definition: matrix_op.h:85
void testATBA_T2(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:401
Expr< TransposeOp< SMatrix< T, D1, D2, R >, T, D1, D2 >, T, D2, D1, typename TranspPolicy< T, D1, D2, R >::RepType > Transpose(const SMatrix< T, D1, D2, R > &rhs)
Matrix Transpose B(i,j) = A(j,i) returning a matrix expression.
TArc * a
Definition: textangle.C:12
void testMscale(const M &m1, double a, double &time, M &result)
Definition: matrix_op.h:98
static double A[]
double testDot_S(const V &v1, const V &v2, double &time)
Definition: matrix_op.h:154
Double_t dot(const TVector2 &v1, const TVector2 &v2)
Definition: CsgOps.cxx:333
#define NLOOP
Definition: piRandom.C:21
void testMM_T(const A &a, const B &b, const C &c, double &time, C &result)
Definition: matrix_op.h:287
void testInv_T(const M &m, double &time, M &result)
Definition: matrix_op.h:350
virtual Double_t Rndm(Int_t i=0)
Machine independent random number generator.
Definition: TRandom3.cxx:94
void testMeq(const M &m, double &time, M &result)
Definition: matrix_op.h:45
void initValues()
Definition: matrix_op.h:18
void testMad(const M &m1, const M &m2, double &time, M &result)
Definition: matrix_op.h:71
double testInnerProd_C(const M &a, const V &v, double &time)
Definition: matrix_op.h:456
void testInv_S(const M &m, double &time, M &result)
Definition: matrix_op.h:194
TThread * t[5]
Definition: threadsh1.C:13
RooCmdArg Timer(Bool_t flag=kTRUE)
static double C[]
void testMscale_T(const M &m1, double a, double &time, M &result)
Definition: matrix_op.h:414
ROOT::R::TRInterface & r
Definition: Object.C:4
SVector< double, 2 > v
Definition: Dict.h:5
void testMT_S(const A &a, double &time, C &result)
Definition: matrix_op.h:243
void testMM_C(const A &a, const B &b, const C &c, double &time, C &result)
Definition: matrix_op.h:510
TMarker * m
Definition: textangle.C:8
TLine * l
Definition: textangle.C:4
double testInnerProd_S(const M &a, const V &v, double &time)
Definition: matrix_op.h:180
double testDot_T(const V &v1, const V &v2, double &time)
Definition: matrix_op.h:324
void testMeq_C(const M &m, double &time, M &result)
Definition: matrix_op.h:472
void testATBA_C(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:537
void testInv_C(const M &a, double &time, M &result)
Definition: matrix_op.h:523
double Double_t
Definition: RtypesCore.h:55
void testMT_C(const A &a, double &time, C &result)
Definition: matrix_op.h:562
void testGMV_T(const M &mat, const V &v1, const V &v2, double &time, V &result)
Definition: matrix_op.h:274
void testATBA_S(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:210
std::vector< float > gV
Definition: matrix_op.h:16
void testGMV(const M &mat, const V &v1, const V &v2, double &time, V &result)
Definition: matrix_op.h:124
void testVad(const V &v1, const V &v2, double &time, V &result)
Definition: matrix_op.h:59
void testVscale_T(const V &v1, double a, double &time, V &result)
Definition: matrix_op.h:387
double testInnerProd_T(const M &a, const V &v, double &time)
Definition: matrix_op.h:337
SMatrix< T, D2, D2, MatRepSym< T, D2 > > SimilarityT(const SMatrix< T, D1, D2, R > &lhs, const SMatrix< T, D1, D1, MatRepSym< T, D1 > > &rhs)
Transpose Similarity Matrix Product : B = U^T * A * U for A symmetric returning a symmetric matrix ex...
void testMad_C(const M &m1, const M &m2, double &time, M &result)
Definition: matrix_op.h:484
void testMT_T(const A &a, double &time, C &result)
Definition: matrix_op.h:427
void testMad_T(const M &m1, const M &m2, double &time, M &result)
Definition: matrix_op.h:300
double testDot_C(const V &v1, const V &v2, double &time)
Definition: matrix_op.h:443
double result[121]
void testMV(const M &mat, const V &v, double &time, V &result)
Definition: matrix_op.h:112
void testVad_T(const V &v1, const V &v2, double &time, V &result)
Definition: matrix_op.h:375
void testMscale_C(const M &m1, double a, double &time, M &result)
Definition: matrix_op.h:497
void testATBA_S2(const A &a, const B &b, double &time, C &result)
Definition: matrix_op.h:227