ROOT  6.06/09
Reference Guide
testKalman.cxx
Go to the documentation of this file.
1 #ifdef USE_VC
2 //using namespace Vc;
3 #include "Vc/Vc"
4 #endif
5 
6 
7 
8 #include <cassert>
9 
10 
11 
12 #include "Math/SVector.h"
13 #include "Math/SMatrix.h"
14 
15 #include "TMatrixD.h"
16 #include "TVectorD.h"
17 
18 #include "TRandom3.h"
19 
20 
21 
22 #include "matrix_util.h"
23 
24 
25 #define TEST_SYM
26 
27 //#define HAVE_CLHEP
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 // #include "SealUtil/SealTimer.h"
35 // #include "SealUtil/SealHRRTChrono.h"
36 // #include "SealUtil/TimingReport.h"
37 
38 #include <iostream>
39 
40 
41 #include <sys/times.h>
42 #include <unistd.h>
43 
44 double cpuTime()
45 {
46  struct tms usage;
47  times(&usage);
48  return ((double) usage.tms_utime) / sysconf(_SC_CLK_TCK);
49 }
50 
51 double clockTime()
52 {
53  struct tms usage;
54  return ((double) times(&usage)) / sysconf(_SC_CLK_TCK);
55 }
56 
57 
58 #ifndef NDIM1
59 #define NDIM1 5
60 #endif
61 #ifndef NDIM2
62 #define NDIM2 5
63 #endif
64 
65 #define NITER 1 // number of iterations
66 
67 #define NLOOP 500000 // number of time the test is repeted
68 #define NLISTSIZE 64 // size of matrix/vector lists
69 
70 using namespace ROOT::Math;
71 
72 #include "TestTimer.h"
73 
74 #ifdef USE_VC
75 typedef Vc::double_v Stype;
76 #else
77 typedef double Stype;
78 #endif
79 
80 #ifdef USE_VC
81 const int NLIST = NLISTSIZE / Vc::double_v::Size;
82 #else
83 const int NLIST = NLISTSIZE;
84 #endif
85 
86 
88 
89  // need to write explicitly the dimensions
90 
91 
92  typedef SMatrix<Stype, NDIM1, NDIM1> MnMatrixNN;
93  typedef SMatrix<Stype, NDIM2, NDIM2> MnMatrixMM;
94  typedef SMatrix<Stype, NDIM1, NDIM2> MnMatrixNM;
95  typedef SMatrix<Stype, NDIM2 , NDIM1> MnMatrixMN;
96  typedef SMatrix<Stype, NDIM1 > MnSymMatrixNN;
97  typedef SMatrix<Stype, NDIM2 > MnSymMatrixMM;
98  typedef SVector<Stype, NDIM1> MnVectorN;
99  typedef SVector<Stype, NDIM2> MnVectorM;
100 
101 
102 
103  int first = NDIM1; //Can change the size of the matrices
104  int second = NDIM2;
105 
106 
107  std::cout << "************************************************\n";
108  std::cout << " SMatrix kalman test " << first << " x " << second << std::endl;
109  std::cout << "************************************************\n";
110 
111 
112 
113 
114  int npass = NITER;
115  TRandom3 r(111);
116  Stype x2sum = 0.0;
117  Stype c2sum = 0.0;
118 
119  for (int ipass = 0; ipass < npass; ipass++) {
120 
121 
122  MnMatrixNM H[NLIST];
123  MnMatrixMN K0[NLIST];
124  MnSymMatrixMM Cp[NLIST];
125  MnSymMatrixNN V[NLIST];
126  MnVectorN m[NLIST];
127  MnVectorM xp[NLIST];
128 
129 
130  // fill matrices with random data
131  for (int j = 0; j < NLIST; j++) fillRandomMat(r,H[j],first,second);
132  for (int j = 0; j < NLIST; j++) fillRandomMat(r,K0[j],second,first);
133  for (int j = 0; j < NLIST; j++) fillRandomSym(r,Cp[j],second);
134  for (int j = 0; j < NLIST; j++) fillRandomSym(r,V[j],first);
135  for (int j = 0; j < NLIST; j++) fillRandomVec(r,m[j],first);
136  for (int j = 0; j < NLIST; j++) fillRandomVec(r,xp[j],second);
137 
138 
139  // MnSymMatrixMM I;
140  // for(int i = 0; i < second; i++)
141  // I(i,i) = 1;
142 
143 #ifdef DEBUG
144  std::cout << "pass " << ipass << std::endl;
145  if (k == 0) {
146  std::cout << " K0 = " << K0[0] << std::endl;
147  std::cout << " H = " << K0[0] << std::endl;
148  std::cout << " Cp = " << Cp[0] << std::endl;
149  std::cout << " V = " << V[0] << std::endl;
150  std::cout << " m = " << m[0] << std::endl;
151  std::cout << " xp = " << xp[0] << std::endl;
152  }
153 #endif
154 
155 
156  {
157  Stype x2 = 0.0,c2 = 0.0;
158  test::Timer t("SMatrix Kalman ");
159 
160  MnVectorM x;
161  MnMatrixMN tmp;
162  MnSymMatrixNN Rinv;
163  MnMatrixMN K;
164  MnSymMatrixMM C;
165  MnVectorN vtmp1;
166  MnVectorN vtmp;
167 
168  for (int l = 0; l < NLOOP; l++) {
169 
170  // loop on the list of matrices
171  for (int k = 0; k < NLIST; k++) {
172 
173 
174 
175  vtmp1 = H[k]*xp[k] -m[k];
176  //x = xp + K0 * (m- H * xp);
177  x = xp[k] - K0[k] * vtmp1;
178  tmp = Cp[k] * Transpose(H[k]);
179  Rinv = V[k]; Rinv += H[k] * tmp;
180 
181  bool test = Rinv.InvertFast();
182  if(!test) {
183  std::cout<<"inversion failed" <<std::endl;
184  std::cout << Rinv << std::endl;
185  }
186 
187  K = tmp * Rinv ;
188  C = Cp[k]; C -= K * Transpose(tmp);
189  //C = ( I - K * H ) * Cp;
190  //x2 = Product(Rinv,m-H*xp); // this does not compile on WIN32
191  vtmp = m[k]-H[k]*xp[k];
192  x2 = Dot(vtmp, Rinv*vtmp);
193 
194 
195  //std::cout << k << " chi2 = " << x2 << std::endl;
196  x2sum += x2;
197  c2 = 0;
198  for (int i=0; i<NDIM2; ++i)
199  for (int j=0; j<NDIM2; ++j)
200  c2 += C(i,j);
201 
202  c2sum += c2;
203  }
204  }
205 
206  //tr.dump();
207  }
208 
209  std::cerr << "SMatrix: x2sum = " << x2sum << "\tc2sum = " << c2sum << std::endl;
210 #ifdef USE_VC
211  double sx2=0; double sc2=0;
212  for (int l=0;l<Stype::Size; ++l) {
213  sx2 += x2sum[l];
214  sc2 += c2sum[l];
215  }
216  std::cerr << "SMatrix: x2sum = " << sx2 << "\tc2sum = " << sc2 << std::endl;
217 #endif
218  }
219  return 0;
220 
221 }
222 
224 
225 #ifdef TEST_SYM
226 
227  // need to write explicitly the dimensions
228 
229 
230  typedef SMatrix<Stype, NDIM1, NDIM1> MnMatrixNN;
231  typedef SMatrix<Stype, NDIM2, NDIM2> MnMatrixMM;
232  typedef SMatrix<Stype, NDIM1, NDIM2> MnMatrixNM;
233  typedef SMatrix<Stype, NDIM2 , NDIM1> MnMatrixMN;
236  typedef SVector<Stype, NDIM1> MnVectorN;
237  typedef SVector<Stype, NDIM2> MnVectorM;
238  typedef SVector<Stype, NDIM1*(NDIM1+1)/2> MnVectorN2;
239  typedef SVector<Stype, NDIM2*(NDIM2+1)/2> MnVectorM2;
240 
241 
242 
243  int first = NDIM1; //Can change the size of the matrices
244  int second = NDIM2;
245 
246 
247  std::cout << "************************************************\n";
248  std::cout << " SMatrix_SyM kalman test " << first << " x " << second << std::endl;
249  std::cout << "************************************************\n";
250 
251 
252 
253  int npass = NITER;
254  TRandom3 r(111);
255  Stype x2sum = 0.0, c2sum = 0.0;
256 
257  for (int ipass = 0; ipass < npass; ipass++) {
258 
259 
260  MnMatrixNM H[NLIST];
261  MnMatrixMN K0[NLIST];
262  MnSymMatrixMM Cp[NLIST];
263  MnSymMatrixNN V[NLIST];
264  MnVectorN m[NLIST];
265  MnVectorM xp[NLIST];
266 
267 
268  // fill matrices with random data
269  for (int j = 0; j < NLIST; j++) fillRandomMat(r,H[j],first,second);
270  for (int j = 0; j < NLIST; j++) fillRandomMat(r,K0[j],second,first);
271  for (int j = 0; j < NLIST; j++) fillRandomSym(r,Cp[j],second);
272  for (int j = 0; j < NLIST; j++) fillRandomSym(r,V[j],first);
273  for (int j = 0; j < NLIST; j++) fillRandomVec(r,m[j],first);
274  for (int j = 0; j < NLIST; j++) fillRandomVec(r,xp[j],second);
275 
276 
277 
278  // MnSymMatrixMM I;
279  // for(int i = 0; i < second; i++)
280  // I(i,i) = 1;
281 
282 #ifdef DEBUG
283  std::cout << "pass " << ipass << std::endl;
284  if (ipass == 0) {
285  std::cout << " K0 = " << K0 << std::endl;
286  std::cout << " H = " << K0 << std::endl;
287  std::cout << " Cp = " << Cp << std::endl;
288  std::cout << " V = " << V << std::endl;
289  std::cout << " m = " << m << std::endl;
290  std::cout << " xp = " << xp << std::endl;
291  }
292 #endif
293 
294 
295  //double clc1 = clockTime();
296 
297  Stype x2 = 0.0,c2 = 0.0;
298  test::Timer t("SMatrix Kalman ");
299 
300  MnVectorM x;
301  MnSymMatrixNN RinvSym;
302  MnMatrixMN K;
303  MnSymMatrixMM C;
304  MnSymMatrixMM Ctmp;
305  MnVectorN vtmp1;
306  MnVectorN vtmp;
307 #define OPTIMIZED_SMATRIX_SYM
308 #ifdef OPTIMIZED_SMATRIX_SYM
309  MnMatrixMN tmp;
310 #endif
311 
312  for (int l = 0; l < NLOOP; l++) {
313  for (int k = 0; k < NLIST; k++) {
314 
315 
316 #ifdef OPTIMIZED_SMATRIX_SYM
317  vtmp1 = H[k]*xp[k] -m[k];
318  //x = xp + K0 * (m- H * xp);
319  x = xp[k] - K0[k] * vtmp1;
320  tmp = Cp[k] * Transpose(H[k]);
321  // we are sure that H*tmp result is symmetric
322  AssignSym::Evaluate(RinvSym,H[k]*tmp);
323  RinvSym += V[k];
324 
325  bool test = RinvSym.InvertFast();
326  if(!test) {
327  std::cout<<"inversion failed" <<std::endl;
328  std::cout << RinvSym << std::endl;
329  }
330 
331  K = tmp * RinvSym ;
332  // we profit from the fact that result of K*tmpT is symmetric
333  AssignSym::Evaluate(Ctmp, K*Transpose(tmp) );
334  C = Cp[k]; C -= Ctmp;
335  //C = ( I - K * H ) * Cp;
336  //x2 = Product(Rinv,m-H*xp); // this does not compile on WIN32
337  vtmp = m[k]-H[k]*xp[k];
338  x2 = Dot(vtmp, RinvSym*vtmp);
339 #else
340  // use similarity function
341  vtmp1 = H[k]*xp[k] -m[k];
342  x = xp[k] - K0[k] * vtmp1;
343  RinvSym = V[k]; RinvSym += Similarity(H[k],Cp[k]);
344 
345  bool test = RinvSym.InvertFast();
346  if(!test) {
347  std::cout<<"inversion failed" <<std::endl;
348  std::cout << RinvSym << std::endl;
349  }
350 
351  Ctmp = SimilarityT(H[k], RinvSym);
352  C = Cp[k]; C -= Similarity(Cp[k], Ctmp);
353  vtmp = m[k]-H[k]*xp[k];
354  x2 = Similarity(vtmp, RinvSym);
355 #endif
356 
357  x2sum += x2;
358  c2 = 0;
359  for (int i=0; i<NDIM2; ++i)
360  for (int j=0; j<NDIM2; ++j)
361  c2 += C(i,j);
362  c2sum += c2;
363 
364  } // end loop on list
365  //std::cout << k << " chi2 = " << x2 << std::endl;
366 
367  } // end loop on trials
368 
369  //tr.dump();
370 
371  // double clc2 = clockTime();
372  // std::cerr << "Time: " << (clc2 - clc1) << std::endl;
373 
374  std::cerr << "SMatrixSym: x2sum = " << x2sum << "\tc2sum = " << c2sum << std::endl;
375 #ifdef USE_VC
376  double sx2=0; double sc2=0;
377  for (int l=0;l<Stype::Size; ++l) {
378  sx2 += x2sum[l];
379  sc2 += c2sum[l];
380  }
381  std::cerr << "SMatrixSym: x2sum = " << sx2 << "\tc2sum = " << sc2 << std::endl;
382 #endif
383 
384 
385 #endif
386  }
387 
388  return 0;
389 
390 }
391 
392 
393 
394 // ROOT test
395 
396 
398 
399 #ifdef USE_TMATRIX
400 
401 
402 
403  typedef TMatrixD MnMatrix;
404  typedef TVectorD MnVector;
405 
406  // typedef boost::numeric::ublas::matrix<Stype> MnMatrix;
407  //typedef HepSymMatrix MnSymMatrixHep;
408 
409 
410  int first = NDIM1; //Can change the size of the matrices
411  int second = NDIM2;
412 
413 
414  std::cout << "************************************************\n";
415  std::cout << " TMatrix Kalman test " << first << " x " << second << std::endl;
416  std::cout << "************************************************\n";
417 
418 
419 
420  int npass = NITER;
421  TRandom3 r(111);
422  gMatrixCheck = 0;
423  Stype x2sum = 0,c2sum = 0;
424 
425  for (int k = 0; k < npass; k++)
426  {
427 
428  MnMatrix H(first,second);
429  MnMatrix K0(second,first);
430  MnMatrix Cp(second,second);
431  MnMatrix V(first,first);
432  MnVector m(first);
433  MnVector xp(second);
434 
435  // fill matrices with random data
436  fillRandomMat(r,H,first,second);
437  fillRandomMat(r,K0,second,first);
438  fillRandomSym(r,Cp,second);
439  fillRandomSym(r,V,first);
440  fillRandomVec(r,m,first);
441  fillRandomVec(r,xp,second);
442 
443 
444  MnMatrix I(second,second);//Identity matrix
445  for (int i = 0; i < second; i++)
446  for(int j = 0; j <second; j++) {
447  I(i,j) = 0.0;
448  if (i==j) I(i,i) = 1.0;
449  }
450 
451  // if (k==0) {
452  // std::cout << " Cp " << std::endl;
453  // Cp.Print();
454  // }
455 
456  {
457  Stype x2 = 0,c2 = 0;
458  TVectorD x(second);
459  TMatrixD Rinv(first,first);
460  TMatrixD Rtmp(first,first);
461  TMatrixDSym RinvSym;
462  TMatrixD K(second,first);
463  TMatrixD C(second,second);
464  TMatrixD Ctmp(second,second);
465  TVectorD tmp1(first);
466  TMatrixD tmp2(second,first);
467 #define OPTIMIZED_TMATRIX
468 #ifndef OPTIMIZED_TMATRIX
469  TMatrixD HT(second,first);
470  TMatrixD tmp2T(first,second);
471 #endif
472 
473  test::Timer t("TMatrix Kalman ");
474  for (Int_t l = 0; l < NLOOP; l++)
475  {
476 #ifdef OPTIMIZED_TMATRIX
477  tmp1 = m; Add(tmp1,-1.0,H,xp);
478  x = xp; Add(x,+1.0,K0,tmp1);
479  tmp2.MultT(Cp,H);
480  Rtmp.Mult(H,tmp2);
481  Rinv.Plus(V,Rtmp);
482  RinvSym.Use(first,Rinv.GetMatrixArray());
483  RinvSym.InvertFast();
484  K.Mult(tmp2,Rinv);
485  Ctmp.MultT(K,tmp2);
486  C.Minus(Cp,Ctmp);
487  x2 = RinvSym.Similarity(tmp1);
488 
489 #else
490  tmp1 = H*xp -m;
491  //x = xp + K0 * (m- H * xp);
492  x = xp - K0 * tmp1;
493  tmp2 = Cp * HT.Transpose(H);
494  Rinv = V; Rinv += H * tmp2;
495  RinvSym.Use(first,Rinv.GetMatrixArray());
496  RinvSym.InvertFast();
497  K= tmp2* Rinv;
498  C = Cp; C -= K*tmp2T.Transpose(tmp2);
499  x2= RinvSym.Similarity(tmp1);
500 #endif
501 
502  }
503  x2sum += x2;
504  c2 = 0;
505  for (int i=0; i<NDIM2; ++i)
506  for (int j=0; j<NDIM2; ++j)
507  c2 += C(i,j);
508  c2sum += c2;
509  }
510 
511  // }
512  }
513  //tr.dump();
514  std::cerr << "TMatrix: x2sum = " << x2sum << "\tc2sum = " << c2sum << std::endl;
515 
516 #endif
517  return 0;
518 
519 }
520 
521 
522 
523 // test CLHEP Kalman
524 
525 #ifdef HAVE_CLHEP
526 int test_clhep_kalman() {
527 
528 
529 
530  typedef HepSymMatrix MnSymMatrix;
531  typedef HepMatrix MnMatrix;
532  typedef HepVector MnVector;
533 
534 
535  // typedef boost::numeric::ublas::matrix<Stype> MnMatrix;
536  //typedef HepSymMatrix MnSymMatrixHep;
537 
538 
539  int first = NDIM1; //Can change the size of the matrices
540  int second = NDIM2;
541 
542 
543  std::cout << "************************************************\n";
544  std::cout << " CLHEP Kalman test " << first << " x " << second << std::endl;
545  std::cout << "************************************************\n";
546 
547 
548 
549  int npass = NITER;
550  TRandom3 r(111);
551  Stype x2sum = 0,c2sum = 0;
552 
553  for (int k = 0; k < npass; k++)
554  {
555 
556  // in CLHEP index starts from 1
557  MnMatrix H(first,second);
558  MnMatrix K0(second,first);
559  MnMatrix Cp(second,second);
560  MnMatrix V(first,first);
561  MnVector m(first);
562  MnVector xp(second);
563 
564  // fill matrices with random data
565  fillRandomMat(r,H,first,second,1);
566  fillRandomMat(r,K0,second,first,1);
567  fillRandomSym(r,Cp,second,1);
568  fillRandomSym(r,V,first,1);
569  fillRandomVec(r,m,first);
570  fillRandomVec(r,xp,second);
571 
572  MnSymMatrix I(second,1);//Identity matrix
573 
574  {
575  Stype x2 = 0,c2 = 0;
576  MnVector x(second);
577  MnMatrix Rinv(first,first);
578  MnSymMatrix RinvSym(first);
579  MnMatrix K(second,first);
580  MnSymMatrix C(second);
581  MnVector vtmp1(first);
582  MnMatrix tmp(second,first);
583 
584  test::Timer t("CLHEP Kalman ");
585  int ifail;
586  for (Int_t l = 0; l < NLOOP; l++)
587  {
588 
589 
590  vtmp1 = H*xp -m;
591  //x = xp + K0 * (m- H * xp);
592  x = xp - K0 * vtmp1;
593  tmp = Cp * H.T();
594  Rinv = V; Rinv += H * tmp;
595  RinvSym.assign(Rinv);
596  RinvSym.invert(ifail);
597  if (ifail !=0) { std::cout << "Error inverting Rinv" << std::endl; break; }
598  K = tmp*RinvSym;
599  //C.assign( (I-K*H)*Cp);
600  //C = (I-K*H)*Cp;
601  C.assign( (I-K*H)*Cp );
602  x2= RinvSym.similarity(vtmp1);
603  if(ifail!=0) { std::cout << "Error inverting Rinv" << std::endl; break; }
604  }
605  // std::cout << k << " chi2 " << x2 << std::endl;
606  x2sum += x2;
607 
608  c2 = 0;
609  for (int i=1; i<=NDIM2; ++i)
610  for (int j=1; j<=NDIM2; ++j)
611  c2 += C(i,j);
612  c2sum += c2;
613  }
614 
615  // }
616  }
617  //tr.dump();
618  std::cerr << "x2sum = " << x2sum << "\tc2sum = " << c2sum << std::endl;
619 
620  return 0;
621 }
622 #endif
623 
624 
625 
626 int testKalman() {
627 
628  //int nlist = NLISTSIZE;
629 #ifdef USE_VC
630  std::cout << "Using VC library - size = " << ROOT::Vc::double_v::Size << " VC_IMPL = " << VC_IMPL << std::endl;
631  //nlist /= Vc::double_v::Size;
632 #endif
633 
634  std::cout << " making vector/matrix lists of size = " << NLIST << std::endl;
635 
636 
637 #ifdef TEST_SYM
639 #endif
640 
643 #ifdef HAVE_CLHEP
644  test_clhep_kalman();
645 #endif
646 
647  return 0;
648 
649 
650 }
651 
652 int main() {
653 
654  return testKalman();
655 }
T Dot(const SVector< T, D > &lhs, const SVector< T, D > &rhs)
Vector dot product.
Definition: Functions.h:166
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 ...
const int npass
Definition: testPermute.cxx:21
Random number generator class based on M.
Definition: TRandom3.h:29
#define VC_IMPL
Definition: global.h:429
#define NITER
Definition: testKalman.cxx:65
void Add(THist< DIMENSION, PRECISIONA > &to, THist< DIMENSION, PRECISIONB > &from)
Definition: THist.h:335
TMatrixT< Element > & Transpose(const TMatrixT< Element > &source)
Transpose matrix source.
Definition: TMatrixT.cxx:1460
TMatrixTSym< Element > & Use(Int_t row_lwb, Int_t row_upb, Element *data)
double Stype
Definition: testKalman.cxx:77
const char * Size
Definition: TXMLSetup.cxx:56
#define H(x, y, z)
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.
void fillRandomVec(TRandom &r, V &v, unsigned int len, unsigned int start=0, double offset=1)
Definition: matrix_util.h:7
void fillRandomSym(TRandom &r, M &m, unsigned int first, unsigned int start=0, double offset=1)
Definition: matrix_util.h:20
int Int_t
Definition: RtypesCore.h:41
#define NDIM1
Definition: testKalman.cxx:59
void Minus(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix summation. Create a matrix C such that C = A - B.
Definition: TMatrixT.cxx:571
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
SMatrix: a generic fixed size D1 x D2 Matrix class.
void MultT(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix multiplication. Create a matrix C such that C = A * B^T.
Definition: TMatrixT.cxx:942
#define NDIM2
Definition: testKalman.cxx:62
int test_smatrix_kalman()
Definition: testKalman.cxx:44
RooCmdArg Timer(Bool_t flag=kTRUE)
static double C[]
ROOT::R::TRInterface & r
Definition: Object.C:4
VECTOR_NAMESPACE::double_v double_v
Definition: vector.h:80
int test_tmatrix_kalman()
Definition: testKalman.cxx:324
TMarker * m
Definition: textangle.C:8
Double_t K()
Definition: TMath.h:95
TLine * l
Definition: textangle.C:4
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:211
R__EXTERN Int_t gMatrixCheck
Definition: TMatrixTBase.h:91
static void Evaluate(SMatrix< T, D, D, MatRepSym< T, D > > &lhs, const Expr< A, T, D, D, R > &rhs)
assign a symmetric matrix from an expression
Definition: HelperOps.h:156
const int NLIST
Definition: testKalman.cxx:83
void fillRandomMat(TRandom &r, M &m, unsigned int first, unsigned int second, unsigned int start=0, double offset=1)
Definition: matrix_util.h:13
return c2
Definition: legend2.C:14
TMatrixTSym< Element > & InvertFast(Double_t *det=0)
Invert the matrix and calculate its determinant.
TMatrixTSym< Element > & Similarity(const TMatrixT< Element > &n)
Calculate B * (*this) * B^T , final matrix will be (nrowsb x nrowsb) This is a similarity transform w...
double clockTime()
Definition: testKalman.cxx:51
int testKalman()
Definition: testKalman.cxx:549
void Plus(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix summation. Create a matrix C such that C = A + B.
Definition: TMatrixT.cxx:503
double cpuTime()
Definition: testKalman.cxx:44
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...
int main()
Definition: testKalman.cxx:566
void Mult(const TMatrixT< Element > &a, const TMatrixT< Element > &b)
General matrix multiplication. Create a matrix C such that C = A * B.
Definition: TMatrixT.cxx:639
#define I(x, y, z)
int test_smatrix_sym_kalman()
Definition: testKalman.cxx:168
#define NLISTSIZE
Definition: testKalman.cxx:68
#define NLOOP
Definition: testKalman.cxx:67
SVector: a generic fixed size Vector class.