Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:02:33

0001 #ifndef RecoTracker_MkFitCore_src_Matriplex_MatriplexVector_h
0002 #define RecoTracker_MkFitCore_src_Matriplex_MatriplexVector_h
0003 
0004 #include "Matriplex.h"
0005 #include "Memory.h"
0006 #include <vector>
0007 #include <cassert>
0008 
0009 namespace Matriplex {
0010 
0011   //------------------------------------------------------------------------------
0012 
0013   template <class MP>
0014   class MatriplexVector {
0015     MP* fV;
0016     const idx_t fN;
0017 
0018     typedef typename MP::value_type T;
0019 
0020   public:
0021     MatriplexVector(idx_t n) : fN(n) { fV = (MP*)aligned_alloc64(sizeof(MP) * fN); }
0022 
0023     ~MatriplexVector() { std::free(fV); }
0024 
0025     idx_t size() const { return fN; }
0026 
0027     MP& mplex(int i) { return fV[i]; }
0028     MP& operator[](int i) { return fV[i]; }
0029 
0030     const MP& mplex(int i) const { return fV[i]; }
0031     const MP& operator[](int i) const { return fV[i]; }
0032 
0033     void setVal(T v) {
0034       for (idx_t i = 0; i < kTotSize; ++i) {
0035         fArray[i] = v;
0036       }
0037     }
0038 
0039     T& At(idx_t n, idx_t i, idx_t j) { return fV[n / fN].At(n % fN, i, j); }
0040 
0041     T& operator()(idx_t n, idx_t i, idx_t j) { return fV[n / fN].At(n % fN, i, j); }
0042 
0043     void copyIn(idx_t n, T* arr) { fV[n / fN].copyIn(n % fN, arr); }
0044     void copyOut(idx_t n, T* arr) { fV[n / fN].copyOut(n % fN, arr); }
0045   };
0046 
0047   template <class MP>
0048   using MPlexVec = MatriplexVector<MP>;
0049 
0050   //==============================================================================
0051 
0052   template <typename T, idx_t D1, idx_t D2, idx_t D3, idx_t N>
0053   void multiply(const MPlexVec<MPlex<T, D1, D2, N>>& A,
0054                 const MPlexVec<MPlex<T, D2, D3, N>>& B,
0055                 MPlexVec<MPlex<T, D1, D3, N>>& C,
0056                 int n_to_process = 0) {
0057     assert(A.size() == B.size());
0058     assert(A.size() == C.size());
0059 
0060     const int np = n_to_process ? n_to_process : A.size();
0061 
0062     for (int i = 0; i < np; ++i) {
0063       multiply(A[i], B[i], C[i]);
0064     }
0065   }
0066 
0067   template <typename T, idx_t D1, idx_t D2, idx_t D3, idx_t N>
0068   void multiplyGeneral(const MPlexVec<MPlex<T, D1, D2, N>>& A,
0069                        const MPlexVec<MPlex<T, D2, D3, N>>& B,
0070                        MPlexVec<MPlex<T, D1, D3, N>>& C,
0071                        int n_to_process = 0) {
0072     assert(A.size() == B.size());
0073     assert(A.size() == C.size());
0074 
0075     const int np = n_to_process ? n_to_process : A.size();
0076 
0077     for (int i = 0; i < np; ++i) {
0078       multiplyGeneral(A[i], B[i], C[i]);
0079     }
0080   }
0081 
0082   template <typename T, idx_t D1, idx_t D2, idx_t D3, idx_t N>
0083   void multiply3in(MPlexVec<MPlex<T, D1, D2, N>>& A,
0084                    MPlexVec<MPlex<T, D2, D3, N>>& B,
0085                    MPlexVec<MPlex<T, D1, D3, N>>& C,
0086                    int n_to_process = 0) {
0087     assert(A.size() == B.size());
0088     assert(A.size() == C.size());
0089 
0090     const int np = n_to_process ? n_to_process : A.size();
0091 
0092     for (int i = 0; i < np; ++i) {
0093       multiply(A[i], B[i], C[i]);
0094       multiply(B[i], C[i], A[i]);
0095       multiply(C[i], A[i], B[i]);
0096     }
0097   }
0098 
0099   template <typename T, idx_t D, idx_t N>
0100   void multiply(const MPlexVec<MPlexSym<T, D, N>>& A,
0101                 const MPlexVec<MPlexSym<T, D, N>>& B,
0102                 MPlexVec<MPlex<T, D, D, N>>& C,
0103                 int n_to_process = 0) {
0104     assert(A.size() == B.size());
0105     assert(A.size() == C.size());
0106 
0107     const int np = n_to_process ? n_to_process : A.size();
0108 
0109     for (int i = 0; i < np; ++i) {
0110       multiply(A[i], B[i], C[i]);
0111     }
0112   }
0113 
0114   //==============================================================================
0115 
0116   template <typename T, idx_t D, idx_t N>
0117   void invertCramer(MPlexVec<MPlex<T, D, D, N>>& A, int n_to_process = 0) {
0118     const int np = n_to_process ? n_to_process : A.size();
0119 
0120     for (int i = 0; i < np; ++i) {
0121       invertCramer(A[i]);
0122     }
0123   }
0124 
0125   template <typename T, idx_t D, idx_t N>
0126   void invertCholesky(MPlexVec<MPlex<T, D, D, N>>& A, int n_to_process = 0) {
0127     const int np = n_to_process ? n_to_process : A.size();
0128 
0129     for (int i = 0; i < np; ++i) {
0130       invertCholesky(A[i]);
0131     }
0132   }
0133 
0134   template <typename T, idx_t D, idx_t N>
0135   void invertCramerSym(MPlexVec<MPlexSym<T, D, N>>& A, int n_to_process = 0) {
0136     const int np = n_to_process ? n_to_process : A.size();
0137 
0138     for (int i = 0; i < np; ++i) {
0139       invertCramerSym(A[i]);
0140     }
0141   }
0142 
0143   template <typename T, idx_t D, idx_t N>
0144   void invertCholeskySym(MPlexVec<MPlexSym<T, D, N>>& A, int n_to_process = 0) {
0145     const int np = n_to_process ? n_to_process : A.size();
0146 
0147     for (int i = 0; i < np; ++i) {
0148       invertCholeskySym(A[i]);
0149     }
0150   }
0151 
0152 }  // namespace Matriplex
0153 
0154 #endif