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