Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 12:54:04

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