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