Merge pull request #68 from davidlt/add-root6-smatrix-patch

Add SMatrix patch to ROOT6

Author: ktf

math/smatrix/inc/Math/MatrixRepresentationsStatic.h

@@ -23,6 +23,11 @@
 #include "Math/StaticCheck.h"
 #endif
 
+#include <cstddef>
+#include <utility>
+#include <type_traits>
+#include <array>
+
 namespace ROOT {
 
 namespace Math {
@@ -142,96 +147,49 @@ namespace Math {
       int fOff[D*D];
    };
 
-// Make the lookup tables available at compile time:
-// Add them to a namespace?
-static const int fOff1x1[] = {0};
-static const int fOff2x2[] = {0, 1, 1, 2};
-static const int fOff3x3[] = {0, 1, 3, 1, 2, 4, 3, 4, 5};
-static const int fOff4x4[] = {0, 1, 3, 6, 1, 2, 4, 7, 3, 4, 5, 8, 6, 7, 8, 9};
-static const int fOff5x5[] = {0, 1, 3, 6, 10, 1, 2, 4, 7, 11, 3, 4, 5, 8, 12, 6, 7, 8, 9, 13, 10, 11, 12, 13, 14};
-static const int fOff6x6[] = {0, 1, 3, 6, 10, 15, 1, 2, 4, 7, 11, 16, 3, 4, 5, 8, 12, 17, 6, 7, 8, 9, 13, 18, 10, 11, 12, 13, 14, 19, 15, 16, 17, 18, 19, 20};
+  namespace rowOffsetsUtils {
 
-static const int fOff7x7[] = {0, 1, 3, 6, 10, 15, 21, 1, 2, 4, 7, 11, 16, 22, 3, 4, 5, 8, 12, 17, 23, 6, 7, 8, 9, 13, 18, 24, 10, 11, 12, 13, 14, 19, 25, 15, 16, 17, 18, 19, 20, 26, 21, 22, 23, 24, 25, 26, 27};
+    ///////////
+    // Some meta template stuff
+    template<int...> struct indices{};
 
-static const int fOff8x8[] = {0, 1, 3, 6, 10, 15, 21, 28, 1, 2, 4, 7, 11, 16, 22, 29, 3, 4, 5, 8, 12, 17, 23, 30, 6, 7, 8, 9, 13, 18, 24, 31, 10, 11, 12, 13, 14, 19, 25, 32, 15, 16, 17, 18, 19, 20, 26, 33, 21, 22, 23, 24, 25, 26, 27, 34, 28, 29, 30, 31, 32, 33, 34, 35};
+    template<int I, class IndexTuple, int N>
+    struct make_indices_impl;
 
-static const int fOff9x9[] = {0, 1, 3, 6, 10, 15, 21, 28, 36, 1, 2, 4, 7, 11, 16, 22, 29, 37, 3, 4, 5, 8, 12, 17, 23, 30, 38, 6, 7, 8, 9, 13, 18, 24, 31, 39, 10, 11, 12, 13, 14, 19, 25, 32, 40, 15, 16, 17, 18, 19, 20, 26, 33, 41, 21, 22, 23, 24, 25, 26, 27, 34, 42, 28, 29, 30, 31, 32, 33, 34, 35, 43, 36, 37, 38, 39, 40, 41, 42, 43, 44};
+    template<int I, int... Indices, int N>
+    struct make_indices_impl<I, indices<Indices...>, N>
+    {
+      typedef typename make_indices_impl<I + 1, indices<Indices..., I>,
+					 N>::type type;
+    };
 
-static const int fOff10x10[] = {0, 1, 3, 6, 10, 15, 21, 28, 36, 45, 1, 2, 4, 7, 11, 16, 22, 29, 37, 46, 3, 4, 5, 8, 12, 17, 23, 30, 38, 47, 6, 7, 8, 9, 13, 18, 24, 31, 39, 48, 10, 11, 12, 13, 14, 19, 25, 32, 40, 49, 15, 16, 17, 18, 19, 20, 26, 33, 41, 50, 21, 22, 23, 24, 25, 26, 27, 34, 42, 51, 28, 29, 30, 31, 32, 33, 34, 35, 43, 52, 36, 37, 38, 39, 40, 41, 42, 43, 44, 53, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54};
+    template<int N, int... Indices>
+    struct make_indices_impl<N, indices<Indices...>, N> {
+      typedef indices<Indices...> type;
+    };
 
-template<>
-   struct RowOffsets<1> {
-      RowOffsets() {}
-      int operator()(unsigned int , unsigned int ) const { return 0; } // Just one element
-      int apply(unsigned int ) const { return 0; }
-   };
+    template<int N>
+    struct make_indices : make_indices_impl<0, indices<>, N> {};
+    // end of stuff
 
-template<>
-   struct RowOffsets<2> {
-      RowOffsets() {}
-      int operator()(unsigned int i, unsigned int j) const { return i+j; /*fOff2x2[i*2+j];*/ }
-      int apply(unsigned int i) const { return fOff2x2[i]; }
-   };
 
-template<>
-   struct RowOffsets<3> {
-      RowOffsets() {}
-      int operator()(unsigned int i, unsigned int j) const { return fOff3x3[i*3+j]; }
-      int apply(unsigned int i) const { return fOff3x3[i]; }
-   };
 
-template<>
-   struct RowOffsets<4> {
-     RowOffsets() {}
-     int operator()(unsigned int i, unsigned int j) const { return fOff4x4[i*4+j]; }
-     int apply(unsigned int i) const { return fOff4x4[i]; }
-   };
+    template<int I0, class F, int... I>
+    constexpr std::array<decltype(std::declval<F>()(std::declval<int>())), sizeof...(I)>
+    do_make(F f, indices<I...>)
+    {
+      return  std::array<decltype(std::declval<F>()(std::declval<int>())),
+			 sizeof...(I)>{{ f(I0 + I)... }};
+    }
 
-   template<>
-   struct RowOffsets<5> {
-     inline RowOffsets() {}
-     inline int operator()(unsigned int i, unsigned int j) const { return fOff5x5[i*5+j]; }
-//   int operator()(unsigned int i, unsigned int j) const {
-//     if(j <= i) return (i * (i + 1)) / 2 + j;
-//      else return (j * (j + 1)) / 2 + i;
-//     }
-   inline int apply(unsigned int i) const { return fOff5x5[i]; }
-   };
+    template<int N, int I0 = 0, class F>
+    constexpr std::array<decltype(std::declval<F>()(std::declval<int>())), N>
+    make(F f) {
+      return do_make<I0>(f, typename make_indices<N>::type());
+    }
 
-template<>
-   struct RowOffsets<6> {
-     RowOffsets() {}
-     int operator()(unsigned int i, unsigned int j) const { return fOff6x6[i*6+j]; }
-     int apply(unsigned int i) const { return fOff6x6[i]; }
-   };
+  } // namespace rowOffsetsUtils
 
-template<>
-   struct RowOffsets<7> {
-     RowOffsets() {}
-     int operator()(unsigned int i, unsigned int j) const { return fOff7x7[i*7+j]; }
-     int apply(unsigned int i) const { return fOff7x7[i]; }
-   };
-
-template<>
-   struct RowOffsets<8> {
-     RowOffsets() {}
-     int operator()(unsigned int i, unsigned int j) const { return fOff8x8[i*8+j]; }
-     int apply(unsigned int i) const { return fOff8x8[i]; }
-   };
-
-template<>
-   struct RowOffsets<9> {
-     RowOffsets() {}
-     int operator()(unsigned int i, unsigned int j) const { return fOff9x9[i*9+j]; }
-     int apply(unsigned int i) const { return fOff9x9[i]; }
-   };
-
-template<>
-   struct RowOffsets<10> {
-     RowOffsets() {}
-     int operator()(unsigned int i, unsigned int j) const { return fOff10x10[i*10+j]; }
-     int apply(unsigned int i) const { return fOff10x10[i]; }
-   };
 
 //_________________________________________________________________________________
    /**
@@ -257,35 +215,32 @@ template<>
 
    public:
 
-      MatRepSym() :fOff(0) { CreateOffsets(); }
+    /* constexpr */ inline MatRepSym(){}
 
-      typedef T  value_type;
+    typedef T  value_type;
 
-      inline const T& operator()(unsigned int i, unsigned int j) const {
-         return fArray[Offsets()(i,j)];
-      }
-      inline T& operator()(unsigned int i, unsigned int j) {
-         return fArray[Offsets()(i,j)];
-      }
 
-      inline T& operator[](unsigned int i) {
-         return fArray[Offsets().apply(i) ];
-//return fArray[Offsets()(i/D, i%D)];
-      }
+    inline T & operator()(unsigned int i, unsigned int j)
+     { return fArray[offset(i, j)]; }
 
-      inline const T& operator[](unsigned int i) const {
-         return fArray[Offsets().apply(i) ];
-//return fArray[Offsets()(i/D, i%D)];
-      }
+     inline /* constexpr */ T const & operator()(unsigned int i, unsigned int j) const
+     { return fArray[offset(i, j)]; }
 
-      inline T apply(unsigned int i) const {
-         return fArray[Offsets().apply(i) ];
-         //return operator()(i/D, i%D);
-      }
+     inline T& operator[](unsigned int i) {
+       return fArray[off(i)];
+     }
 
-      inline T* Array() { return fArray; }
+     inline /* constexpr */ T const & operator[](unsigned int i) const {
+       return fArray[off(i)];
+     }
 
-      inline const T* Array() const { return fArray; }
+     inline /* constexpr */ T apply(unsigned int i) const {
+       return fArray[off(i)];
+     }
+
+     inline T* Array() { return fArray; }
+
+     inline const T* Array() const { return fArray; }
 
       /**
          assignment : only symmetric to symmetric allowed
@@ -346,22 +301,26 @@ template<>
          kSize = D*(D+1)/2
       };
 
+     static constexpr int off0(int i) { return i==0 ? 0 : off0(i-1)+i;}
+     static constexpr int off2(int i, int j) { return j<i ? off0(i)+j : off0(j)+i; }
+     static constexpr int off1(int i) { return off2(i/D, i%D);}
 
-      void CreateOffsets() {
-         const static RowOffsets<D> off;
-         fOff = &off;
-      }
+     static int off(int i) {
+       static constexpr auto v = rowOffsetsUtils::make<D*D>(off1);
+       return v[i];
+     }
 
-      inline const RowOffsets<D> & Offsets() const {
-         return *fOff;
-      }
+     static inline constexpr unsigned int
+     offset(unsigned int i, unsigned int j)
+     {
+       //if (j > i) std::swap(i, j);
+       return off(i*D+j);
+       // return (i>j) ? (i * (i+1) / 2) + j :  (j * (j+1) / 2) + i;
+     }
 
    private:
       //T __attribute__ ((aligned (16))) fArray[kSize];
       T fArray[kSize];
-
-      const RowOffsets<D> * fOff;   //! transient
-
    };
 
 

math/smatrix/inc/Math/SMatrix.h

@@ -97,6 +97,7 @@ namespace Math {
 template <class T, unsigned int D> class SVector;
 
 struct SMatrixIdentity { };
+struct SMatrixNoInit { };
 
 //__________________________________________________________________________
 /**
@@ -149,6 +150,11 @@ class SMatrix {
    */
    SMatrix();
    ///
+  /**
+       construct from without initialization
+   */
+   inline SMatrix( SMatrixNoInit ){}
+
    /**
        construct from an identity matrix
    */