Logo ROOT  
Reference Guide
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(0) {}
36
37public:
38 typedef vec Type;
39
40 // LAVector() : fSize(0), fData(0) {}
41
42 LAVector(unsigned int n) : fSize(n), fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * n))
43 {
44 // assert(fSize>0);
45 std::memset(fData, 0, size() * sizeof(double));
46 // std::cout<<"LAVector(unsigned int n), n= "<<n<<std::endl;
47 }
48
50 {
51 // std::cout<<"~LAVector()"<<std::endl;
52 // if(fData) std::cout<<"deleting "<<fSize<<std::endl;
53 // else std::cout<<"no delete"<<std::endl;
54 // if(fData) delete [] fData;
55 if (fData)
57 }
58
60 : fSize(v.size()), fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * v.size()))
61 {
62 // std::cout<<"LAVector(const LAVector& v)"<<std::endl;
63 std::memcpy(fData, v.Data(), fSize * sizeof(double));
64 }
65
67 {
68 // std::cout<<"LAVector& operator=(const LAVector& v)"<<std::endl;
69 // std::cout<<"fSize= "<<fSize<<std::endl;
70 // std::cout<<"v.size()= "<<v.size()<<std::endl;
71 assert(fSize == v.size());
72 std::memcpy(fData, v.Data(), fSize * sizeof(double));
73 return *this;
74 }
75
76 template <class T>
78 : fSize(v.Obj().size()), fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * v.Obj().size()))
79 {
80 // std::cout<<"LAVector(const ABObj<LAVector, T>& v)"<<std::endl;
81 // std::cout<<"allocate "<<fSize<<std::endl;
82 std::memcpy(fData, v.Obj().Data(), fSize * sizeof(T));
83 (*this) *= T(v.f());
84 // std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
85 }
86
87 template <class A, class B, class T>
89 {
90 // std::cout<<"template<class A, class B, class T> LAVector(const ABObj<ABSum<ABObj<A, T>, ABObj<B, T> > >&
91 // sum)"<<std::endl;
92 (*this) = sum.Obj().A();
93 (*this) += sum.Obj().B();
94 (*this) *= double(sum.f());
95 }
96
97 template <class A, class T>
99 {
100 // std::cout<<"template<class A, class T> LAVector(const ABObj<ABSum<ABObj<LAVector, T>, ABObj<A, T> >,T>&
101 // sum)"<<std::endl;
102
103 // recursive construction
104 // std::cout<<"(*this)=sum.Obj().B();"<<std::endl;
105 (*this) = sum.Obj().B();
106 // std::cout<<"(*this)+=sum.Obj().A();"<<std::endl;
107 (*this) += sum.Obj().A();
108 (*this) *= double(sum.f());
109 // std::cout<<"leaving template<class A, class T> LAVector(const ABObj<ABSum<ABObj<LAVector,.."<<std::endl;
110 }
111
112 template <class A, class T>
113 LAVector(const ABObj<vec, ABObj<vec, A, T>, T> &something) : fSize(0), fData(0)
114 {
115 // std::cout<<"template<class A, class T> LAVector(const ABObj<ABObj<A, T>, T>& something)"<<std::endl;
116 (*this) = something.Obj();
117 (*this) *= something.f();
118 }
119
120 //
121 template <class T>
123 : fSize(prod.Obj().B().Obj().size()),
124 fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * prod.Obj().B().Obj().size()))
125 {
126 // std::cout<<"template<class T> LAVector(const ABObj<vec, ABProd<ABObj<sym, LASymMatrix, T>, ABObj<vec,
127 // LAVector, T> >, T>& prod)"<<std::endl;
128
129 Mndspmv("U", fSize, prod.f() * prod.Obj().A().f() * prod.Obj().B().f(), prod.Obj().A().Obj().Data(),
130 prod.Obj().B().Obj().Data(), 1, 0., fData, 1);
131 }
132
133 //
134 template <class T>
136 vec,
138 T> &prod)
139 : fSize(0), fData(0)
140 {
141 (*this) = prod.Obj().B();
142 (*this) += prod.Obj().A();
143 (*this) *= double(prod.f());
144 }
145
146 //
148 {
149 // std::cout<<"LAVector& operator+=(const LAVector& m)"<<std::endl;
150 assert(fSize == m.size());
151 Mndaxpy(fSize, 1., m.Data(), 1, fData, 1);
152 return *this;
153 }
154
156 {
157 // std::cout<<"LAVector& operator-=(const LAVector& m)"<<std::endl;
158 assert(fSize == m.size());
159 Mndaxpy(fSize, -1., m.Data(), 1, fData, 1);
160 return *this;
161 }
162
163 template <class T>
165 {
166 // std::cout<<"template<class T> LAVector& operator+=(const ABObj<LAVector, T>& m)"<<std::endl;
167 assert(fSize == m.Obj().size());
168 if (m.Obj().Data() == fData) {
169 Mndscal(fSize, 1. + double(m.f()), fData, 1);
170 } else {
171 Mndaxpy(fSize, double(m.f()), m.Obj().Data(), 1, fData, 1);
172 }
173 // std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
174 return *this;
175 }
176
177 template <class A, class T>
179 {
180 // std::cout<<"template<class A, class T> LAVector& operator+=(const ABObj<A,T>& m)"<<std::endl;
181 (*this) += LAVector(m);
182 return *this;
183 }
184
185 template <class T>
187 {
188 Mndspmv("U", fSize, prod.f() * prod.Obj().A().f() * prod.Obj().B().f(), prod.Obj().A().Obj().Data(),
189 prod.Obj().B().Data(), 1, 1., fData, 1);
190 return *this;
191 }
192
193 LAVector &operator*=(double scal)
194 {
195 Mndscal(fSize, scal, fData, 1);
196 return *this;
197 }
198
199 double operator()(unsigned int i) const
200 {
201 assert(i < fSize);
202 return fData[i];
203 }
204
205 double &operator()(unsigned int i)
206 {
207 assert(i < fSize);
208 return fData[i];
209 }
210
211 double operator[](unsigned int i) const
212 {
213 assert(i < fSize);
214 return fData[i];
215 }
216
217 double &operator[](unsigned int i)
218 {
219 assert(i < fSize);
220 return fData[i];
221 }
222
223 const double *Data() const { return fData; }
224
225 double *Data() { return fData; }
226
227 unsigned int size() const { return fSize; }
228
229private:
230 unsigned int fSize;
231 double *fData;
232
233public:
234 template <class T>
236 {
237 // std::cout<<"template<class T> LAVector& operator=(ABObj<LAVector, T>& v)"<<std::endl;
238 if (fSize == 0 && fData == 0) {
239 fSize = v.Obj().size();
240 fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
241 } else {
242 assert(fSize == v.Obj().size());
243 }
244 std::memcpy(fData, v.Obj().Data(), fSize * sizeof(double));
245 (*this) *= T(v.f());
246 return *this;
247 }
248
249 template <class A, class T>
251 {
252 // std::cout<<"template<class A, class T> LAVector& operator=(const ABObj<ABObj<A, T>, T>&
253 // something)"<<std::endl;
254 if (fSize == 0 && fData == 0) {
255 (*this) = something.Obj();
256 } else {
257 LAVector tmp(something.Obj());
258 assert(fSize == tmp.size());
259 std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
260 }
261 (*this) *= something.f();
262 return *this;
263 }
264
265 template <class A, class B, class T>
267 {
268 if (fSize == 0 && fData == 0) {
269 (*this) = sum.Obj().A();
270 (*this) += sum.Obj().B();
271 } else {
272 LAVector tmp(sum.Obj().A());
273 tmp += sum.Obj().B();
274 assert(fSize == tmp.size());
275 std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
276 }
277 (*this) *= sum.f();
278 return *this;
279 }
280
281 template <class A, class T>
283 {
284 if (fSize == 0 && fData == 0) {
285 (*this) = sum.Obj().B();
286 (*this) += sum.Obj().A();
287 } else {
288 LAVector tmp(sum.Obj().A());
289 tmp += sum.Obj().B();
290 assert(fSize == tmp.size());
291 std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
292 }
293 (*this) *= sum.f();
294 return *this;
295 }
296
297 //
298 template <class T>
300 {
301 if (fSize == 0 && fData == 0) {
302 fSize = prod.Obj().B().Obj().size();
303 fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
304 Mndspmv("U", fSize, double(prod.f() * prod.Obj().A().f() * prod.Obj().B().f()), prod.Obj().A().Obj().Data(),
305 prod.Obj().B().Obj().Data(), 1, 0., fData, 1);
306 } else {
307 LAVector tmp(prod.Obj().B());
308 assert(fSize == tmp.size());
309 Mndspmv("U", fSize, double(prod.f() * prod.Obj().A().f()), prod.Obj().A().Obj().Data(), tmp.Data(), 1, 0.,
310 fData, 1);
311 }
312 return *this;
313 }
314
315 //
316 template <class T>
317 LAVector &
319 vec,
321 T> &prod)
322 {
323 if (fSize == 0 && fData == 0) {
324 (*this) = prod.Obj().B();
325 (*this) += prod.Obj().A();
326 } else {
327 // std::cout<<"creating tmp variable"<<std::endl;
328 LAVector tmp(prod.Obj().B());
329 tmp += prod.Obj().A();
330 assert(fSize == tmp.size());
331 std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
332 }
333 (*this) *= prod.f();
334 return *this;
335 }
336};
337
338} // namespace Minuit2
339
340} // namespace ROOT
341
342#endif // ROOT_Minuit2_LAVector
double
Definition: Converters.cxx:939
LAVector & operator+=(const ABObj< vec, ABProd< ABObj< sym, LASymMatrix, T >, ABObj< vec, LAVector, T >>, T > &prod)
Definition: LAVector.h:186
LAVector & operator=(const ABObj< vec, ABProd< ABObj< sym, LASymMatrix, T >, ABObj< vec, LAVector, T >>, T > &prod)
Definition: LAVector.h:299
unsigned int size() const
Definition: LAVector.h:227
const double * Data() const
Definition: LAVector.h:223
double operator[](unsigned int i) const
Definition: LAVector.h:211
LAVector & operator*=(double scal)
Definition: LAVector.h:193
LAVector(unsigned int n)
Definition: LAVector.h:42
double operator()(unsigned int i) const
Definition: LAVector.h:199
LAVector(const ABObj< vec, ABObj< vec, A, T >, T > &something)
Definition: LAVector.h:113
LAVector & operator-=(const LAVector &m)
Definition: LAVector.h:155
LAVector(const ABObj< vec, ABProd< ABObj< sym, LASymMatrix, T >, ABObj< vec, LAVector, T >>, T > &prod)
Definition: LAVector.h:122
double & operator()(unsigned int i)
Definition: LAVector.h:205
LAVector & operator+=(const ABObj< vec, A, T > &m)
Definition: LAVector.h:178
LAVector & operator=(const ABObj< vec, ABSum< ABObj< vec, A, T >, ABObj< vec, B, T >>, T > &sum)
Definition: LAVector.h:266
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:135
LAVector & operator=(const ABObj< vec, ABObj< vec, A, T >, T > &something)
Definition: LAVector.h:250
unsigned int fSize
Definition: LAVector.h:230
LAVector & operator=(const ABObj< vec, LAVector, T > &v)
Definition: LAVector.h:235
double & operator[](unsigned int i)
Definition: LAVector.h:217
LAVector(const ABObj< vec, LAVector, T > &v)
Definition: LAVector.h:77
LAVector & operator=(const LAVector &v)
Definition: LAVector.h:66
LAVector & operator+=(const ABObj< vec, LAVector, T > &m)
Definition: LAVector.h:164
LAVector(const LAVector &v)
Definition: LAVector.h:59
LAVector(const ABObj< vec, ABSum< ABObj< vec, A, T >, ABObj< vec, B, T >>, T > &sum)
Definition: LAVector.h:88
LAVector(const ABObj< vec, ABSum< ABObj< vec, LAVector, T >, ABObj< vec, A, T >>, T > &sum)
Definition: LAVector.h:98
LAVector & operator=(const ABObj< vec, ABSum< ABObj< vec, LAVector, T >, ABObj< vec, A, T >>, T > &sum)
Definition: LAVector.h:282
LAVector & operator+=(const LAVector &m)
Definition: LAVector.h:147
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:318
static StackAllocator & Get()
void * Allocate(size_t nBytes)
const Int_t n
Definition: legend1.C:16
RPY_EXPORTED TCppObject_t Allocate(TCppType_t type)
static double B[]
double T(double x)
Definition: ChebyshevPol.h:34
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...
Definition: RNumpyDS.hxx:30
auto * m
Definition: textangle.C:8
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345