Logo ROOT   master
Reference Guide
LASymMatrix.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_LASymMatrix
11 #define ROOT_Minuit2_LASymMatrix
12 
13 #include "Minuit2/MnConfig.h"
14 #include "Minuit2/ABSum.h"
16 #include "Minuit2/MatrixInverse.h"
17 
18 #include <cassert>
19 #include <memory>
20 
21 
22 // #include <iostream>
23 
24 #include "Minuit2/StackAllocator.h"
25 //extern StackAllocator StackAllocatorHolder::Get();
26 
27 // for memcopy
28 #include <string.h>
29 
30 namespace ROOT {
31 
32  namespace Minuit2 {
33 
34 
35 int Mndaxpy(unsigned int, double, const double*, int, double*, int);
36 int Mndscal(unsigned int, double, double*, int);
37 
38 class LAVector;
39 
40 int Invert ( LASymMatrix & );
41 
42 /**
43  Class describing a symmetric matrix of size n.
44  The size is specified as a run-time argument passed in the
45  constructor.
46  The class uses expression templates for the operations and functions.
47  Only the independent data are kept in the fdata array of size n*(n+1)/2
48  containing the lower triangular data
49  */
50 
51 class LASymMatrix {
52 
53 private:
54 
55  LASymMatrix() : fSize(0), fNRow(0), fData(0) {}
56 
57 public:
58 
59  typedef sym Type;
60 
61  LASymMatrix(unsigned int n) : fSize(n*(n+1)/2), fNRow(n), fData((double*)StackAllocatorHolder::Get().Allocate(sizeof(double)*n*(n+1)/2)) {
62 // assert(fSize>0);
63  memset(fData, 0, fSize*sizeof(double));
64 // std::cout<<"LASymMatrix(unsigned int n), n= "<<n<<std::endl;
65  }
66 
68 // std::cout<<"~LASymMatrix()"<<std::endl;
69 // if(fData) std::cout<<"deleting "<<fSize<<std::endl;
70 // else std::cout<<"no delete"<<std::endl;
71 // if(fData) delete [] fData;
73  }
74 
76  fSize(v.size()), fNRow(v.Nrow()), fData((double*)StackAllocatorHolder::Get().Allocate(sizeof(double)*v.size())) {
77 // std::cout<<"LASymMatrix(const LASymMatrix& v)"<<std::endl;
78  memcpy(fData, v.Data(), fSize*sizeof(double));
79  }
80 
82 // std::cout<<"LASymMatrix& operator=(const LASymMatrix& v)"<<std::endl;
83 // std::cout<<"fSize= "<<fSize<<std::endl;
84 // std::cout<<"v.size()= "<<v.size()<<std::endl;
85  assert(fSize == v.size());
86  memcpy(fData, v.Data(), fSize*sizeof(double));
87  return *this;
88  }
89 
90  template<class T>
92  fSize(v.Obj().size()), fNRow(v.Obj().Nrow()), fData((double*)StackAllocatorHolder::Get().Allocate(sizeof(double)*v.Obj().size())) {
93 // std::cout<<"LASymMatrix(const ABObj<sym, LASymMatrix, T>& v)"<<std::endl;
94  //std::cout<<"allocate "<<fSize<<std::endl;
95  memcpy(fData, v.Obj().Data(), fSize*sizeof(double));
96  Mndscal(fSize, double(v.f()), fData, 1);
97  //std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
98  }
99 
100  template<class A, class B, class T>
102 // std::cout<<"template<class A, class B, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, A, T>, ABObj<sym, B, T> > >& sum)"<<std::endl;
103 // recursive construction
104  (*this) = sum.Obj().A();
105  (*this) += sum.Obj().B();
106  //std::cout<<"leaving template<class A, class B, class T> LASymMatrix(const ABObj..."<<std::endl;
107  }
108 
109  template<class A, class T>
111 // std::cout<<"template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix, T>, ABObj<sym, A, T> >,T>& sum)"<<std::endl;
112 
113  // recursive construction
114  //std::cout<<"(*this)=sum.Obj().B();"<<std::endl;
115  (*this)=sum.Obj().B();
116  //std::cout<<"(*this)+=sum.Obj().A();"<<std::endl;
117  (*this)+=sum.Obj().A();
118  //std::cout<<"leaving template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix,.."<<std::endl;
119  }
120 
121  template<class A, class T>
122  LASymMatrix(const ABObj<sym, ABObj<sym, A, T>, T>& something) : fSize(0), fNRow(0), fData(0) {
123 // std::cout<<"template<class A, class T> LASymMatrix(const ABObj<sym, ABObj<sym, A, T>, T>& something)"<<std::endl;
124  (*this) = something.Obj();
125  (*this) *= something.f();
126  //std::cout<<"leaving template<class A, class T> LASymMatrix(const ABObj<sym, ABObj<sym, A, T>, T>& something)"<<std::endl;
127  }
128 
129  template<class T>
130  LASymMatrix(const ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T>& inv) : fSize(inv.Obj().Obj().Obj().size()), fNRow(inv.Obj().Obj().Obj().Nrow()), fData((double*)StackAllocatorHolder::Get().Allocate(sizeof(double)*inv.Obj().Obj().Obj().size())) {
131  memcpy(fData, inv.Obj().Obj().Obj().Data(), fSize*sizeof(double));
132  Mndscal(fSize, double(inv.Obj().Obj().f()), fData, 1);
133  Invert(*this);
134  Mndscal(fSize, double(inv.f()), fData, 1);
135  }
136 
137  template<class A, class T>
139 // std::cout<<"template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T>, ABObj<sym, A, T> >, T>& sum)"<<std::endl;
140 
141  // recursive construction
142  (*this)=sum.Obj().B();
143  (*this)+=sum.Obj().A();
144  //std::cout<<"leaving template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix,.."<<std::endl;
145  }
146 
148 
149  template<class A, class T>
151 // std::cout<<"template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, T>, T>, T> ABObj<sym, A, T> >,T>& sum)"<<std::endl;
152 
153  // recursive construction
154  (*this)=sum.Obj().B();
155  (*this)+=sum.Obj().A();
156  //std::cout<<"leaving template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix,.."<<std::endl;
157  }
158 
160 // std::cout<<"LASymMatrix& operator+=(const LASymMatrix& m)"<<std::endl;
161  assert(fSize==m.size());
162  Mndaxpy(fSize, 1., m.Data(), 1, fData, 1);
163  return *this;
164  }
165 
167 // std::cout<<"LASymMatrix& operator-=(const LASymMatrix& m)"<<std::endl;
168  assert(fSize==m.size());
169  Mndaxpy(fSize, -1., m.Data(), 1, fData, 1);
170  return *this;
171  }
172 
173  template<class T>
175 // std::cout<<"template<class T> LASymMatrix& operator+=(const ABObj<sym, LASymMatrix, T>& m)"<<std::endl;
176  assert(fSize==m.Obj().size());
177  if(m.Obj().Data()==fData) {
178  Mndscal(fSize, 1.+double(m.f()), fData, 1);
179  } else {
180  Mndaxpy(fSize, double(m.f()), m.Obj().Data(), 1, fData, 1);
181  }
182  //std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
183  return *this;
184  }
185 
186  template<class A, class T>
188 // std::cout<<"template<class A, class T> LASymMatrix& operator+=(const ABObj<sym, A,T>& m)"<<std::endl;
189  (*this) += LASymMatrix(m);
190  return *this;
191  }
192 
193  template<class T>
195 // std::cout<<"template<class T> LASymMatrix& operator+=(const ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T>& m)"<<std::endl;
196  assert(fNRow > 0);
197  LASymMatrix tmp(m.Obj().Obj());
198  Invert(tmp);
199  tmp *= double(m.f());
200  (*this) += tmp;
201  return *this;
202  }
203 
204  template<class T>
206 // std::cout<<"template<class T> LASymMatrix& operator+=(const ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, T>, T>, T>&"<<std::endl;
207  assert(fNRow > 0);
208  Outer_prod(*this, m.Obj().Obj().Obj(), m.f()*m.Obj().Obj().f()*m.Obj().Obj().f());
209  return *this;
210  }
211 
212  LASymMatrix& operator*=(double scal) {
213  Mndscal(fSize, scal, fData, 1);
214  return *this;
215  }
216 
217  double operator()(unsigned int row, unsigned int col) const {
218  assert(row<fNRow && col < fNRow);
219  if(row > col)
220  return fData[col+row*(row+1)/2];
221  else
222  return fData[row+col*(col+1)/2];
223  }
224 
225  double& operator()(unsigned int row, unsigned int col) {
226  assert(row<fNRow && col < fNRow);
227  if(row > col)
228  return fData[col+row*(row+1)/2];
229  else
230  return fData[row+col*(col+1)/2];
231  }
232 
233  const double* Data() const {return fData;}
234 
235  double* Data() {return fData;}
236 
237  unsigned int size() const {return fSize;}
238 
239  unsigned int Nrow() const {return fNRow;}
240 
241  unsigned int Ncol() const {return Nrow();}
242 
243 private:
244 
245  unsigned int fSize;
246  unsigned int fNRow;
247  double* fData;
248 
249 public:
250 
251  template<class T>
253  //std::cout<<"template<class T> LASymMatrix& operator=(ABObj<sym, LASymMatrix, T>& v)"<<std::endl;
254  if(fSize == 0 && fData == 0) {
255  fSize = v.Obj().size();
256  fNRow = v.Obj().Nrow();
257  fData = (double*)StackAllocatorHolder::Get().Allocate(sizeof(double)*fSize);
258  } else {
259  assert(fSize == v.Obj().size());
260  }
261  //std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
262  memcpy(fData, v.Obj().Data(), fSize*sizeof(double));
263  (*this) *= v.f();
264  return *this;
265  }
266 
267  template<class A, class T>
269  //std::cout<<"template<class A, class T> LASymMatrix& operator=(const ABObj<sym, ABObj<sym, A, T>, T>& something)"<<std::endl;
270  if(fSize == 0 && fData == 0) {
271  (*this) = something.Obj();
272  (*this) *= something.f();
273  } else {
274  LASymMatrix tmp(something.Obj());
275  tmp *= something.f();
276  assert(fSize == tmp.size());
277  memcpy(fData, tmp.Data(), fSize*sizeof(double));
278  }
279  //std::cout<<"template<class A, class T> LASymMatrix& operator=(const ABObj<sym, ABObj<sym, A, T>, T>& something)"<<std::endl;
280  return *this;
281  }
282 
283  template<class A, class B, class T>
285  //std::cout<<"template<class A, class B, class T> LASymMatrix& operator=(const ABObj<sym, ABSum<ABObj<sym, A, T>, ABObj<sym, B, T> >,T>& sum)"<<std::endl;
286  // recursive construction
287  if(fSize == 0 && fData == 0) {
288  (*this) = sum.Obj().A();
289  (*this) += sum.Obj().B();
290  (*this) *= sum.f();
291  } else {
292  LASymMatrix tmp(sum.Obj().A());
293  tmp += sum.Obj().B();
294  tmp *= sum.f();
295  assert(fSize == tmp.size());
296  memcpy(fData, tmp.Data(), fSize*sizeof(double));
297  }
298  return *this;
299  }
300 
301  template<class A, class T>
303  //std::cout<<"template<class A, class T> LASymMatrix& operator=(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix, T>, ABObj<sym, A, T> >,T>& sum)"<<std::endl;
304 
305  if(fSize == 0 && fData == 0) {
306  //std::cout<<"fSize == 0 && fData == 0"<<std::endl;
307  (*this) = sum.Obj().B();
308  (*this) += sum.Obj().A();
309  (*this) *= sum.f();
310  } else {
311  //std::cout<<"creating tmp variable"<<std::endl;
312  LASymMatrix tmp(sum.Obj().B());
313  tmp += sum.Obj().A();
314  tmp *= sum.f();
315  assert(fSize == tmp.size());
316  memcpy(fData, tmp.Data(), fSize*sizeof(double));
317  }
318  //std::cout<<"leaving LASymMatrix& operator=(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix..."<<std::endl;
319  return *this;
320  }
321 
322  template<class T>
324  if(fSize == 0 && fData == 0) {
325  fSize = inv.Obj().Obj().Obj().size();
326  fNRow = inv.Obj().Obj().Obj().Nrow();
327  fData = (double*)StackAllocatorHolder::Get().Allocate(sizeof(double)*fSize);
328  memcpy(fData, inv.Obj().Obj().Obj().Data(), fSize*sizeof(double));
329  (*this) *= inv.Obj().Obj().f();
330  Invert(*this);
331  (*this) *= inv.f();
332  } else {
333  LASymMatrix tmp(inv.Obj().Obj());
334  Invert(tmp);
335  tmp *= double(inv.f());
336  assert(fSize == tmp.size());
337  memcpy(fData, tmp.Data(), fSize*sizeof(double));
338  }
339  return *this;
340  }
341 
343 };
344 
345  } // namespace Minuit2
346 
347 } // namespace ROOT
348 
349 #endif // ROOT_Minuit2_LASymMatrix
static long int sum(long int i)
Definition: Factory.cxx:2272
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
double
Definition: Converters.cxx:921
auto * m
Definition: textangle.C:8
Returns the available number of logical cores.
Definition: RNumpyDS.hxx:30
LASymMatrix(unsigned int n)
Definition: LASymMatrix.h:61
int Mndaxpy(unsigned int, double, const double *, int, double *, int)
Definition: mndaxpy.cxx:20
double T(double x)
Definition: ChebyshevPol.h:34
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
LASymMatrix & operator+=(const LASymMatrix &m)
Definition: LASymMatrix.h:159
LASymMatrix(const ABObj< sym, LASymMatrix, T > &v)
Definition: LASymMatrix.h:91
RPY_EXPORTED TCppObject_t Allocate(TCppType_t type)
double operator()(unsigned int row, unsigned int col) const
Definition: LASymMatrix.h:217
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T >, ABObj< sym, A, T > >, T > &sum)
Definition: LASymMatrix.h:138
LASymMatrix & operator+=(const ABObj< sym, VectorOuterProduct< ABObj< vec, LAVector, T >, T >, T > &m)
Definition: LASymMatrix.h:205
LASymMatrix & operator+=(const ABObj< sym, LASymMatrix, T > &m)
Definition: LASymMatrix.h:174
LASymMatrix & operator=(const ABObj< sym, ABObj< sym, A, T >, T > &something)
Definition: LASymMatrix.h:268
void Outer_prod(LASymMatrix &, const LAVector &, double f=1.)
LASymMatrix(const ABObj< sym, ABObj< sym, A, T >, T > &something)
Definition: LASymMatrix.h:122
const double * Data() const
Definition: LASymMatrix.h:233
LASymMatrix(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &inv)
Definition: LASymMatrix.h:130
LASymMatrix & operator=(const ABObj< sym, ABSum< ABObj< sym, A, T >, ABObj< sym, B, T > >, T > &sum)
Definition: LASymMatrix.h:284
LASymMatrix(const LASymMatrix &v)
Definition: LASymMatrix.h:75
LASymMatrix & operator*=(double scal)
Definition: LASymMatrix.h:212
unsigned int Ncol() const
Definition: LASymMatrix.h:241
void * Allocate(size_t nBytes)
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition: rsaaux.cxx:949
unsigned int Nrow() const
Definition: LASymMatrix.h:239
LASymMatrix & operator-=(const LASymMatrix &m)
Definition: LASymMatrix.h:166
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, VectorOuterProduct< ABObj< vec, LAVector, T >, T >, T >, ABObj< sym, A, T > >, T > &sum)
Definition: LASymMatrix.h:150
LASymMatrix & operator=(const ABObj< sym, ABSum< ABObj< sym, LASymMatrix, T >, ABObj< sym, A, T > >, T > &sum)
Definition: LASymMatrix.h:302
LASymMatrix & operator=(const LASymMatrix &v)
Definition: LASymMatrix.h:81
int Mndscal(unsigned int, double, double *, int)
Definition: mndscal.cxx:20
LASymMatrix & operator=(const ABObj< sym, LASymMatrix, T > &v)
Definition: LASymMatrix.h:252
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, LASymMatrix, T >, ABObj< sym, A, T > >, T > &sum)
Definition: LASymMatrix.h:110
double & operator()(unsigned int row, unsigned int col)
Definition: LASymMatrix.h:225
unsigned int size() const
Definition: LASymMatrix.h:237
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, A, T >, ABObj< sym, B, T > >, T > &sum)
Definition: LASymMatrix.h:101
LASymMatrix & operator+=(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &m)
Definition: LASymMatrix.h:194
static StackAllocator & Get()
const Int_t n
Definition: legend1.C:16
LASymMatrix & operator+=(const ABObj< sym, A, T > &m)
Definition: LASymMatrix.h:187
LASymMatrix & operator=(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &inv)
Definition: LASymMatrix.h:323