ROOT   6.10/09 Reference Guide
MatrixRepresentationsStatic.h
Go to the documentation of this file.
1 // @(#)root/smatrix:$Id$
2 // Author: L. Moneta, J. Palacios 2006
3
4 #ifndef ROOT_Math_MatrixRepresentationsStatic
5 #define ROOT_Math_MatrixRepresentationsStatic 1
6
7 // Include files
8
9 /**
10  @defgroup MatRep SMatrix Storage Representation
11  @ingroup SMatrixGroup
12
13  @author Juan Palacios
14  @date 2006-01-15
15
16  Classes MatRepStd and MatRepSym for generic and symmetric matrix
17  data storage and manipulation. Define data storage and access, plus
18  operators =, +=, -=, ==.
19
20  */
21
22 #include "Math/StaticCheck.h"
23
24 #include <cstddef>
25 #include <utility>
26 #include <type_traits>
27 #include <array>
28
29 namespace ROOT {
30
31 namespace Math {
32
33  //________________________________________________________________________________
34  /**
35  MatRepStd
36  Standard Matrix representation for a general D1 x D2 matrix.
37  This class is itself a template on the contained type T, the number of rows and the number of columns.
38  Its data member is an array T[nrows*ncols] containing the matrix data.
39  The data are stored in the row-major C convention.
40  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 are stored in the following order:
41  \f[
42  M = \left( \begin{array}{ccc}
43  a_0 & a_1 & a_2 \\
44  a_3 & a_4 & a_5 \\
45  a_6 & a_7 & a_8 \end{array} \right)
46  \f]
47
48  @ingroup MatRep
49  */
50
51
52  template <class T, unsigned int D1, unsigned int D2=D1>
53  class MatRepStd {
54
55  public:
56
57  typedef T value_type;
58
59  inline const T& operator()(unsigned int i, unsigned int j) const {
60  return fArray[i*D2+j];
61  }
62  inline T& operator()(unsigned int i, unsigned int j) {
63  return fArray[i*D2+j];
64  }
65  inline T& operator[](unsigned int i) { return fArray[i]; }
66
67  inline const T& operator[](unsigned int i) const { return fArray[i]; }
68
69  inline T apply(unsigned int i) const { return fArray[i]; }
70
71  inline T* Array() { return fArray; }
72
73  inline const T* Array() const { return fArray; }
74
75  template <class R>
76  inline MatRepStd<T, D1, D2>& operator+=(const R& rhs) {
77  for(unsigned int i=0; i<kSize; ++i) fArray[i] += rhs[i];
78  return *this;
79  }
80
81  template <class R>
82  inline MatRepStd<T, D1, D2>& operator-=(const R& rhs) {
83  for(unsigned int i=0; i<kSize; ++i) fArray[i] -= rhs[i];
84  return *this;
85  }
86
87  template <class R>
88  inline MatRepStd<T, D1, D2>& operator=(const R& rhs) {
89  for(unsigned int i=0; i<kSize; ++i) fArray[i] = rhs[i];
90  return *this;
91  }
92
93  template <class R>
94  inline bool operator==(const R& rhs) const {
95  bool rc = true;
96  for(unsigned int i=0; i<kSize; ++i) {
97  rc = rc && (fArray[i] == rhs[i]);
98  }
99  return rc;
100  }
101
102  enum {
103  /// return no. of matrix rows
104  kRows = D1,
105  /// return no. of matrix columns
106  kCols = D2,
107  /// return no of elements: rows*columns
108  kSize = D1*D2
109  };
110
111  private:
112  //T __attribute__ ((aligned (16))) fArray[kSize];
114  };
115
116
117 // template<unigned int D>
118 // struct Creator {
119 // static const RowOffsets<D> & Offsets() {
120 // static RowOffsets<D> off;
121 // return off;
122 // }
123
124  /**
125  Static structure to keep the conversion from (i,j) to offsets in the storage data for a
126  symmetric matrix
127  */
128
129  template<unsigned int D>
130  struct RowOffsets {
131  inline RowOffsets() {
132  int v[D];
133  v[0]=0;
134  for (unsigned int i=1; i<D; ++i)
135  v[i]=v[i-1]+i;
136  for (unsigned int i=0; i<D; ++i) {
137  for (unsigned int j=0; j<=i; ++j)
138  fOff[i*D+j] = v[i]+j;
139  for (unsigned int j=i+1; j<D; ++j)
140  fOff[i*D+j] = v[j]+i ;
141  }
142  }
143  inline int operator()(unsigned int i, unsigned int j) const { return fOff[i*D+j]; }
144  inline int apply(unsigned int i) const { return fOff[i]; }
145  int fOff[D*D];
146  };
147
148  namespace rowOffsetsUtils {
149
150  ///////////
151  // Some meta template stuff
152  template<int...> struct indices{};
153
154  template<int I, class IndexTuple, int N>
156
157  template<int I, int... Indices, int N>
158  struct make_indices_impl<I, indices<Indices...>, N>
159  {
160  typedef typename make_indices_impl<I + 1, indices<Indices..., I>,
162  };
163
164  template<int N, int... Indices>
165  struct make_indices_impl<N, indices<Indices...>, N> {
166  typedef indices<Indices...> type;
167  };
168
169  template<int N>
170  struct make_indices : make_indices_impl<0, indices<>, N> {};
171  // end of stuff
172
173
174
175  template<int I0, class F, int... I>
176  constexpr std::array<decltype(std::declval<F>()(std::declval<int>())), sizeof...(I)>
178  {
179  return std::array<decltype(std::declval<F>()(std::declval<int>())),
180  sizeof...(I)>{{ f(I0 + I)... }};
181  }
182
183  template<int N, int I0 = 0, class F>
184  constexpr std::array<decltype(std::declval<F>()(std::declval<int>())), N>
185  make(F f) {
186  return do_make<I0>(f, typename make_indices<N>::type());
187  }
188
189  } // namespace rowOffsetsUtils
190
191
192 //_________________________________________________________________________________
193  /**
194  MatRepSym
195  Matrix storage representation for a symmetric matrix of dimension NxN
196  This class is a template on the contained type and on the symmetric matrix size, N.
197  It has as data member an array of type T of size N*(N+1)/2,
198  containing the lower diagonal block of the matrix.
199  The order follows the lower diagonal block, still in a row-major convention.
200  For example for a symmetric 3x3 matrix the order of the 6 elements
201  \f$\left[a_0,a_1.....a_5 \right]\f$ is:
202  \f[
203  M = \left( \begin{array}{ccc}
204  a_0 & a_1 & a_3 \\
205  a_1 & a_2 & a_4 \\
206  a_3 & a_4 & a_5 \end{array} \right)
207  \f]
208
209  @ingroup MatRep
210  */
211  template <class T, unsigned int D>
212  class MatRepSym {
213
214  public:
215
216  /* constexpr */ inline MatRepSym(){}
217
218  typedef T value_type;
219
220
221  inline T & operator()(unsigned int i, unsigned int j)
222  { return fArray[offset(i, j)]; }
223
224  inline /* constexpr */ T const & operator()(unsigned int i, unsigned int j) const
225  { return fArray[offset(i, j)]; }
226
227  inline T& operator[](unsigned int i) {
228  return fArray[off(i)];
229  }
230
231  inline /* constexpr */ T const & operator[](unsigned int i) const {
232  return fArray[off(i)];
233  }
234
235  inline /* constexpr */ T apply(unsigned int i) const {
236  return fArray[off(i)];
237  }
238
239  inline T* Array() { return fArray; }
240
241  inline const T* Array() const { return fArray; }
242
243  /**
244  assignment : only symmetric to symmetric allowed
245  */
246  template <class R>
247  inline MatRepSym<T, D>& operator=(const R&) {
248  STATIC_CHECK(0==1,
249  Cannot_assign_general_to_symmetric_matrix_representation);
250  return *this;
251  }
252  inline MatRepSym<T, D>& operator=(const MatRepSym& rhs) {
253  for(unsigned int i=0; i<kSize; ++i) fArray[i] = rhs.Array()[i];
254  return *this;
255  }
256
257  /**
258  self addition : only symmetric to symmetric allowed
259  */
260  template <class R>
261  inline MatRepSym<T, D>& operator+=(const R&) {
262  STATIC_CHECK(0==1,
264  return *this;
265  }
266  inline MatRepSym<T, D>& operator+=(const MatRepSym& rhs) {
267  for(unsigned int i=0; i<kSize; ++i) fArray[i] += rhs.Array()[i];
268  return *this;
269  }
270
271  /**
272  self subtraction : only symmetric to symmetric allowed
273  */
274  template <class R>
275  inline MatRepSym<T, D>& operator-=(const R&) {
276  STATIC_CHECK(0==1,
277  Cannot_substract_general_to_symmetric_matrix_representation);
278  return *this;
279  }
280  inline MatRepSym<T, D>& operator-=(const MatRepSym& rhs) {
281  for(unsigned int i=0; i<kSize; ++i) fArray[i] -= rhs.Array()[i];
282  return *this;
283  }
284  template <class R>
285  inline bool operator==(const R& rhs) const {
286  bool rc = true;
287  for(unsigned int i=0; i<D*D; ++i) {
288  rc = rc && (operator[](i) == rhs[i]);
289  }
290  return rc;
291  }
292
293  enum {
294  /// return no. of matrix rows
295  kRows = D,
296  /// return no. of matrix columns
297  kCols = D,
298  /// return no of elements: rows*columns
299  kSize = D*(D+1)/2
300  };
301
302  static constexpr int off0(int i) { return i==0 ? 0 : off0(i-1)+i;}
303  static constexpr int off2(int i, int j) { return j<i ? off0(i)+j : off0(j)+i; }
304  static constexpr int off1(int i) { return off2(i/D, i%D);}
305
306  static int off(int i) {
307  static constexpr auto v = rowOffsetsUtils::make<D*D>(off1);
308  return v[i];
309  }
310
311  static inline constexpr unsigned int
312  offset(unsigned int i, unsigned int j)
313  {
314  //if (j > i) std::swap(i, j);
315  return off(i*D+j);
316  // return (i>j) ? (i * (i+1) / 2) + j : (j * (j+1) / 2) + i;
317  }
318
319  private:
320  //T __attribute__ ((aligned (16))) fArray[kSize];
322  };
323
324
325
326 } // namespace Math
327 } // namespace ROOT
328
329
330 #endif // MATH_MATRIXREPRESENTATIONSSTATIC_H
int apply(unsigned int i) const
T & operator()(unsigned int i, unsigned int j)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
MatRepStd< T, D1, D2 > & operator=(const R &rhs)
static constexpr int off2(int i, int j)
T & operator()(unsigned int i, unsigned int j)
double T(double x)
Definition: ChebyshevPol.h:34
MatRepSym< T, D > & operator-=(const R &)
self subtraction : only symmetric to symmetric allowed
MatRepSym Matrix storage representation for a symmetric matrix of dimension NxN This class is a templ...
Definition: HelperOps.h:35
MatRepSym< T, D > & operator=(const R &)
assignment : only symmetric to symmetric allowed
#define N
bool operator==(const R &rhs) const
constexpr std::array< decltype(std::declval< F >)(std::declval< int >))), N > make(F f)
return no. of matrix columns
#define STATIC_CHECK(expr, msg)
Definition: StaticCheck.h:56
Static structure to keep the conversion from (i,j) to offsets in the storage data for a symmetric mat...
MatRepSym< T, D > & operator+=(const MatRepSym &rhs)
T const & operator[](unsigned int i) const
int operator()(unsigned int i, unsigned int j) const
Expression wrapper class for Matrix objects.
Definition: Expression.h:134
T apply(unsigned int i) const
static constexpr int off0(int i)
#define F(x, y, z)
SVector< double, 2 > v
Definition: Dict.h:5
T apply(unsigned int i) const
MatRepSym< T, D > & operator=(const MatRepSym &rhs)
return no of elements: rows*columns
MatRepStd< T, D1, D2 > & operator+=(const R &rhs)
double f(double x)
int type
Definition: TGX11.cxx:120
static constexpr unsigned int offset(unsigned int i, unsigned int j)
Namespace for new Math classes and functions.
MatRepSym< T, D > & operator+=(const R &)
self addition : only symmetric to symmetric allowed
T const & operator()(unsigned int i, unsigned int j) const
const T & operator()(unsigned int i, unsigned int j) const
constexpr std::array< decltype(std::declval< F >)(std::declval< int >))), sizeof...(I)> do_make(F f, indices< I... >)
static constexpr int off1(int i)
const T & operator[](unsigned int i) const
#define I(x, y, z)
bool operator==(const R &rhs) const
MatRepSym< T, D > & operator-=(const MatRepSym &rhs)
TRandom3 R
a TMatrixD.
Definition: testIO.cxx:28
MatRepStd< T, D1, D2 > & operator-=(const R &rhs)