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