Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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(nullptr) {}
49
50public:
51 typedef sym Type;
52
53 LASymMatrix(unsigned int n)
54 : fSize(n * (n + 1) / 2),
55 fNRow(n),
56 fData( (n> 0) ? (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * n * (n + 1) / 2)
57 : nullptr )
58 {
59 // assert(fSize>0);
60 if (fData)
61 std::memset(fData, 0, fSize * sizeof(double));
62 }
63
65 {
66 if (fData)
68 }
69
71 : fSize(v.size()), fNRow(v.Nrow()),
72 fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * v.size()))
73 {
74 std::memcpy(fData, v.Data(), fSize * sizeof(double));
75 }
76
78 {
79 if (fSize < v.size()) {
80 if (fData)
82 fSize = v.size();
83 fNRow = v.Nrow();
84 fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
85 }
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(nullptr)
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(nullptr)
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(nullptr)
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(nullptr)
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) {
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 == nullptr) {
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 == nullptr) {
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 == nullptr) {
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 == nullptr) {
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
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)
LASymMatrix & operator=(const ABObj< sym, LASymMatrix, T > &v)
LASymMatrix & operator-=(const LASymMatrix &m)
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, LASymMatrix, T >, ABObj< sym, A, T > >, T > &sum)
double & operator()(unsigned int row, unsigned int col)
const double * Data() const
unsigned int Ncol() const
double operator()(unsigned int row, unsigned int col) const
LASymMatrix & operator=(const LASymMatrix &v)
Definition LASymMatrix.h:77
LASymMatrix & operator+=(const ABObj< sym, LASymMatrix, T > &m)
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T >, ABObj< sym, A, T > >, T > &sum)
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, A, T >, ABObj< sym, B, T > >, T > &sum)
LASymMatrix & operator*=(double scal)
LASymMatrix & operator+=(const ABObj< sym, A, T > &m)
LASymMatrix & operator=(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &inv)
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)
LASymMatrix(const LASymMatrix &v)
Definition LASymMatrix.h:70
unsigned int Nrow() const
LASymMatrix(const ABObj< sym, ABObj< sym, A, T >, T > &something)
LASymMatrix & operator+=(const ABObj< sym, MatrixInverse< sym, ABObj< sym, LASymMatrix, T >, T >, T > &m)
LASymMatrix & operator=(const ABObj< sym, ABSum< ABObj< sym, LASymMatrix, T >, ABObj< sym, A, T > >, T > &sum)
LASymMatrix(const ABObj< sym, ABSum< ABObj< sym, VectorOuterProduct< ABObj< vec, LAVector, T >, T >, T >, ABObj< sym, A, T > >, T > &sum)
LASymMatrix & operator+=(const LASymMatrix &m)
LASymMatrix & operator+=(const ABObj< sym, VectorOuterProduct< ABObj< vec, LAVector, T >, T >, T > &m)
unsigned int size() const
LASymMatrix & operator=(const ABObj< sym, ABObj< sym, A, T >, T > &something)
static StackAllocator & Get()
void * Allocate(size_t nBytes)
const Int_t n
Definition legend1.C:16
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.)
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.
void inv(rsa_NUMBER *, rsa_NUMBER *, rsa_NUMBER *)
Definition rsaaux.cxx:949
TMarker m
Definition textangle.C:8
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345