Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 #ifndef DataFormat_Math_ProjectMatrix_H
0002 #define DataFormat_Math_ProjectMatrix_H
0003 
0004 #define SMATRIX_USE_CONSTEXPR
0005 #include "Math/SMatrix.h"
0006 
0007 template <typename T, unsigned int N, unsigned int D>
0008 struct ProjectMatrix {
0009   typedef ROOT::Math::SMatrix<T, D, D, ROOT::Math::MatRepSym<T, D> > SMatDD;
0010   typedef ROOT::Math::SMatrix<T, N, N> SMatNN;
0011   typedef ROOT::Math::SMatrix<T, N, D> SMatND;
0012   typedef ROOT::Math::SMatrix<T, D, N> SMatDN;
0013 
0014   // no constructor
0015 
0016   // H
0017   SMatDN matrix() const {
0018     SMatDN r;
0019     for (unsigned int i = 0; i < D; i++)
0020       r(i, index[i]) = T(1.);
0021     return r;
0022   }
0023 
0024   // constuct given H
0025   void fromH(SMatDN const& H) {
0026     for (unsigned int i = 0; i < D; i++) {
0027       index[i] = N;
0028       for (unsigned int j = 0; j < N; j++)
0029         if (H(i, j) > 0) {
0030           index[i] = j;
0031           break;
0032         }
0033     }
0034   }
0035 
0036   // H*S
0037   SMatND project(SMatDD const& s) const {
0038     SMatND r;
0039     for (unsigned int i = 0; i < D; i++)
0040       for (unsigned int j = 0; j < D; j++)
0041         r(index[i], j) = s(i, j);
0042     return r;
0043   }
0044 
0045   // K*H
0046   SMatNN project(SMatND const& k) const {
0047     SMatNN s;
0048     for (unsigned int i = 0; i < N; i++)
0049       for (unsigned int j = 0; j < D; j++)
0050         s(i, index[j]) = k(i, j);
0051     return s;
0052   }
0053 
0054   // S-K*H
0055   void projectAndSubtractFrom(SMatNN& __restrict__ s, SMatND const& __restrict__ k) const {
0056     for (unsigned int i = 0; i < N; i++)
0057       for (unsigned int j = 0; j < D; j++)
0058         s(i, index[j]) -= k(i, j);
0059   }
0060 
0061   // only H(i,index[i])=1.
0062   unsigned int index[D];
0063 };
0064 
0065 #endif