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
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
29namespace ROOT {
30
31namespace 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>
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 }
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 }
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
#define f(i)
Definition: RSha256.hxx:104
#define STATIC_CHECK(expr, msg)
Definition: StaticCheck.h:56
#define N
int type
Definition: TGX11.cxx:120
Expression wrapper class for Matrix objects.
bool operator==(const R &rhs) const
@ kCols
return no. of matrix columns
@ kSize
return no of elements: rows*columns
@ kRows
return no. of matrix rows
const T & operator()(unsigned int i, unsigned int j) const
T apply(unsigned int i) const
MatRepStd< T, D1, D2 > & operator-=(const R &rhs)
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)
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)
@ kCols
return no. of matrix columns
@ kRows
return no. of matrix rows
@ kSize
return no of elements: rows*columns
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...
Definition: StringConv.hxx:21
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