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 
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>
194  LASymMatrix& operator+=(const ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T>& m) {
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>
268  LASymMatrix& operator=(const ABObj<sym, ABObj<sym, A, T>, T>& something) {
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>
284  LASymMatrix& operator=(const ABObj<sym, ABSum<ABObj<sym, A, T>, ABObj<sym, B, T> >,T>& sum) {
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>
302  LASymMatrix& operator=(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix, T>, ABObj<sym, A, T> >,T>& sum) {
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>
323  LASymMatrix& operator=(const ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T>& inv) {
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 
342  LASymMatrix& operator=(const ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, double>, double>, double>&);
343 };
344 
345  } // namespace Minuit2
346 
347 } // namespace ROOT
348 
349 #endif // ROOT_Minuit2_LASymMatrix
ROOT::Minuit2::LASymMatrix::LASymMatrix
LASymMatrix(const ABObj< sym, ABObj< sym, A, T >, T > &something)
Definition: LASymMatrix.h:130
m
auto * m
Definition: textangle.C:8
n
const Int_t n
Definition: legend1.C:16
ROOT::Minuit2::LASymMatrix::~LASymMatrix
~LASymMatrix()
Definition: LASymMatrix.h:75
ROOT::Minuit2::Invert
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:30
ROOT::Minuit2::LASymMatrix::Data
const double * Data() const
Definition: LASymMatrix.h:241
ROOT::Minuit2::StackAllocatorHolder
Definition: StackAllocator.h:224
sym
#define sym(otri1, otri2)
Definition: triangle.c:932
ROOT::Minuit2::LASymMatrix::LASymMatrix
LASymMatrix()
Definition: LASymMatrix.h:63
ROOT::Minuit2::LASymMatrix::operator-=
LASymMatrix & operator-=(const LASymMatrix &m)
Definition: LASymMatrix.h:174
ROOT::Minuit2::LASymMatrix::operator=
LASymMatrix & operator=(const LASymMatrix &v)
Definition: LASymMatrix.h:89
ROOT::Minuit2::ABObj
Definition: ABObj.h:29
ROOT::Minuit2::Mndaxpy
int Mndaxpy(unsigned int, double, const double *, int, double *, int)
Definition: mndaxpy.cxx:28
ROOT::Minuit2::LASymMatrix
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:59
ROOT::Minuit2::LASymMatrix::operator+=
LASymMatrix & operator+=(const LASymMatrix &m)
Definition: LASymMatrix.h:167
ROOT::Minuit2::StackAllocator::Deallocate
void Deallocate(void *p)
Definition: StackAllocator.h:112
inv
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition: rsaaux.cxx:949
ROOT::Minuit2::LASymMatrix::size
unsigned int size() const
Definition: LASymMatrix.h:245
ROOT::Minuit2::LASymMatrix::LASymMatrix
LASymMatrix(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &inv)
Definition: LASymMatrix.h:138
ROOT::Minuit2::MatrixInverse
Definition: MatrixInverse.h:30
ROOT::Minuit2::Outer_prod
void Outer_prod(LASymMatrix &, const LAVector &, double f=1.)
Definition: LaOuterProduct.cxx:58
v
@ v
Definition: rootcling_impl.cxx:3635
ROOT::Minuit2::Mndscal
int Mndscal(unsigned int, double, double *, int)
Definition: mndscal.cxx:28
ROOT::Minuit2::ABObj< vec, LAVector, double >
Definition: ABObj.h:73
ROOT::Minuit2::VectorOuterProduct
Definition: VectorOuterProduct.h:30
Cppyy::Allocate
RPY_EXPORTED TCppObject_t Allocate(TCppType_t type)
Definition: clingwrapper.cxx:662
ROOT::Minuit2::LASymMatrix::fData
double * fData
Definition: LASymMatrix.h:255
MatrixInverse.h
ABSum.h
ROOT::Minuit2::LASymMatrix::Type
sym Type
Definition: LASymMatrix.h:67
ROOT::Minuit2::LASymMatrix::fNRow
unsigned int fNRow
Definition: LASymMatrix.h:254
double
double
Definition: Converters.cxx:921
StackAllocator.h
ROOT::Minuit2::StackAllocator::Allocate
void * Allocate(size_t nBytes)
Definition: StackAllocator.h:83
ROOT::Minuit2::LASymMatrix::Nrow
unsigned int Nrow() const
Definition: LASymMatrix.h:247
ROOT::Minuit2::StackAllocatorHolder::Get
static StackAllocator & Get()
Definition: StackAllocator.h:232
ROOT::Minuit2::sym
Definition: ABTypes.h:27
ROOT::Minuit2::LASymMatrix::operator()
double operator()(unsigned int row, unsigned int col) const
Definition: LASymMatrix.h:225
VectorOuterProduct.h
sum
static long int sum(long int i)
Definition: Factory.cxx:2272
ROOT::Minuit2::ABSum
Definition: ABSum.h:29
ROOT::Math::Chebyshev::T
double T(double x)
Definition: ChebyshevPol.h:52
ROOT::Minuit2::LASymMatrix::Ncol
unsigned int Ncol() const
Definition: LASymMatrix.h:249
ROOT::Minuit2::LASymMatrix::fSize
unsigned int fSize
Definition: LASymMatrix.h:253
ROOT::Minuit2::LASymMatrix::operator*=
LASymMatrix & operator*=(double scal)
Definition: LASymMatrix.h:220
ROOT
VSD Structures.
Definition: StringConv.hxx:21
MnConfig.h