Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-04 22:54:37 

0001 #ifdef NEWSMATRIX
0002 // #include "DataFormats/Math/interface/MulSymMatrix.h"
0003 #define SMATRIX_USE_CONSTEXPR
0004 #include "Math/SMatrix.h"
0005 
0006 typedef unsigned int IndexType;
0007 //typedef unsigned long long IndexType;
0008 
0009 namespace ROOT {
0010 
0011   namespace Math {
0012 
0013     constexpr int fOff15x15[] = {
0014         0,   1,   3,   6,   10,  15,  21,  28,  36,  45,  55,  66,  78,  91,  105, 1,   2,   4,   7,   11,  16,
0015         22,  29,  37,  46,  56,  67,  79,  92,  106, 3,   4,   5,   8,   12,  17,  23,  30,  38,  47,  57,  68,
0016         80,  93,  107, 6,   7,   8,   9,   13,  18,  24,  31,  39,  48,  58,  69,  81,  94,  108, 10,  11,  12,
0017         13,  14,  19,  25,  32,  40,  49,  59,  70,  82,  95,  109, 15,  16,  17,  18,  19,  20,  26,  33,  41,
0018         50,  60,  71,  83,  96,  110, 21,  22,  23,  24,  25,  26,  27,  34,  42,  51,  61,  72,  84,  97,  111,
0019         28,  29,  30,  31,  32,  33,  34,  35,  43,  52,  62,  73,  85,  98,  112, 36,  37,  38,  39,  40,  41,
0020         42,  43,  44,  53,  63,  74,  86,  99,  113, 45,  46,  47,  48,  49,  50,  51,  52,  53,  54,  64,  75,
0021         87,  100, 114, 55,  56,  57,  58,  59,  60,  61,  62,  63,  64,  65,  76,  88,  101, 115, 66,  67,  68,
0022         69,  70,  71,  72,  73,  74,  75,  76,  77,  89,  102, 116, 78,  79,  80,  81,  82,  83,  84,  85,  86,
0023         87,  88,  89,  90,  103, 117, 91,  92,  93,  94,  95,  96,  97,  98,  99,  100, 101, 102, 103, 104, 118,
0024         105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119};
0025 
0026     constexpr int const* pfOff15x15[] = {fOff15x15 + 0,
0027                                          fOff15x15 + 15,
0028                                          fOff15x15 + 2 * 15,
0029                                          fOff15x15 + 3 * 15,
0030                                          fOff15x15 + 4 * 15,
0031                                          fOff15x15 + 5 * 15,
0032                                          fOff15x15 + 6 * 15,
0033                                          fOff15x15 + 7 * 15,
0034                                          fOff15x15 + 8 * 15,
0035                                          fOff15x15 + 9 * 15,
0036                                          fOff15x15 + 10 * 15,
0037                                          fOff15x15 + 11 * 15,
0038                                          fOff15x15 + 12 * 15,
0039                                          fOff15x15 + 13 * 15,
0040                                          fOff15x15 + 14 * 15};
0041     template <>
0042     struct RowOffsets<15> {
0043       RowOffsets() {}
0044       int operator()(unsigned int i, unsigned int j) const { return pfOff15x15[i][j]; }  // fOff15x15[i*15+j]; }
0045       int apply(unsigned int i) const { return fOff15x15[i]; }
0046     };
0047 
0048     /** 
0049        Force Expression evaluation from general to symmetric. 
0050        To be used when is known (like in similarity products) that the result 
0051        is symmetric
0052        Note this is function used in the simmilarity product: no check for temporary is 
0053        done since in that case is not needed
0054    */
0055     struct AssignAsSym {
0056       /// assign a symmetric matrix from an expression
0057       template <class T, IndexType D, class A, class R>
0058       static void Evaluate(SMatrix<T, D, D, MatRepStd<T, D> >& lhs, const Expr<A, T, D, D, R>& rhs) {
0059         //for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = rhs.apply(i);
0060         for (IndexType i = 0; i < D; ++i)
0061           // storage of symmetric matrix is in lower block
0062           for (IndexType j = 0; j <= i; ++j) {
0063             lhs(i, j) = rhs(i, j);
0064           }
0065       }
0066       /// assign the symmetric matric  from a general matrix
0067       template <class T, IndexType D, class R>
0068       static void Evaluate(SMatrix<T, D, D, MatRepStd<T, D> >& lhs, const SMatrix<T, D, D, R>& rhs) {
0069         //for(unsigned int i=0; i<D1*D2; ++i) lhs.fRep[i] = rhs.apply(i);
0070         for (IndexType i = 0; i < D; ++i)
0071           // storage of symmetric matrix is in lower block
0072           for (IndexType j = 0; j <= i; ++j) {
0073             lhs(i, j) = rhs(i, j);
0074           }
0075       }
0076     };  // struct AssignSym
0077 
0078     template <class T, IndexType D1, IndexType D2, class R>
0079     inline SMatrix<T, D1, D1, MatRepSym<T, D1> > Similarity(const SMatrix<T, D1, D2, R>& lhs,
0080                                                             const SMatrix<T, D2, D2, MatRepStd<T, D2> >& rhs) {
0081       SMatrix<T, D1, D2, MatRepStd<T, D1, D2> > tmp = lhs * rhs;
0082       typedef SMatrix<T, D1, D1, MatRepSym<T, D1> > SMatrixSym;
0083       SMatrixSym mret;
0084       AssignSym::Evaluate(mret, tmp * Transpose(lhs));
0085       return mret;
0086     }
0087 
0088   }  // namespace Math
0089 }  // namespace ROOT
0090 
0091 template <typename T, IndexType N>
0092 inline void mult(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> >& a,
0093                  ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> > const& rh,
0094                  ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> > const& lh) {
0095   // a(i,j) = r(i,k)*l(k,j)
0096   for (IndexType i = 0; i != N; ++i) {
0097     IndexType off_i = a.fRep.offset(i, 0);
0098     for (IndexType k = 0; k != N; ++k) {
0099       IndexType off_k = a.fRep.offset(k, 0);
0100       if (k < i) {
0101         for (IndexType j = 0; j != (k + 1); ++j)
0102           a.Array()[off_i + j] += rh(i, k) * lh.Array()[off_k + j];
0103         for (IndexType j = k + 1; j != (i + 1); ++j)
0104           a.Array()[off_i + j] += rh(i, k) * lh(k, j);
0105       } else
0106         for (IndexType j = 0; j != (i + 1); ++j)
0107           a.Array()[off_i + j] += rh(i, k) * lh.Array()[off_k + j];
0108     }
0109   }
0110 }
0111 
0112 template <typename T, IndexType N>
0113 inline void mult(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> >& a,
0114                  ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& rh,
0115                  ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& lh) {
0116   // a(i,j) = r(i,k)*l(k,j)
0117   for (IndexType i = 0; i != N; ++i) {
0118     for (IndexType j = 0; j <= i; ++j) {
0119       // a(i,j)=0;
0120       //    a(i,j) += rh(i,k)*lh(k,j);
0121       for (IndexType k = 0; k != N; ++k)
0122         a(i, j) += rh(i, k) * lh(j, k);
0123     }
0124   }
0125 
0126   for (IndexType i = 0; i != N - 1; ++i)
0127     for (IndexType j = i + 1; j != N; ++j)
0128       a(i, j) = a(j, i);
0129 }
0130 
0131 // U(i,k) * A(k,l) * U(j,l)
0132 template <typename T, IndexType N>
0133 inline void similarity(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> >& b,
0134                        ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& u,
0135                        ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& a) {
0136   for (IndexType i = 0; i != N; ++i)
0137     for (IndexType j = 0; j <= i; ++j)
0138       for (IndexType k = 0; k != N; ++k)
0139         for (IndexType l = 0; l != N; ++l)
0140           b(i, j) += u(i, k) * a(k, l) * u(j, l);
0141   /*
0142   T s[N];
0143   for (IndexType i=0; i!=N; ++i) 
0144     for (IndexType j=0; j<=i; ++j) { 
0145       for (IndexType k=0; k!=N; ++k) {
0146     s[k]=0;
0147     for (IndexType l=0; l!=N; ++l) 
0148       s[k] += a(k,l)*u(j,l);
0149       }
0150       for (IndexType k=0; k!=N; ++k)
0151     b(i,j) += u(i,k)*s[k];
0152     }
0153   */
0154 
0155   for (IndexType i = 0; i != N - 1; ++i)
0156     for (IndexType j = i + 1; j != N; ++j)
0157       b(i, j) = b(j, i);
0158 }
0159 
0160 // U(k,i) * A(k,l) * U(l,j)
0161 template <typename T, IndexType N>
0162 inline void similarityT(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> >& b,
0163                         ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& u,
0164                         ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > const& a) {
0165   for (IndexType i = 0; i != N; ++i)
0166     for (IndexType j = 0; j <= i; ++j)
0167       for (IndexType k = 0; k != N; ++k)
0168         for (IndexType l = 0; l != N; ++l)
0169           b(i, j) += u(k, i) * a(k, l) * u(l, j);
0170 
0171   for (IndexType i = 0; i != N - 1; ++i)
0172     for (IndexType j = i + 1; j != N; ++j)
0173       b(i, j) = b(j, i);
0174 }
0175 
0176 template <typename M1, typename M2>
0177 double eps(M1 const& m1, M2 const& m2) {
0178   IndexType N = M1::kRows;
0179   double ret = 0.;
0180   for (IndexType i = 0; i != N; ++i)
0181     for (IndexType j = 0; j != N; ++j)
0182       ret = std::max(ret, std::abs(m1(i, j) - m2(i, j)));
0183   return ret;
0184 }
0185 
0186 template <typename M>
0187 bool isSym(M const& m) {
0188   IndexType N = M::kRows;
0189   for (IndexType i = 0; i != N; ++i)
0190     for (IndexType j = 0; j <= i; ++j)
0191       if (m(i, j) != m(j, i))
0192         return false;
0193   return true;
0194 }
0195 
0196 #include <iostream>
0197 #include "FWCore/Utilities/interface/HRRealTime.h"
0198 
0199 #include <random>
0200 
0201 namespace {
0202   std::mt19937 eng;
0203   std::uniform_real_distribution<double> rgen(-5., 5.);
0204   template <typename T, IndexType N>
0205   inline void fillRandom(ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> >& a) {
0206     for (IndexType k = 0; k != N; ++k) {
0207       a(k, k) = std::abs(rgen(eng));
0208       for (IndexType j = 0; j < k; ++j)
0209         a(k, j) = rgen(eng);
0210     }
0211     for (IndexType i = 0; i != N - 1; ++i)
0212       for (IndexType j = i + 1; j != N; ++j)
0213         a(i, j) = -a(j, i);
0214   }
0215 }  // namespace
0216 
0217 template <typename T, IndexType N>
0218 void go(edm::HRTimeType& s1,
0219         edm::HRTimeType& s2,
0220         edm::HRTimeType& s3,
0221         edm::HRTimeType& s4,
0222         edm::HRTimeType& s5,
0223         edm::HRTimeType& s6,
0224         edm::HRTimeType& s7,
0225         edm::HRTimeType& s8,
0226         bool print) {
0227   typedef ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepStd<T, N> > Matrix;
0228   typedef ROOT::Math::SMatrix<T, N, N, ROOT::Math::MatRepSym<T, N> > SymMatrix;
0229 
0230   ROOT::Math::SMatrixIdentity id;
0231 
0232   Matrix im1(id);
0233   Matrix im2(id);
0234   /*
0235   fillRandom(im1);
0236   im1 = im1*ROOT::Math::Transpose(im1);
0237   im1 *=im1;
0238 
0239   Matrix im2; fillRandom(im2);
0240   im2 = im2*ROOT::Math::Transpose(im2);
0241   im2 *=im2;
0242   */
0243 
0244   SymMatrix is1;
0245   ROOT::Math::AssignSym::Evaluate(is1, im1);
0246   SymMatrix is2;
0247   ROOT::Math::AssignSym::Evaluate(is2, im2);
0248 
0249   Matrix j1;
0250   fillRandom(j1);
0251   Matrix j2;
0252   fillRandom(j2);
0253 
0254   SymMatrix rh;
0255   s6 = edm::hrRealTime();
0256   rh = ROOT::Math::Similarity(j1, is1);
0257   s6 = edm::hrRealTime() - s6;
0258 
0259   SymMatrix rh2;
0260   s7 = edm::hrRealTime();
0261   rh2 = ROOT::Math::Similarity(j1, im1);
0262   s7 = edm::hrRealTime() - s7;
0263 
0264   Matrix mrh;
0265   s8 = edm::hrRealTime();
0266   similarity(mrh, j1, im1);
0267   s8 = edm::hrRealTime() - s8;
0268 
0269   SymMatrix lh = ROOT::Math::SimilarityT(j2, is1);
0270 
0271   Matrix mlh;
0272   similarityT(mlh, j2, im1);
0273 
0274   SymMatrix a;
0275   s1 = edm::hrRealTime();
0276   mult(a, lh, rh);
0277   s1 = edm::hrRealTime() - s1;
0278 
0279   SymMatrix b;
0280   s2 = edm::hrRealTime();
0281   ROOT::Math::AssignSym::Evaluate(b, lh * rh);
0282   s2 = edm::hrRealTime() - s2;
0283 
0284   Matrix m0 = b;
0285 
0286   Matrix m;
0287   s3 = edm::hrRealTime();
0288   m = mlh * mrh;
0289   s3 = edm::hrRealTime() - s3;
0290 
0291   SymMatrix sm;
0292   s4 = edm::hrRealTime();
0293   ROOT::Math::AssignSym::Evaluate(sm, mlh * mrh);
0294   s4 = edm::hrRealTime() - s4;
0295 
0296   Matrix m2;
0297   s5 = edm::hrRealTime();
0298   mult(m2, mlh, mrh);
0299   s5 = edm::hrRealTime() - s5;
0300 
0301   SymMatrix smm;
0302   ROOT::Math::AssignSym::Evaluate(smm, m);
0303   SymMatrix smm2;
0304   ROOT::Math::AssignSym::Evaluate(smm2, m2);
0305 
0306   if (print) {
0307     if (!isSym(im1))
0308       std::cout << " im is not sym" << std::endl;
0309     if (!isSym(mrh))
0310       std::cout << " rh is not sym" << std::endl;
0311     if (!isSym(mlh))
0312       std::cout << " lh is not sym" << std::endl;
0313     if (!isSym(m))
0314       std::cout << " m is not sym" << std::endl;
0315     if (!isSym(m2))
0316       std::cout << "m2 is not sym" << std::endl;
0317 
0318     std::cout << "eps sim  " << eps(rh, mrh) << std::endl;
0319     std::cout << "eps sim  " << eps(rh, rh2) << std::endl;
0320     std::cout << "eps simT " << eps(lh, mlh) << std::endl;
0321     std::cout << "eps s m  " << eps(m, b) << std::endl;
0322     std::cout << "eps s sm " << eps(m2, b) << std::endl;
0323 
0324     // std::cout << b << std::endl;
0325     // std::cout << m << std::endl;
0326 
0327     if (m != m0)
0328       std::cout << "problem with SMatrix Assign " << eps(m, m0) << std::endl;
0329     if (smm != b)
0330       std::cout << "problem with SMatrix * " << eps(smm, b) << std::endl;
0331     if (sm != b)
0332       std::cout << "problem with SMatrix evaluate " << eps(sm, b) << std::endl;
0333     if (a != b)
0334       std::cout << "problem with MulSymMatrix " << eps(a, b) << std::endl;
0335     if (a != smm)
0336       std::cout << "problem with MulSymMatrix twice " << eps(a, smm) << std::endl;
0337     if (m != m2)
0338       std::cout << "problem with MulMatrix " << eps(m, m2) << std::endl;
0339     if (smm != smm2)
0340       std::cout << "problem with MulMatrix twice " << eps(smm, smm2) << std::endl;
0341 
0342     std::cout << "sym mult   " << s1 << std::endl;
0343     std::cout << "sym    *   " << s2 << std::endl;
0344     std::cout << "mat    *   " << s3 << std::endl;
0345     std::cout << "mat as sym " << s4 << std::endl;
0346     std::cout << "mat mult   " << s5 << std::endl;
0347     std::cout << "sym  sim   " << s6 << std::endl;
0348     std::cout << "std  sim   " << s7 << std::endl;
0349     std::cout << "loop sim   " << s8 << std::endl;
0350   }
0351 }
0352 
0353 int main() {
0354   edm::HRTimeType s1 = 0, s2 = 0, s3 = 0, s4 = 0, s5 = 0, s6 = 0, s7 = 0, s8 = 0;
0355   go<double, 3>(s1, s2, s3, s4, s5, s6, s7, s8, true);
0356 
0357   go<double, 15>(s1, s2, s3, s4, s5, s6, s7, s8, true);
0358 
0359   edm::HRTimeType t1 = 0;
0360   edm::HRTimeType t2 = 0;
0361   edm::HRTimeType t3 = 0;
0362   edm::HRTimeType t4 = 0;
0363   edm::HRTimeType t5 = 0;
0364   edm::HRTimeType t6 = 0;
0365   edm::HRTimeType t7 = 0;
0366   edm::HRTimeType t8 = 0;
0367 
0368   for (int i = 0; i != 50000; ++i) {
0369     go<double, 15>(s1, s2, s3, s4, s5, s6, s7, s8, false);
0370     t1 += s1;
0371     t2 += s2;
0372     t3 += s3;
0373     t4 += s4;
0374     t5 += s5;
0375     t6 += s6;
0376     t7 += s7;
0377     t8 += s8;
0378   }
0379   std::cout << "sym mult   " << t1 / 50000 << std::endl;
0380   std::cout << "sym    *   " << t2 / 50000 << std::endl;
0381   std::cout << "mat    *   " << t3 / 50000 << std::endl;
0382   std::cout << "mat as sym " << t4 / 50000 << std::endl;
0383   std::cout << "mat mult   " << t5 / 50000 << std::endl;
0384   std::cout << "sym  sim " << t6 / 50000 << std::endl;
0385   std::cout << "std  sim " << t7 / 50000 << std::endl;
0386   std::cout << "loop sim " << t8 / 50000 << std::endl;
0387 
0388   return 0;
0389 }
0390 #else
0391 int main() { return 0; }
0392 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!