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