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 #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:
LASymMatrix()48    LASymMatrix() : fSize(0), fNRow(0), fData(0) {}
49 
50 public:
51    typedef sym Type;
52 
LASymMatrix(unsigned int n)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 
~LASymMatrix()62    ~LASymMatrix()
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)
69          StackAllocatorHolder::Get().Deallocate(fData);
70    }
71 
LASymMatrix(const LASymMatrix & v)72    LASymMatrix(const LASymMatrix &v)
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 
80    LASymMatrix &operator=(const LASymMatrix &v)
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>
LASymMatrix(const ABObj<sym,LASymMatrix,T> & v)91    LASymMatrix(const ABObj<sym, LASymMatrix, T> &v)
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>
LASymMatrix(const ABObj<sym,ABSum<ABObj<sym,A,T>,ABObj<sym,B,T>>,T> & sum)103    LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, A, T>, ABObj<sym, B, T>>, T> &sum) : fSize(0), fNRow(0), fData(0)
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>
LASymMatrix(const ABObj<sym,ABSum<ABObj<sym,LASymMatrix,T>,ABObj<sym,A,T>>,T> & sum)113    LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix, T>, ABObj<sym, A, T>>, T> &sum)
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>
LASymMatrix(const ABObj<sym,ABObj<sym,A,T>,T> & something)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>
LASymMatrix(const ABObj<sym,MatrixInverse<sym,ABObj<sym,LASymMatrix,T>,T>,T> & inv)140    LASymMatrix(const ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T> &inv)
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>
LASymMatrix(const ABObj<sym,ABSum<ABObj<sym,MatrixInverse<sym,ABObj<sym,LASymMatrix,T>,T>,T>,ABObj<sym,A,T>>,T> & sum)151    LASymMatrix(
152       const ABObj<sym, ABSum<ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T>, ABObj<sym, A, T>>, 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 
166    LASymMatrix(const ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, double>, double>, double> &);
167 
168    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)169    LASymMatrix(
170       const ABObj<sym, ABSum<ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, T>, T>, T>, ABObj<sym, A, T>>, T> &sum)
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 
183    LASymMatrix &operator+=(const LASymMatrix &m)
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 
191    LASymMatrix &operator-=(const LASymMatrix &m)
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>
200    LASymMatrix &operator+=(const ABObj<sym, LASymMatrix, T> &m)
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>
214    LASymMatrix &operator+=(const ABObj<sym, A, T> &m)
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>
222    LASymMatrix &operator+=(const ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T> &m)
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>
235    LASymMatrix &operator+=(const ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, T>, T>, T> &m)
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 
operator()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 
operator()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 
Data()268    const double *Data() const { return fData; }
269 
Data()270    double *Data() { return fData; }
271 
size()272    unsigned int size() const { return fSize; }
273 
Nrow()274    unsigned int Nrow() const { return fNRow; }
275 
Ncol()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>
285    LASymMatrix &operator=(const ABObj<sym, LASymMatrix, T> &v)
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>
302    LASymMatrix &operator=(const ABObj<sym, ABObj<sym, A, T>, T> &something)
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>
321    LASymMatrix &operator=(const ABObj<sym, ABSum<ABObj<sym, A, T>, ABObj<sym, B, T>>, T> &sum)
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>
341    LASymMatrix &operator=(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix, T>, ABObj<sym, A, T>>, T> &sum)
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>
364    LASymMatrix &operator=(const ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T> &inv)
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 
384    LASymMatrix &operator=(const ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, double>, double>, double> &);
385 };
386 
387 } // namespace Minuit2
388 
389 } // namespace ROOT
390 
391 #endif // ROOT_Minuit2_LASymMatrix
392