Math/test/Matrices_t.cpp

0001 // test of Matrices
0002
0003 #define SMATRIX_USE_CONSTEXPR
0004 #include "DataFormats/Math/interface/invertPosDefMatrix.h"
0005 #include "Math/SMatrix.h"
0006 #include "Math/SVector.h"
0007 #include <memory>
0008 #include "FWCore/Utilities/interface/HRRealTime.h"
0009 #include <iostream>
0010
0011 typedef ROOT::Math::SMatrix<double, 2, 2, ROOT::Math::MatRepSym<double, 2> > Matrix2;
0012 typedef ROOT::Math::SMatrix<double, 3, 3, ROOT::Math::MatRepSym<double, 3> > Matrix3;
0013
0014 void finvert(Matrix2& mm) {
0015   auto m = mm.Array();
0016
0017   auto c0 = 1 / m[0];
0018   auto c1 = m[1] * m[1] * c0;
0019   auto c2 = 1 / (m[2] - c1);
0020
0021   auto li21 = c1 * c0 * c2;
0022   m[0] = li21 + c0;
0023   m[1] = -m[1] * c0 * c2;
0024   m[2] = c2;
0025 }
0026
0027 void st() {}
0028 void en() {}
0029
0030 int main(int argc, char* argv[]) {
0031   double v[3] = {1., -0.2, 0.5};
0032   Matrix2& m = *(new Matrix2(v, 3));
0033
0034   std::cout << m << std::endl;
0035   invertPosDefMatrix(m);
0036   std::cout << m << std::endl;
0037   invertPosDefMatrix(m);
0038   std::cout << m << std::endl;
0039   m.Invert();
0040   std::cout << m << std::endl;
0041   m.Invert();
0042   std::cout << m << std::endl;
0043
0044   finvert(m);
0045   std::cout << m << std::endl;
0046   finvert(m);
0047   std::cout << m << std::endl;
0048
0049   if (argc > 1) {
0050     {
0051       edm::HRTimeType s = edm::hrRealTime();
0052       st();
0053       invertPosDefMatrix(m);
0054       en();
0055       edm::HRTimeType e = edm::hrRealTime();
0056       std::cout << e - s << std::endl;
0057     }
0058
0059     {
0060       edm::HRTimeType s = edm::hrRealTime();
0061       st();
0062       invertPosDefMatrix(m);
0063       en();
0064       edm::HRTimeType e = edm::hrRealTime();
0065       std::cout << e - s << std::endl;
0066     }
0067   } else {
0068     {
0069       edm::HRTimeType s = edm::hrRealTime();
0070       st();
0071       m.Invert();
0072       en();
0073       edm::HRTimeType e = edm::hrRealTime();
0074       std::cout << e - s << std::endl;
0075     }
0076
0077     {
0078       edm::HRTimeType s = edm::hrRealTime();
0079       st();
0080       m.Invert();
0081       en();
0082       edm::HRTimeType e = edm::hrRealTime();
0083       std::cout << e - s << std::endl;
0084     }
0085   }
0086   return 0;
0087 }