Logo ROOT   6.16/01
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"
17
18#include <cassert>
19#include <memory>
20
21
22// #include <iostream>
23
25//extern StackAllocator StackAllocatorHolder::Get();
26
27// for memcopy
28#include <string.h>
29
30namespace ROOT {
31
32 namespace Minuit2 {
33
34
35int Mndaxpy(unsigned int, double, const double*, int, double*, int);
36int Mndscal(unsigned int, double, double*, int);
37
38class LAVector;
39
40int 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
52
53private:
54
55 LASymMatrix() : fSize(0), fNRow(0), fData(0) {}
56
57public:
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
243private:
244
245 unsigned int fSize;
246 unsigned int fNRow;
247 double* fData;
248
249public:
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
SVector< double, 2 > v
Definition: Dict.h:5
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
LASymMatrix(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &inv)
Definition: LASymMatrix.h:130
LASymMatrix & operator=(const ABObj< sym, LASymMatrix, T > &v)
Definition: LASymMatrix.h:252
LASymMatrix & operator-=(const LASymMatrix &m)
Definition: LASymMatrix.h:166
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
const double * Data() const
Definition: LASymMatrix.h:233
unsigned int Ncol() const
Definition: LASymMatrix.h:241
double operator()(unsigned int row, unsigned int col) const
Definition: LASymMatrix.h:217
LASymMatrix & operator=(const LASymMatrix &v)
Definition: LASymMatrix.h:81
LASymMatrix & operator+=(const ABObj< sym, LASymMatrix, T > &m)
Definition: LASymMatrix.h:174
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(const ABObj< sym, ABSum< ABObj< sym, A, T >, ABObj< sym, B, T > >, T > &sum)
Definition: LASymMatrix.h:101
LASymMatrix & operator*=(double scal)
Definition: LASymMatrix.h:212
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
LASymMatrix(unsigned int n)
Definition: LASymMatrix.h:61
LASymMatrix(const ABObj< sym, LASymMatrix, T > &v)
Definition: LASymMatrix.h:91
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
unsigned int Nrow() const
Definition: LASymMatrix.h:239
LASymMatrix(const ABObj< sym, ABObj< sym, A, T >, T > &something)
Definition: LASymMatrix.h:122
LASymMatrix & operator+=(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &m)
Definition: LASymMatrix.h:194
LASymMatrix & operator=(const ABObj< sym, ABSum< ABObj< sym, LASymMatrix, T >, ABObj< sym, A, T > >, T > &sum)
Definition: LASymMatrix.h:302
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 LASymMatrix &m)
Definition: LASymMatrix.h:159
LASymMatrix & operator+=(const ABObj< sym, VectorOuterProduct< ABObj< vec, LAVector, T >, T >, T > &m)
Definition: LASymMatrix.h:205
unsigned int size() const
Definition: LASymMatrix.h:237
LASymMatrix & operator=(const ABObj< sym, ABObj< sym, A, T >, T > &something)
Definition: LASymMatrix.h:268
static StackAllocator & Get()
void * Allocate(size_t nBytes)
const Int_t n
Definition: legend1.C:16
TCppObject_t Allocate(TCppType_t type)
Definition: Cppyy.cxx:271
double T(double x)
Definition: ChebyshevPol.h:34
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
int Mndaxpy(unsigned int, double, const double *, int, double *, int)
Definition: mndaxpy.cxx:20
int Mndscal(unsigned int, double, double *, int)
Definition: mndscal.cxx:20
void Outer_prod(LASymMatrix &, const LAVector &, double f=1.)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition: rsaaux.cxx:949
auto * m
Definition: textangle.C:8
static long int sum(long int i)
Definition: Factory.cxx:2258