Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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"
15 #include "Minuit2/VectorOuterProduct.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 
67  ~LASymMatrix() {
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;
72  if(fData) StackAllocatorHolder::Get().Deallocate(fData);
73  }
74 
75  LASymMatrix(const LASymMatrix& v) :
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 
81  LASymMatrix& operator=(const LASymMatrix& v) {
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>
91  LASymMatrix(const ABObj<sym, LASymMatrix, T>& v) :
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>
101  LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, A, T>, ABObj<sym, B, T> >,T>& sum) : fSize(0), fNRow(0), fData(0) {
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>
110  LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix, T>, ABObj<sym, A, T> >,T>& sum) : fSize(0), fNRow(0), fData(0) {
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>
138  LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T>, ABObj<sym, A, T> >, T>& sum) : fSize(0), fNRow(0), fData(0) {
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 
147  LASymMatrix(const ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, double>, double>, double>&);
148 
149  template<class A, class T>
150  LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, T>, T>, T>, ABObj<sym, A, T> >, T>& sum) : fSize(0), fNRow(0), fData(0) {
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 
159  LASymMatrix& operator+=(const LASymMatrix& m) {
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 
166  LASymMatrix& operator-=(const LASymMatrix& m) {
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>
174  LASymMatrix& operator+=(const ABObj<sym, LASymMatrix, T>& m) {
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>
187  LASymMatrix& operator+=(const ABObj<sym, A, T>& m) {
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>
205  LASymMatrix& operator+=(const ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, T>, T>, T>& m) {
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>
252  LASymMatrix& operator=(const ABObj<sym, LASymMatrix, T>& v) {
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