src/Matriplex/MatriplexVector.h

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