Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
LAVector.h
Go to the documentation of this file.
1// @(#)root/minuit2:$Id$
2// Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3
4/**********************************************************************
5 * *
6 * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7 * *
8 **********************************************************************/
9
10#ifndef ROOT_Minuit2_LAVector
11#define ROOT_Minuit2_LAVector
12
13#include "Minuit2/ABSum.h"
14#include "Minuit2/ABProd.h"
15#include "Minuit2/LASymMatrix.h"
16
17#include <cassert>
18#include <memory>
19
21
22namespace ROOT {
23
24namespace Minuit2 {
25
26// extern StackAllocator StackAllocatorHolder::Get();
27
28int Mndaxpy(unsigned int, double, const double *, int, double *, int);
29int Mndscal(unsigned int, double, double *, int);
30int Mndspmv(const char *, unsigned int, double, const double *, const double *, int, double, double *, int);
31
32class LAVector {
33
34private:
35 LAVector() : fSize(0), fData(nullptr) {}
36
37public:
38 typedef vec Type;
39
40 LAVector(unsigned int n) :
41 fSize(n),
42 fData( (n>0) ? (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * n)
43 : nullptr)
44 {
45 // assert(fSize>0);
46 if (fData)
47 std::memset(fData, 0, size() * sizeof(double));
48 // std::cout<<"LAVector(unsigned int n), n= "<<n<<std::endl;
49 }
50
52 {
53 // std::cout<<"~LAVector()"<<std::endl;
54 // if(fData) std::cout<<"deleting "<<fSize<<std::endl;
55 // else std::cout<<"no delete"<<std::endl;
56 // if(fData) delete [] fData;
57 if (fData)
59 }
60
62 : fSize(v.size()), fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * v.size()))
63 {
64 // std::cout<<"LAVector(const LAVector& v)"<<std::endl;
65 std::memcpy(fData, v.Data(), fSize * sizeof(double));
66 }
67
69 {
70 // implement proper copy constructor in case size is different
71 if (v.size() > fSize) {
73 fSize = v.size();
74 fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
75 }
76 std::memcpy(fData, v.Data(), fSize * sizeof(double));
77 return *this;
78 }
79
80 template <class T>
82 : fSize(v.Obj().size()), fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * v.Obj().size()))
83 {
84 // std::cout<<"LAVector(const ABObj<LAVector, T>& v)"<<std::endl;
85 // std::cout<<"allocate "<<fSize<<std::endl;
86 std::memcpy(fData, v.Obj().Data(), fSize * sizeof(T));
87 (*this) *= T(v.f());
88 // std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
89 }
90
91 template <class A, class B, class T>
93 {
94 // std::cout<<"template<class A, class B, class T> LAVector(const ABObj<ABSum<ABObj<A, T>, ABObj<B, T> > >&
95 // sum)"<<std::endl;
96 (*this) = sum.Obj().A();
97 (*this) += sum.Obj().B();
98 (*this) *= double(sum.f());
99 }
100
101 template <class A, class T>
103 {
104 // std::cout<<"template<class A, class T> LAVector(const ABObj<ABSum<ABObj<LAVector, T>, ABObj<A, T> >,T>&
105 // sum)"<<std::endl;
106
107 // recursive construction
108 // std::cout<<"(*this)=sum.Obj().B();"<<std::endl;
109 (*this) = sum.Obj().B();
110 // std::cout<<"(*this)+=sum.Obj().A();"<<std::endl;
111 (*this) += sum.Obj().A();
112 (*this) *= double(sum.f());
113 // std::cout<<"leaving template<class A, class T> LAVector(const ABObj<ABSum<ABObj<LAVector,.."<<std::endl;
114 }
115
116 template <class A, class T>
117 LAVector(const ABObj<vec, ABObj<vec, A, T>, T> &something) : fSize(0), fData(nullptr)
118 {
119 // std::cout<<"template<class A, class T> LAVector(const ABObj<ABObj<A, T>, T>& something)"<<std::endl;
120 (*this) = something.Obj();
121 (*this) *= something.f();
122 }
123
124 //
125 template <class T>
127 : fSize(prod.Obj().B().Obj().size()),
128 fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * prod.Obj().B().Obj().size()))
129 {
130 // std::cout<<"template<class T> LAVector(const ABObj<vec, ABProd<ABObj<sym, LASymMatrix, T>, ABObj<vec,
131 // LAVector, T> >, T>& prod)"<<std::endl;
132
133 Mndspmv("U", fSize, prod.f() * prod.Obj().A().f() * prod.Obj().B().f(), prod.Obj().A().Obj().Data(),
134 prod.Obj().B().Obj().Data(), 1, 0., fData, 1);
135 }
136
137 //
138 template <class T>
140 vec,
142 T> &prod)
143 : fSize(0), fData(nullptr)
144 {
145 (*this) = prod.Obj().B();
146 (*this) += prod.Obj().A();
147 (*this) *= double(prod.f());
148 }
149
150 //
152 {
153 // std::cout<<"LAVector& operator+=(const LAVector& m)"<<std::endl;
154 assert(fSize == m.size());
155 Mndaxpy(fSize, 1., m.Data(), 1, fData, 1);
156 return *this;
157 }
158
160 {
161 // std::cout<<"LAVector& operator-=(const LAVector& m)"<<std::endl;
162 assert(fSize == m.size());
163 Mndaxpy(fSize, -1., m.Data(), 1, fData, 1);
164 return *this;
165 }
166
167 template <class T>
169 {
170 // std::cout<<"template<class T> LAVector& operator+=(const ABObj<LAVector, T>& m)"<<std::endl;
171 assert(fSize == m.Obj().size());
172 if (m.Obj().Data() == fData) {
173 Mndscal(fSize, 1. + double(m.f()), fData, 1);
174 } else {
175 Mndaxpy(fSize, double(m.f()), m.Obj().Data(), 1, fData, 1);
176 }
177 // std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
178 return *this;
179 }
180
181 template <class A, class T>
183 {
184 // std::cout<<"template<class A, class T> LAVector& operator+=(const ABObj<A,T>& m)"<<std::endl;
185 (*this) += LAVector(m);
186 return *this;
187 }
188
189 template <class T>
191 {
192 Mndspmv("U", fSize, prod.f() * prod.Obj().A().f() * prod.Obj().B().f(), prod.Obj().A().Obj().Data(),
193 prod.Obj().B().Data(), 1, 1., fData, 1);
194 return *this;
195 }
196
197 LAVector &operator*=(double scal)
198 {
199 Mndscal(fSize, scal, fData, 1);
200 return *this;
201 }
202
203 double operator()(unsigned int i) const
204 {
205 assert(i < fSize);
206 return fData[i];
207 }
208
209 double &operator()(unsigned int i)
210 {
211 assert(i < fSize);
212 return fData[i];
213 }
214
215 double operator[](unsigned int i) const
216 {
217 assert(i < fSize);
218 return fData[i];
219 }
220
221 double &operator[](unsigned int i)
222 {
223 assert(i < fSize);
224 return fData[i];
225 }
226
227 const double *Data() const { return fData; }
228
229 double *Data() { return fData; }
230
231 unsigned int size() const { return fSize; }
232
233private:
234 unsigned int fSize;
235 double *fData;
236
237public:
238 template <class T>
240 {
241 // std::cout<<"template<class T> LAVector& operator=(ABObj<LAVector, T>& v)"<<std::endl;
242 if (fSize == 0 && !fData) {
243 fSize = v.Obj().size();
244 fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
245 } else {
246 assert(fSize == v.Obj().size());
247 }
248 std::memcpy(fData, v.Obj().Data(), fSize * sizeof(double));
249 (*this) *= T(v.f());
250 return *this;
251 }
252
253 template <class A, class T>
255 {
256 // std::cout<<"template<class A, class T> LAVector& operator=(const ABObj<ABObj<A, T>, T>&
257 // something)"<<std::endl;
258 if (fSize == 0 && !fData) {
259 (*this) = something.Obj();
260 } else {
261 LAVector tmp(something.Obj());
262 assert(fSize == tmp.size());
263 std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
264 }
265 (*this) *= something.f();
266 return *this;
267 }
268
269 template <class A, class B, class T>
271 {
272 if (fSize == 0 && !fData) {
273 (*this) = sum.Obj().A();
274 (*this) += sum.Obj().B();
275 } else {
276 LAVector tmp(sum.Obj().A());
277 tmp += sum.Obj().B();
278 assert(fSize == tmp.size());
279 std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
280 }
281 (*this) *= sum.f();
282 return *this;
283 }
284
285 template <class A, class T>
287 {
288 if (fSize == 0 && !fData) {
289 (*this) = sum.Obj().B();
290 (*this) += sum.Obj().A();
291 } else {
292 LAVector tmp(sum.Obj().A());
293 tmp += sum.Obj().B();
294 assert(fSize == tmp.size());
295 std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
296 }
297 (*this) *= sum.f();
298 return *this;
299 }
300
301 //
302 template <class T>
304 {
305 if (fSize == 0 && !fData) {
306 fSize = prod.Obj().B().Obj().size();
307 fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
308 Mndspmv("U", fSize, double(prod.f() * prod.Obj().A().f() * prod.Obj().B().f()), prod.Obj().A().Obj().Data(),
309 prod.Obj().B().Obj().Data(), 1, 0., fData, 1);
310 } else {
311 LAVector tmp(prod.Obj().B());
312 assert(fSize == tmp.size());
313 Mndspmv("U", fSize, double(prod.f() * prod.Obj().A().f()), prod.Obj().A().Obj().Data(), tmp.Data(), 1, 0.,
314 fData, 1);
315 }
316 return *this;
317 }
318
319 //
320 template <class T>
321 LAVector &
323 vec,
325 T> &prod)
326 {
327 if (fSize == 0 && !fData) {
328 (*this) = prod.Obj().B();
329 (*this) += prod.Obj().A();
330 } else {
331 // std::cout<<"creating tmp variable"<<std::endl;
332 LAVector tmp(prod.Obj().B());
333 tmp += prod.Obj().A();
334 assert(fSize == tmp.size());
335 std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
336 }
337 (*this) *= prod.f();
338 return *this;
339 }
340};
341
342} // namespace Minuit2
343
344} // namespace ROOT
345
346#endif // ROOT_Minuit2_LAVector
LAVector & operator=(const ABObj< vec, ABSum< ABObj< vec, A, T >, ABObj< vec, B, T > >, T > &sum)
Definition LAVector.h:270
const double * Data() const
Definition LAVector.h:227
LAVector & operator*=(double scal)
Definition LAVector.h:197
unsigned int size() const
Definition LAVector.h:231
double operator[](unsigned int i) const
Definition LAVector.h:215
LAVector(unsigned int n)
Definition LAVector.h:40
double operator()(unsigned int i) const
Definition LAVector.h:203
LAVector & operator+=(const LAVector &m)
Definition LAVector.h:151
LAVector(const ABObj< vec, ABProd< ABObj< sym, LASymMatrix, T >, ABObj< vec, LAVector, T > >, T > &prod)
Definition LAVector.h:126
LAVector(const ABObj< vec, ABObj< vec, A, T >, T > &something)
Definition LAVector.h:117
LAVector & operator=(const ABObj< vec, LAVector, T > &v)
Definition LAVector.h:239
LAVector(const ABObj< vec, ABSum< ABObj< vec, A, T >, ABObj< vec, B, T > >, T > &sum)
Definition LAVector.h:92
double & operator[](unsigned int i)
Definition LAVector.h:221
LAVector & operator=(const ABObj< vec, ABObj< vec, A, T >, T > &something)
Definition LAVector.h:254
LAVector & operator=(const LAVector &v)
Definition LAVector.h:68
LAVector(const ABObj< vec, ABSum< ABObj< vec, LAVector, T >, ABObj< vec, A, T > >, T > &sum)
Definition LAVector.h:102
LAVector & operator+=(const ABObj< vec, A, T > &m)
Definition LAVector.h:182
LAVector(const ABObj< vec, ABSum< ABObj< vec, ABProd< ABObj< sym, LASymMatrix, T >, ABObj< vec, LAVector, T > >, T >, ABObj< vec, LAVector, T > >, T > &prod)
Definition LAVector.h:139
LAVector(const ABObj< vec, LAVector, T > &v)
Definition LAVector.h:81
LAVector & operator-=(const LAVector &m)
Definition LAVector.h:159
LAVector(const LAVector &v)
Definition LAVector.h:61
LAVector & operator=(const ABObj< vec, ABSum< ABObj< vec, ABProd< ABObj< sym, LASymMatrix, T >, ABObj< vec, LAVector, T > >, T >, ABObj< vec, LAVector, T > >, T > &prod)
Definition LAVector.h:322
double & operator()(unsigned int i)
Definition LAVector.h:209
LAVector & operator=(const ABObj< vec, ABSum< ABObj< vec, LAVector, T >, ABObj< vec, A, T > >, T > &sum)
Definition LAVector.h:286
LAVector & operator=(const ABObj< vec, ABProd< ABObj< sym, LASymMatrix, T >, ABObj< vec, LAVector, T > >, T > &prod)
Definition LAVector.h:303
LAVector & operator+=(const ABObj< vec, ABProd< ABObj< sym, LASymMatrix, T >, ABObj< vec, LAVector, T > >, T > &prod)
Definition LAVector.h:190
LAVector & operator+=(const ABObj< vec, LAVector, T > &m)
Definition LAVector.h:168
static StackAllocator & Get()
void * Allocate(size_t nBytes)
const Int_t n
Definition legend1.C:16
int Mndaxpy(unsigned int, double, const double *, int, double *, int)
Definition mndaxpy.cxx:19
int Mndscal(unsigned int, double, double *, int)
Definition mndscal.cxx:19
int Mndspmv(const char *, unsigned int, double, const double *, const double *, int, double, double *, int)
Definition mndspmv.cxx:22
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
TMarker m
Definition textangle.C:8
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345