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