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"
18
19#include <cassert>
20#include <memory>
21#include <cstring> // for memcopy
22
23// extern StackAllocator StackAllocatorHolder::Get();
24
25namespace ROOT {
26
27namespace Minuit2 {
28
29int Mndaxpy(unsigned int, double, const double *, int, double *, int);
30int Mndscal(unsigned int, double, double *, int);
31
32class LAVector;
33
34int 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
46
47private:
48 LASymMatrix() : fSize(0), fNRow(0), fData(0) {}
49
50public:
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
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
278private:
279 unsigned int fSize;
280 unsigned int fNRow;
281 double *fData;
282
283public:
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
double
Definition: Converters.cxx:939
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:45
LASymMatrix(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &inv)
Definition: LASymMatrix.h:140
LASymMatrix & operator=(const ABObj< sym, LASymMatrix, T > &v)
Definition: LASymMatrix.h:285
LASymMatrix & operator-=(const LASymMatrix &m)
Definition: LASymMatrix.h:191
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, LASymMatrix, T >, ABObj< sym, A, T > >, T > &sum)
Definition: LASymMatrix.h:113
double & operator()(unsigned int row, unsigned int col)
Definition: LASymMatrix.h:259
const double * Data() const
Definition: LASymMatrix.h:268
unsigned int Ncol() const
Definition: LASymMatrix.h:276
double operator()(unsigned int row, unsigned int col) const
Definition: LASymMatrix.h:250
LASymMatrix & operator=(const LASymMatrix &v)
Definition: LASymMatrix.h:80
LASymMatrix & operator+=(const ABObj< sym, LASymMatrix, T > &m)
Definition: LASymMatrix.h:200
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
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, A, T >, ABObj< sym, B, T > >, T > &sum)
Definition: LASymMatrix.h:103
LASymMatrix & operator*=(double scal)
Definition: LASymMatrix.h:244
LASymMatrix & operator+=(const ABObj< sym, A, T > &m)
Definition: LASymMatrix.h:214
LASymMatrix & operator=(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &inv)
Definition: LASymMatrix.h:364
LASymMatrix(unsigned int n)
Definition: LASymMatrix.h:53
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:321
LASymMatrix(const LASymMatrix &v)
Definition: LASymMatrix.h:72
unsigned int Nrow() const
Definition: LASymMatrix.h:274
LASymMatrix(const ABObj< sym, ABObj< sym, A, T >, T > &something)
Definition: LASymMatrix.h:129
LASymMatrix & operator+=(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &m)
Definition: LASymMatrix.h:222
LASymMatrix & operator=(const ABObj< sym, ABSum< ABObj< sym, LASymMatrix, T >, ABObj< sym, A, T > >, T > &sum)
Definition: LASymMatrix.h:341
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, VectorOuterProduct< ABObj< vec, LAVector, T >, T >, T >, ABObj< sym, A, T > >, T > &sum)
Definition: LASymMatrix.h:169
LASymMatrix & operator+=(const LASymMatrix &m)
Definition: LASymMatrix.h:183
LASymMatrix & operator+=(const ABObj< sym, VectorOuterProduct< ABObj< vec, LAVector, T >, T >, T > &m)
Definition: LASymMatrix.h:235
unsigned int size() const
Definition: LASymMatrix.h:272
LASymMatrix & operator=(const ABObj< sym, ABObj< sym, A, T >, T > &something)
Definition: LASymMatrix.h:302
static StackAllocator & Get()
void * Allocate(size_t nBytes)
const Int_t n
Definition: legend1.C:16
RPY_EXPORTED TCppObject_t Allocate(TCppType_t type)
double T(double x)
Definition: ChebyshevPol.h:34
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:21
int Mndaxpy(unsigned int, double, const double *, int, double *, int)
Definition: mndaxpy.cxx:19
int Mndscal(unsigned int, double, double *, int)
Definition: mndscal.cxx:19
void Outer_prod(LASymMatrix &, const LAVector &, double f=1.)
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition: rsaaux.cxx:949
auto * m
Definition: textangle.C:8