root/Minuit2/LASymMatrix.h

0001 // @(#)root/minuit2:$Id$
0002 // Authors: M. Winkler, F. James, L. Moneta, A. Zsenei   2003-2005
0003
0004 /**********************************************************************
0005  *                                                                    *
0006  * Copyright (c) 2005 LCG ROOT Math team,  CERN/PH-SFT                *
0007  *                                                                    *
0008  **********************************************************************/
0009
0010 #ifndef ROOT_Minuit2_LASymMatrix
0011 #define ROOT_Minuit2_LASymMatrix
0012
0013 #include "Minuit2/MnConfig.h"
0014 #include "Minuit2/ABSum.h"
0015 #include "Minuit2/VectorOuterProduct.h"
0016 #include "Minuit2/MatrixInverse.h"
0017 #include "Minuit2/StackAllocator.h"
0018
0019 #include <cassert>
0020 #include <memory>
0021 #include <cstring> // for memcopy
0022
0023 // extern StackAllocator StackAllocatorHolder::Get();
0024
0025 namespace ROOT {
0026
0027 namespace Minuit2 {
0028
0029 int Mndaxpy(unsigned int, double, const double *, int, double *, int);
0030 int Mndscal(unsigned int, double, double *, int);
0031
0032 class LAVector;
0033
0034 int Invert(LASymMatrix &);
0035
0036 /**
0037    Class describing a symmetric matrix of size n.
0038    The size is specified as a run-time argument passed in the
0039    constructor.
0040    The class uses expression templates for the operations and functions.
0041    Only the independent data are kept in the fdata array of size n*(n+1)/2
0042    containing the lower triangular data
0043  */
0044
0045 class LASymMatrix {
0046
0047 private:
0048    LASymMatrix() : fSize(0), fNRow(0), fData(nullptr) {}
0049
0050 public:
0051    typedef sym Type;
0052
0053    LASymMatrix(unsigned int n)
0054       : fSize(n * (n + 1) / 2),
0055         fNRow(n),
0056         fData( (n> 0) ? (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * n * (n + 1) / 2)
0057                       : nullptr )
0058    {
0059       //     assert(fSize>0);
0060       if (fData)
0061          std::memset(fData, 0, fSize * sizeof(double));
0062    }
0063
0064    ~LASymMatrix()
0065    {
0066       if (fData)
0067          StackAllocatorHolder::Get().Deallocate(fData);
0068    }
0069
0070    LASymMatrix(const LASymMatrix &v)
0071       : fSize(v.size()), fNRow(v.Nrow()),
0072         fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * v.size()))
0073    {
0074       std::memcpy(fData, v.Data(), fSize * sizeof(double));
0075    }
0076
0077    LASymMatrix &operator=(const LASymMatrix &v)
0078    {
0079       if (fSize < v.size()) {
0080          if (fData)
0081           StackAllocatorHolder::Get().Deallocate(fData);
0082          fSize = v.size();
0083          fNRow = v.Nrow();
0084          fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
0085       }
0086       std::memcpy(fData, v.Data(), fSize * sizeof(double));
0087       return *this;
0088    }
0089
0090    template <class T>
0091    LASymMatrix(const ABObj<sym, LASymMatrix, T> &v)
0092       : fSize(v.Obj().size()), fNRow(v.Obj().Nrow()),
0093         fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * v.Obj().size()))
0094    {
0095       //     std::cout<<"LASymMatrix(const ABObj<sym, LASymMatrix, T>& v)"<<std::endl;
0096       // std::cout<<"allocate "<<fSize<<std::endl;
0097       std::memcpy(fData, v.Obj().Data(), fSize * sizeof(double));
0098       Mndscal(fSize, double(v.f()), fData, 1);
0099       // std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
0100    }
0101
0102    template <class A, class B, class T>
0103    LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, A, T>, ABObj<sym, B, T>>, T> &sum) : fSize(0), fNRow(0), fData(nullptr)
0104    {
0105       //     std::cout<<"template<class A, class B, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, A, T>,
0106       //     ABObj<sym, B, T> > >& sum)"<<std::endl; recursive construction
0107       (*this) = sum.Obj().A();
0108       (*this) += sum.Obj().B();
0109       // std::cout<<"leaving template<class A, class B, class T> LASymMatrix(const ABObj..."<<std::endl;
0110    }
0111
0112    template <class A, class T>
0113    LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix, T>, ABObj<sym, A, T>>, T> &sum)
0114       : fSize(0), fNRow(0), fData(nullptr)
0115    {
0116       //     std::cout<<"template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix, T>,
0117       //     ABObj<sym, A, T> >,T>& sum)"<<std::endl;
0118
0119       // recursive construction
0120       // std::cout<<"(*this)=sum.Obj().B();"<<std::endl;
0121       (*this) = sum.Obj().B();
0122       // std::cout<<"(*this)+=sum.Obj().A();"<<std::endl;
0123       (*this) += sum.Obj().A();
0124       // std::cout<<"leaving template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym,
0125       // LASymMatrix,.."<<std::endl;
0126    }
0127
0128    template <class A, class T>
0129    LASymMatrix(const ABObj<sym, ABObj<sym, A, T>, T> &something) : fSize(0), fNRow(0), fData(nullptr)
0130    {
0131       //     std::cout<<"template<class A, class T> LASymMatrix(const ABObj<sym, ABObj<sym, A, T>, T>&
0132       //     something)"<<std::endl;
0133       (*this) = something.Obj();
0134       (*this) *= something.f();
0135       // std::cout<<"leaving template<class A, class T> LASymMatrix(const ABObj<sym, ABObj<sym, A, T>, T>&
0136       // something)"<<std::endl;
0137    }
0138
0139    template <class T>
0140    LASymMatrix(const ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T> &inv)
0141       : fSize(inv.Obj().Obj().Obj().size()), fNRow(inv.Obj().Obj().Obj().Nrow()),
0142         fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * inv.Obj().Obj().Obj().size()))
0143    {
0144       std::memcpy(fData, inv.Obj().Obj().Obj().Data(), fSize * sizeof(double));
0145       Mndscal(fSize, double(inv.Obj().Obj().f()), fData, 1);
0146       Invert(*this);
0147       Mndscal(fSize, double(inv.f()), fData, 1);
0148    }
0149
0150    template <class A, class T>
0151    LASymMatrix(
0152       const ABObj<sym, ABSum<ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T>, ABObj<sym, A, T>>, T>
0153          &sum)
0154       : fSize(0), fNRow(0), fData(nullptr)
0155    {
0156       //     std::cout<<"template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym, MatrixInverse<sym,
0157       //     ABObj<sym, LASymMatrix, T>, T>, T>, ABObj<sym, A, T> >, T>& sum)"<<std::endl;
0158
0159       // recursive construction
0160       (*this) = sum.Obj().B();
0161       (*this) += sum.Obj().A();
0162       // std::cout<<"leaving template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym,
0163       // LASymMatrix,.."<<std::endl;
0164    }
0165
0166    LASymMatrix(const ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, double>, double>, double> &);
0167
0168    template <class A, class T>
0169    LASymMatrix(
0170       const ABObj<sym, ABSum<ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, T>, T>, T>, ABObj<sym, A, T>>, T> &sum)
0171       : fSize(0), fNRow(0), fData(nullptr)
0172    {
0173       //     std::cout<<"template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym,
0174       //     VectorOuterProduct<ABObj<vec, LAVector, T>, T>, T> ABObj<sym, A, T> >,T>& sum)"<<std::endl;
0175
0176       // recursive construction
0177       (*this) = sum.Obj().B();
0178       (*this) += sum.Obj().A();
0179       // std::cout<<"leaving template<class A, class T> LASymMatrix(const ABObj<sym, ABSum<ABObj<sym,
0180       // LASymMatrix,.."<<std::endl;
0181    }
0182
0183    LASymMatrix &operator+=(const LASymMatrix &m)
0184    {
0185       //     std::cout<<"LASymMatrix& operator+=(const LASymMatrix& m)"<<std::endl;
0186       assert(fSize == m.size());
0187       Mndaxpy(fSize, 1., m.Data(), 1, fData, 1);
0188       return *this;
0189    }
0190
0191    LASymMatrix &operator-=(const LASymMatrix &m)
0192    {
0193       //     std::cout<<"LASymMatrix& operator-=(const LASymMatrix& m)"<<std::endl;
0194       assert(fSize == m.size());
0195       Mndaxpy(fSize, -1., m.Data(), 1, fData, 1);
0196       return *this;
0197    }
0198
0199    template <class T>
0200    LASymMatrix &operator+=(const ABObj<sym, LASymMatrix, T> &m)
0201    {
0202       //     std::cout<<"template<class T> LASymMatrix& operator+=(const ABObj<sym, LASymMatrix, T>& m)"<<std::endl;
0203       assert(fSize == m.Obj().size());
0204       if (m.Obj().Data() == fData) {
0205          Mndscal(fSize, 1. + double(m.f()), fData, 1);
0206       } else {
0207          Mndaxpy(fSize, double(m.f()), m.Obj().Data(), 1, fData, 1);
0208       }
0209       // std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
0210       return *this;
0211    }
0212
0213    template <class A, class T>
0214    LASymMatrix &operator+=(const ABObj<sym, A, T> &m)
0215    {
0216       //     std::cout<<"template<class A, class T> LASymMatrix& operator+=(const ABObj<sym, A,T>& m)"<<std::endl;
0217       (*this) += LASymMatrix(m);
0218       return *this;
0219    }
0220
0221    template <class T>
0222    LASymMatrix &operator+=(const ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T> &m)
0223    {
0224       //     std::cout<<"template<class T> LASymMatrix& operator+=(const ABObj<sym, MatrixInverse<sym, ABObj<sym,
0225       //     LASymMatrix, T>, T>, T>& m)"<<std::endl;
0226       assert(fNRow > 0);
0227       LASymMatrix tmp(m.Obj().Obj());
0228       Invert(tmp);
0229       tmp *= double(m.f());
0230       (*this) += tmp;
0231       return *this;
0232    }
0233
0234    template <class T>
0235    LASymMatrix &operator+=(const ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, T>, T>, T> &m)
0236    {
0237       //     std::cout<<"template<class T> LASymMatrix& operator+=(const ABObj<sym, VectorOuterProduct<ABObj<vec,
0238       //     LAVector, T>, T>, T>&"<<std::endl;
0239       assert(fNRow > 0);
0240       Outer_prod(*this, m.Obj().Obj().Obj(), m.f() * m.Obj().Obj().f() * m.Obj().Obj().f());
0241       return *this;
0242    }
0243
0244    LASymMatrix &operator*=(double scal)
0245    {
0246       Mndscal(fSize, scal, fData, 1);
0247       return *this;
0248    }
0249
0250    double operator()(unsigned int row, unsigned int col) const
0251    {
0252       assert(row < fNRow && col < fNRow);
0253       if (row > col)
0254          return fData[col + row * (row + 1) / 2];
0255       else
0256          return fData[row + col * (col + 1) / 2];
0257    }
0258
0259    double &operator()(unsigned int row, unsigned int col)
0260    {
0261       assert(row < fNRow && col < fNRow);
0262       if (row > col)
0263          return fData[col + row * (row + 1) / 2];
0264       else
0265          return fData[row + col * (col + 1) / 2];
0266    }
0267
0268    const double *Data() const { return fData; }
0269
0270    double *Data() { return fData; }
0271
0272    unsigned int size() const { return fSize; }
0273
0274    unsigned int Nrow() const { return fNRow; }
0275
0276    unsigned int Ncol() const { return Nrow(); }
0277
0278 private:
0279    unsigned int fSize;
0280    unsigned int fNRow;
0281    double *fData;
0282
0283 public:
0284    template <class T>
0285    LASymMatrix &operator=(const ABObj<sym, LASymMatrix, T> &v)
0286    {
0287       // std::cout<<"template<class T> LASymMatrix& operator=(ABObj<sym, LASymMatrix, T>& v)"<<std::endl;
0288       if (fSize == 0 && !fData) {
0289          fSize = v.Obj().size();
0290          fNRow = v.Obj().Nrow();
0291          fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
0292       } else {
0293          assert(fSize == v.Obj().size());
0294       }
0295       // std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
0296       std::memcpy(fData, v.Obj().Data(), fSize * sizeof(double));
0297       (*this) *= v.f();
0298       return *this;
0299    }
0300
0301    template <class A, class T>
0302    LASymMatrix &operator=(const ABObj<sym, ABObj<sym, A, T>, T> &something)
0303    {
0304       // std::cout<<"template<class A, class T> LASymMatrix& operator=(const ABObj<sym, ABObj<sym, A, T>, T>&
0305       // something)"<<std::endl;
0306       if (fSize == 0 && fData == nullptr) {
0307          (*this) = something.Obj();
0308          (*this) *= something.f();
0309       } else {
0310          LASymMatrix tmp(something.Obj());
0311          tmp *= something.f();
0312          assert(fSize == tmp.size());
0313          std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
0314       }
0315       // std::cout<<"template<class A, class T> LASymMatrix& operator=(const ABObj<sym, ABObj<sym, A, T>, T>&
0316       // something)"<<std::endl;
0317       return *this;
0318    }
0319
0320    template <class A, class B, class T>
0321    LASymMatrix &operator=(const ABObj<sym, ABSum<ABObj<sym, A, T>, ABObj<sym, B, T>>, T> &sum)
0322    {
0323       // std::cout<<"template<class A, class B, class T> LASymMatrix& operator=(const ABObj<sym, ABSum<ABObj<sym, A, T>,
0324       // ABObj<sym, B, T> >,T>& sum)"<<std::endl;
0325       // recursive construction
0326       if (fSize == 0 && fData == nullptr) {
0327          (*this) = sum.Obj().A();
0328          (*this) += sum.Obj().B();
0329          (*this) *= sum.f();
0330       } else {
0331          LASymMatrix tmp(sum.Obj().A());
0332          tmp += sum.Obj().B();
0333          tmp *= sum.f();
0334          assert(fSize == tmp.size());
0335          std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
0336       }
0337       return *this;
0338    }
0339
0340    template <class A, class T>
0341    LASymMatrix &operator=(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix, T>, ABObj<sym, A, T>>, T> &sum)
0342    {
0343       // std::cout<<"template<class A, class T> LASymMatrix& operator=(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix,
0344       // T>, ABObj<sym, A, T> >,T>& sum)"<<std::endl;
0345
0346       if (fSize == 0 && fData == nullptr) {
0347          // std::cout<<"fSize == 0 && fData == 0"<<std::endl;
0348          (*this) = sum.Obj().B();
0349          (*this) += sum.Obj().A();
0350          (*this) *= sum.f();
0351       } else {
0352          // std::cout<<"creating tmp variable"<<std::endl;
0353          LASymMatrix tmp(sum.Obj().B());
0354          tmp += sum.Obj().A();
0355          tmp *= sum.f();
0356          assert(fSize == tmp.size());
0357          std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
0358       }
0359       // std::cout<<"leaving LASymMatrix& operator=(const ABObj<sym, ABSum<ABObj<sym, LASymMatrix..."<<std::endl;
0360       return *this;
0361    }
0362
0363    template <class T>
0364    LASymMatrix &operator=(const ABObj<sym, MatrixInverse<sym, ABObj<sym, LASymMatrix, T>, T>, T> &inv)
0365    {
0366       if (fSize == 0 && fData == nullptr) {
0367          fSize = inv.Obj().Obj().Obj().size();
0368          fNRow = inv.Obj().Obj().Obj().Nrow();
0369          fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
0370          std::memcpy(fData, inv.Obj().Obj().Obj().Data(), fSize * sizeof(double));
0371          (*this) *= inv.Obj().Obj().f();
0372          Invert(*this);
0373          (*this) *= inv.f();
0374       } else {
0375          LASymMatrix tmp(inv.Obj().Obj());
0376          Invert(tmp);
0377          tmp *= double(inv.f());
0378          assert(fSize == tmp.size());
0379          std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
0380       }
0381       return *this;
0382    }
0383
0384    LASymMatrix &operator=(const ABObj<sym, VectorOuterProduct<ABObj<vec, LAVector, double>, double>, double> &);
0385 };
0386
0387 } // namespace Minuit2
0388
0389 } // namespace ROOT
0390
0391 #endif // ROOT_Minuit2_LASymMatrix