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