root/Math/MatrixRepresentationsStatic.h

0001 // @(#)root/smatrix:$Id$
0002 // Author: L. Moneta, J. Palacios    2006
0003
0004 #ifndef ROOT_Math_MatrixRepresentationsStatic
0005 #define ROOT_Math_MatrixRepresentationsStatic 1
0006
0007 // Include files
0008
0009 /**
0010 \defgroup MatRep SMatrix Storage Representation
0011 \ingroup SMatrixGroup
0012
0013 Classes MatRepStd and MatRepSym for generic and symmetric matrix
0014 data storage and manipulation. Define data storage and access, plus
0015 operators =, +=, -=, ==.
0016
0017 \author Juan Palacios
0018 \date   2006-01-15
0019  */
0020
0021 #include "Math/StaticCheck.h"
0022
0023 #include <cstddef>
0024 #include <utility>
0025 #include <type_traits>
0026 #include <array>
0027
0028 namespace ROOT {
0029
0030 namespace Math {
0031
0032 /**
0033 \defgroup MatRepStd Standard Matrix representation
0034 \ingroup MatRep
0035
0036 Standard Matrix representation for a general D1 x D2 matrix.
0037 This class is itself a template on the contained type T, the number of rows and the number of columns.
0038 Its data member is an array T[nrows*ncols] containing the matrix data.
0039 The data are stored in the row-major C convention.
0040 For example, for a matrix, M, of size 3x3, the data \f$ \left[a_0,a_1,a_2,.......,a_7,a_8 \right] \f$d
0041 are stored in the following order:
0042
0043 \f[
0044 M = \left( \begin{array}{ccc}
0045 a_0 & a_1 & a_2  \\
0046 a_3 & a_4  & a_5  \\
0047 a_6 & a_7  & a_8   \end{array} \right)
0048 \f]
0049
0050 */
0051
0052
0053    template <class T, unsigned int D1, unsigned int D2=D1>
0054    class MatRepStd {
0055
0056    public:
0057
0058       typedef T  value_type;
0059
0060       inline const T& operator()(unsigned int i, unsigned int j) const {
0061          return fArray[i*D2+j];
0062       }
0063       inline T& operator()(unsigned int i, unsigned int j) {
0064          return fArray[i*D2+j];
0065       }
0066       inline T& operator[](unsigned int i) { return fArray[i]; }
0067
0068       inline const T& operator[](unsigned int i) const { return fArray[i]; }
0069
0070       inline T apply(unsigned int i) const { return fArray[i]; }
0071
0072       inline T* Array() { return fArray; }
0073
0074       inline const T* Array() const { return fArray; }
0075
0076       template <class R>
0077       inline MatRepStd<T, D1, D2>& operator+=(const R& rhs) {
0078          for(unsigned int i=0; i<kSize; ++i) fArray[i] += rhs[i];
0079          return *this;
0080       }
0081
0082       template <class R>
0083       inline MatRepStd<T, D1, D2>& operator-=(const R& rhs) {
0084          for(unsigned int i=0; i<kSize; ++i) fArray[i] -= rhs[i];
0085          return *this;
0086       }
0087
0088       template <class R>
0089       inline MatRepStd<T, D1, D2>& operator=(const R& rhs) {
0090          for(unsigned int i=0; i<kSize; ++i) fArray[i] = rhs[i];
0091          return *this;
0092       }
0093
0094       template <class R>
0095       inline bool operator==(const R& rhs) const {
0096          bool rc = true;
0097          for(unsigned int i=0; i<kSize; ++i) {
0098             rc = rc && (fArray[i] == rhs[i]);
0099          }
0100          return rc;
0101       }
0102
0103       enum {
0104          /// return no. of matrix rows
0105          kRows = D1,
0106          /// return no. of matrix columns
0107          kCols = D2,
0108          /// return no of elements: rows*columns
0109          kSize = D1*D2
0110       };
0111
0112    private:
0113       //T __attribute__ ((aligned (16))) fArray[kSize];
0114       T  fArray[kSize];
0115    };
0116
0117
0118 //     template<unsigned int D>
0119 //     struct Creator {
0120 //       static const RowOffsets<D> & Offsets() {
0121 //          static RowOffsets<D> off;
0122 //           return off;
0123 //       }
0124
0125    /**
0126       Static structure to keep the conversion from (i,j) to offsets in the storage data for a
0127       symmetric matrix
0128    */
0129
0130    template<unsigned int D>
0131    struct RowOffsets {
0132       inline RowOffsets() {
0133          int v[D];
0134          v[0]=0;
0135          for (unsigned int i=1; i<D; ++i)
0136             v[i]=v[i-1]+i;
0137          for (unsigned int i=0; i<D; ++i) {
0138             for (unsigned int j=0; j<=i; ++j)
0139                fOff[i*D+j] = v[i]+j;
0140             for (unsigned int j=i+1; j<D; ++j)
0141                fOff[i*D+j] = v[j]+i ;
0142          }
0143       }
0144       inline int operator()(unsigned int i, unsigned int j) const { return fOff[i*D+j]; }
0145       inline int apply(unsigned int i) const { return fOff[i]; }
0146       int fOff[D*D];
0147    };
0148
0149   namespace rowOffsetsUtils {
0150
0151     ///////////
0152     // Some meta template stuff
0153     template<int...> struct indices{};
0154
0155     template<int I, class IndexTuple, int N>
0156     struct make_indices_impl;
0157
0158     template<int I, int... Indices, int N>
0159     struct make_indices_impl<I, indices<Indices...>, N>
0160     {
0161       typedef typename make_indices_impl<I + 1, indices<Indices..., I>,
0162                      N>::type type;
0163     };
0164
0165     template<int N, int... Indices>
0166     struct make_indices_impl<N, indices<Indices...>, N> {
0167       typedef indices<Indices...> type;
0168     };
0169
0170     template<int N>
0171     struct make_indices : make_indices_impl<0, indices<>, N> {};
0172     // end of stuff
0173
0174
0175
0176     template<int I0, class F, int... I>
0177     constexpr std::array<decltype(std::declval<F>()(std::declval<int>())), sizeof...(I)>
0178     do_make(F f, indices<I...>)
0179     {
0180       return  std::array<decltype(std::declval<F>()(std::declval<int>())),
0181              sizeof...(I)>{{ f(I0 + I)... }};
0182     }
0183
0184     template<int N, int I0 = 0, class F>
0185     constexpr std::array<decltype(std::declval<F>()(std::declval<int>())), N>
0186     make(F f) {
0187       return do_make<I0>(f, typename make_indices<N>::type());
0188     }
0189
0190   } // namespace rowOffsetsUtils
0191
0192
0193 //_________________________________________________________________________________
0194    /**
0195       MatRepSym
0196       Matrix storage representation for a symmetric matrix of dimension NxN
0197       This class is a template on the contained type and on the symmetric matrix size, N.
0198       It has as data member an array of type T of size N*(N+1)/2,
0199       containing the lower diagonal block of the matrix.
0200       The order follows the lower diagonal block, still in a row-major convention.
0201       For example for a symmetric 3x3 matrix the order of the 6 elements
0202       \f$ \left[a_0,a_1.....a_5 \right]\f$ is:
0203       \f[
0204       M = \left( \begin{array}{ccc}
0205       a_0 & a_1  & a_3  \\
0206       a_1 & a_2  & a_4  \\
0207       a_3 & a_4 & a_5   \end{array} \right)
0208       \f]
0209
0210       @ingroup MatRep
0211    */
0212    template <class T, unsigned int D>
0213    class MatRepSym {
0214
0215    public:
0216
0217     /* constexpr */ inline MatRepSym(){}
0218
0219     typedef T  value_type;
0220
0221
0222     inline T & operator()(unsigned int i, unsigned int j)
0223      { return fArray[offset(i, j)]; }
0224
0225      inline /* constexpr */ T const & operator()(unsigned int i, unsigned int j) const
0226      { return fArray[offset(i, j)]; }
0227
0228      inline T& operator[](unsigned int i) {
0229        return fArray[off(i)];
0230      }
0231
0232      inline /* constexpr */ T const & operator[](unsigned int i) const {
0233        return fArray[off(i)];
0234      }
0235
0236      inline /* constexpr */ T apply(unsigned int i) const {
0237        return fArray[off(i)];
0238      }
0239
0240      inline T* Array() { return fArray; }
0241
0242      inline const T* Array() const { return fArray; }
0243
0244       /**
0245          assignment : only symmetric to symmetric allowed
0246        */
0247       template <class R>
0248       inline MatRepSym<T, D>& operator=(const R&) {
0249          STATIC_CHECK(0==1,
0250                       Cannot_assign_general_to_symmetric_matrix_representation);
0251          return *this;
0252       }
0253       inline MatRepSym<T, D>& operator=(const MatRepSym& rhs) {
0254          for(unsigned int i=0; i<kSize; ++i) fArray[i] = rhs.Array()[i];
0255          return *this;
0256       }
0257
0258       /**
0259          self addition : only symmetric to symmetric allowed
0260        */
0261       template <class R>
0262       inline MatRepSym<T, D>& operator+=(const R&) {
0263          STATIC_CHECK(0==1,
0264                       Cannot_add_general_to_symmetric_matrix_representation);
0265          return *this;
0266       }
0267       inline MatRepSym<T, D>& operator+=(const MatRepSym& rhs) {
0268          for(unsigned int i=0; i<kSize; ++i) fArray[i] += rhs.Array()[i];
0269          return *this;
0270       }
0271
0272       /**
0273          self subtraction : only symmetric to symmetric allowed
0274        */
0275       template <class R>
0276       inline MatRepSym<T, D>& operator-=(const R&) {
0277          STATIC_CHECK(0==1,
0278                       Cannot_substract_general_to_symmetric_matrix_representation);
0279          return *this;
0280       }
0281       inline MatRepSym<T, D>& operator-=(const MatRepSym& rhs) {
0282          for(unsigned int i=0; i<kSize; ++i) fArray[i] -= rhs.Array()[i];
0283          return *this;
0284       }
0285       template <class R>
0286       inline bool operator==(const R& rhs) const {
0287          bool rc = true;
0288          for(unsigned int i=0; i<D*D; ++i) {
0289             rc = rc && (operator[](i) == rhs[i]);
0290          }
0291          return rc;
0292       }
0293
0294       enum {
0295          /// return no. of matrix rows
0296          kRows = D,
0297          /// return no. of matrix columns
0298          kCols = D,
0299          /// return no of elements: rows*columns
0300          kSize = D*(D+1)/2
0301       };
0302
0303      static constexpr int off0(int i) { return i==0 ? 0 : off0(i-1)+i;}
0304      static constexpr int off2(int i, int j) { return j<i ? off0(i)+j : off0(j)+i; }
0305      static constexpr int off1(int i) { return off2(i/D, i%D);}
0306
0307      static int off(int i) {
0308        static constexpr auto v = rowOffsetsUtils::make<D*D>(off1);
0309        return v[i];
0310      }
0311
0312      static inline constexpr unsigned int
0313      offset(unsigned int i, unsigned int j)
0314      {
0315        //if (j > i) std::swap(i, j);
0316        return off(i*D+j);
0317        // return (i>j) ? (i * (i+1) / 2) + j :  (j * (j+1) / 2) + i;
0318      }
0319
0320    private:
0321       //T __attribute__ ((aligned (16))) fArray[kSize];
0322       T fArray[kSize];
0323    };
0324
0325
0326
0327 } // namespace Math
0328 } // namespace ROOT
0329
0330
0331 #endif // MATH_MATRIXREPRESENTATIONSSTATIC_H