Math/test/MulSymMatrix_t.cpp

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