Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-30 10:22:12 

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_LAVector
0011 #define ROOT_Minuit2_LAVector
0012 
0013 #include "Minuit2/ABSum.h"
0014 #include "Minuit2/ABProd.h"
0015 #include "Minuit2/LASymMatrix.h"
0016 
0017 #include <cassert>
0018 #include <memory>
0019 
0020 #include "Minuit2/StackAllocator.h"
0021 
0022 namespace ROOT {
0023 
0024 namespace Minuit2 {
0025 
0026 // extern StackAllocator StackAllocatorHolder::Get();
0027 
0028 int Mndaxpy(unsigned int, double, const double *, int, double *, int);
0029 int Mndscal(unsigned int, double, double *, int);
0030 int Mndspmv(const char *, unsigned int, double, const double *, const double *, int, double, double *, int);
0031 
0032 class LAVector {
0033 
0034 private:
0035    LAVector() : fSize(0), fData(nullptr) {}
0036 
0037 public:
0038    typedef vec Type;
0039 
0040    LAVector(unsigned int n) :
0041    fSize(n),
0042    fData( (n>0) ? (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * n)
0043                 : nullptr)
0044    {
0045       //     assert(fSize>0);
0046       if (fData)
0047          std::memset(fData, 0, size() * sizeof(double));
0048       //     std::cout<<"LAVector(unsigned int n), n= "<<n<<std::endl;
0049    }
0050 
0051    ~LAVector()
0052    {
0053       //     std::cout<<"~LAVector()"<<std::endl;
0054       //    if(fData) std::cout<<"deleting "<<fSize<<std::endl;
0055       //     else std::cout<<"no delete"<<std::endl;
0056       //     if(fData) delete [] fData;
0057       if (fData)
0058          StackAllocatorHolder::Get().Deallocate(fData);
0059    }
0060 
0061    LAVector(const LAVector &v)
0062       : fSize(v.size()), fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * v.size()))
0063    {
0064       //     std::cout<<"LAVector(const LAVector& v)"<<std::endl;
0065       std::memcpy(fData, v.Data(), fSize * sizeof(double));
0066    }
0067 
0068    LAVector &operator=(const LAVector &v)
0069    {
0070       // implement proper copy constructor in case size is different
0071       if (v.size() > fSize) {
0072          if (fData) StackAllocatorHolder::Get().Deallocate(fData);
0073          fSize = v.size();
0074          fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
0075       }
0076       std::memcpy(fData, v.Data(), fSize * sizeof(double));
0077       return *this;
0078    }
0079 
0080    template <class T>
0081    LAVector(const ABObj<vec, LAVector, T> &v)
0082       : fSize(v.Obj().size()), fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * v.Obj().size()))
0083    {
0084       //     std::cout<<"LAVector(const ABObj<LAVector, T>& v)"<<std::endl;
0085       //     std::cout<<"allocate "<<fSize<<std::endl;
0086       std::memcpy(fData, v.Obj().Data(), fSize * sizeof(T));
0087       (*this) *= T(v.f());
0088       //     std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
0089    }
0090 
0091    template <class A, class B, class T>
0092    LAVector(const ABObj<vec, ABSum<ABObj<vec, A, T>, ABObj<vec, B, T>>, T> &sum) : fSize(0), fData(nullptr)
0093    {
0094       //     std::cout<<"template<class A, class B, class T> LAVector(const ABObj<ABSum<ABObj<A, T>, ABObj<B, T> > >&
0095       //     sum)"<<std::endl;
0096       (*this) = sum.Obj().A();
0097       (*this) += sum.Obj().B();
0098       (*this) *= double(sum.f());
0099    }
0100 
0101    template <class A, class T>
0102    LAVector(const ABObj<vec, ABSum<ABObj<vec, LAVector, T>, ABObj<vec, A, T>>, T> &sum) : fSize(0), fData(nullptr)
0103    {
0104       //     std::cout<<"template<class A, class T> LAVector(const ABObj<ABSum<ABObj<LAVector, T>, ABObj<A, T> >,T>&
0105       //     sum)"<<std::endl;
0106 
0107       // recursive construction
0108       //     std::cout<<"(*this)=sum.Obj().B();"<<std::endl;
0109       (*this) = sum.Obj().B();
0110       //     std::cout<<"(*this)+=sum.Obj().A();"<<std::endl;
0111       (*this) += sum.Obj().A();
0112       (*this) *= double(sum.f());
0113       //     std::cout<<"leaving template<class A, class T> LAVector(const ABObj<ABSum<ABObj<LAVector,.."<<std::endl;
0114    }
0115 
0116    template <class A, class T>
0117    LAVector(const ABObj<vec, ABObj<vec, A, T>, T> &something) : fSize(0), fData(nullptr)
0118    {
0119       //     std::cout<<"template<class A, class T> LAVector(const ABObj<ABObj<A, T>, T>& something)"<<std::endl;
0120       (*this) = something.Obj();
0121       (*this) *= something.f();
0122    }
0123 
0124    //
0125    template <class T>
0126    LAVector(const ABObj<vec, ABProd<ABObj<sym, LASymMatrix, T>, ABObj<vec, LAVector, T>>, T> &prod)
0127       : fSize(prod.Obj().B().Obj().size()),
0128         fData((double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * prod.Obj().B().Obj().size()))
0129    {
0130       //     std::cout<<"template<class T> LAVector(const ABObj<vec, ABProd<ABObj<sym, LASymMatrix, T>, ABObj<vec,
0131       //     LAVector, T> >, T>& prod)"<<std::endl;
0132 
0133       Mndspmv("U", fSize, prod.f() * prod.Obj().A().f() * prod.Obj().B().f(), prod.Obj().A().Obj().Data(),
0134               prod.Obj().B().Obj().Data(), 1, 0., fData, 1);
0135    }
0136 
0137    //
0138    template <class T>
0139    LAVector(const ABObj<
0140             vec,
0141             ABSum<ABObj<vec, ABProd<ABObj<sym, LASymMatrix, T>, ABObj<vec, LAVector, T>>, T>, ABObj<vec, LAVector, T>>,
0142             T> &prod)
0143       : fSize(0), fData(nullptr)
0144    {
0145       (*this) = prod.Obj().B();
0146       (*this) += prod.Obj().A();
0147       (*this) *= double(prod.f());
0148    }
0149 
0150    //
0151    LAVector &operator+=(const LAVector &m)
0152    {
0153       //     std::cout<<"LAVector& operator+=(const LAVector& m)"<<std::endl;
0154       assert(fSize == m.size());
0155       Mndaxpy(fSize, 1., m.Data(), 1, fData, 1);
0156       return *this;
0157    }
0158 
0159    LAVector &operator-=(const LAVector &m)
0160    {
0161       //     std::cout<<"LAVector& operator-=(const LAVector& m)"<<std::endl;
0162       assert(fSize == m.size());
0163       Mndaxpy(fSize, -1., m.Data(), 1, fData, 1);
0164       return *this;
0165    }
0166 
0167    template <class T>
0168    LAVector &operator+=(const ABObj<vec, LAVector, T> &m)
0169    {
0170       //     std::cout<<"template<class T> LAVector& operator+=(const ABObj<LAVector, T>& m)"<<std::endl;
0171       assert(fSize == m.Obj().size());
0172       if (m.Obj().Data() == fData) {
0173          Mndscal(fSize, 1. + double(m.f()), fData, 1);
0174       } else {
0175          Mndaxpy(fSize, double(m.f()), m.Obj().Data(), 1, fData, 1);
0176       }
0177       //     std::cout<<"fData= "<<fData[0]<<" "<<fData[1]<<std::endl;
0178       return *this;
0179    }
0180 
0181    template <class A, class T>
0182    LAVector &operator+=(const ABObj<vec, A, T> &m)
0183    {
0184       //     std::cout<<"template<class A, class T> LAVector& operator+=(const ABObj<A,T>& m)"<<std::endl;
0185       (*this) += LAVector(m);
0186       return *this;
0187    }
0188 
0189    template <class T>
0190    LAVector &operator+=(const ABObj<vec, ABProd<ABObj<sym, LASymMatrix, T>, ABObj<vec, LAVector, T>>, T> &prod)
0191    {
0192       Mndspmv("U", fSize, prod.f() * prod.Obj().A().f() * prod.Obj().B().f(), prod.Obj().A().Obj().Data(),
0193               prod.Obj().B().Data(), 1, 1., fData, 1);
0194       return *this;
0195    }
0196 
0197    LAVector &operator*=(double scal)
0198    {
0199       Mndscal(fSize, scal, fData, 1);
0200       return *this;
0201    }
0202 
0203    double operator()(unsigned int i) const
0204    {
0205       assert(i < fSize);
0206       return fData[i];
0207    }
0208 
0209    double &operator()(unsigned int i)
0210    {
0211       assert(i < fSize);
0212       return fData[i];
0213    }
0214 
0215    double operator[](unsigned int i) const
0216    {
0217       assert(i < fSize);
0218       return fData[i];
0219    }
0220 
0221    double &operator[](unsigned int i)
0222    {
0223       assert(i < fSize);
0224       return fData[i];
0225    }
0226 
0227    const double *Data() const { return fData; }
0228 
0229    double *Data() { return fData; }
0230 
0231    unsigned int size() const { return fSize; }
0232 
0233 private:
0234    unsigned int fSize;
0235    double *fData;
0236 
0237 public:
0238    template <class T>
0239    LAVector &operator=(const ABObj<vec, LAVector, T> &v)
0240    {
0241       //     std::cout<<"template<class T> LAVector& operator=(ABObj<LAVector, T>& v)"<<std::endl;
0242       if (fSize == 0 && !fData) {
0243          fSize = v.Obj().size();
0244          fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
0245       } else {
0246          assert(fSize == v.Obj().size());
0247       }
0248       std::memcpy(fData, v.Obj().Data(), fSize * sizeof(double));
0249       (*this) *= T(v.f());
0250       return *this;
0251    }
0252 
0253    template <class A, class T>
0254    LAVector &operator=(const ABObj<vec, ABObj<vec, A, T>, T> &something)
0255    {
0256       //     std::cout<<"template<class A, class T> LAVector& operator=(const ABObj<ABObj<A, T>, T>&
0257       //     something)"<<std::endl;
0258       if (fSize == 0 && !fData) {
0259          (*this) = something.Obj();
0260       } else {
0261          LAVector tmp(something.Obj());
0262          assert(fSize == tmp.size());
0263          std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
0264       }
0265       (*this) *= something.f();
0266       return *this;
0267    }
0268 
0269    template <class A, class B, class T>
0270    LAVector &operator=(const ABObj<vec, ABSum<ABObj<vec, A, T>, ABObj<vec, B, T>>, T> &sum)
0271    {
0272       if (fSize == 0 && !fData) {
0273          (*this) = sum.Obj().A();
0274          (*this) += sum.Obj().B();
0275       } else {
0276          LAVector tmp(sum.Obj().A());
0277          tmp += sum.Obj().B();
0278          assert(fSize == tmp.size());
0279          std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
0280       }
0281       (*this) *= sum.f();
0282       return *this;
0283    }
0284 
0285    template <class A, class T>
0286    LAVector &operator=(const ABObj<vec, ABSum<ABObj<vec, LAVector, T>, ABObj<vec, A, T>>, T> &sum)
0287    {
0288       if (fSize == 0 && !fData) {
0289          (*this) = sum.Obj().B();
0290          (*this) += sum.Obj().A();
0291       } else {
0292          LAVector tmp(sum.Obj().A());
0293          tmp += sum.Obj().B();
0294          assert(fSize == tmp.size());
0295          std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
0296       }
0297       (*this) *= sum.f();
0298       return *this;
0299    }
0300 
0301    //
0302    template <class T>
0303    LAVector &operator=(const ABObj<vec, ABProd<ABObj<sym, LASymMatrix, T>, ABObj<vec, LAVector, T>>, T> &prod)
0304    {
0305       if (fSize == 0 && !fData) {
0306          fSize = prod.Obj().B().Obj().size();
0307          fData = (double *)StackAllocatorHolder::Get().Allocate(sizeof(double) * fSize);
0308          Mndspmv("U", fSize, double(prod.f() * prod.Obj().A().f() * prod.Obj().B().f()), prod.Obj().A().Obj().Data(),
0309                  prod.Obj().B().Obj().Data(), 1, 0., fData, 1);
0310       } else {
0311          LAVector tmp(prod.Obj().B());
0312          assert(fSize == tmp.size());
0313          Mndspmv("U", fSize, double(prod.f() * prod.Obj().A().f()), prod.Obj().A().Obj().Data(), tmp.Data(), 1, 0.,
0314                  fData, 1);
0315       }
0316       return *this;
0317    }
0318 
0319    //
0320    template <class T>
0321    LAVector &
0322    operator=(const ABObj<
0323              vec,
0324              ABSum<ABObj<vec, ABProd<ABObj<sym, LASymMatrix, T>, ABObj<vec, LAVector, T>>, T>, ABObj<vec, LAVector, T>>,
0325              T> &prod)
0326    {
0327       if (fSize == 0 && !fData) {
0328          (*this) = prod.Obj().B();
0329          (*this) += prod.Obj().A();
0330       } else {
0331          //       std::cout<<"creating tmp variable"<<std::endl;
0332          LAVector tmp(prod.Obj().B());
0333          tmp += prod.Obj().A();
0334          assert(fSize == tmp.size());
0335          std::memcpy(fData, tmp.Data(), fSize * sizeof(double));
0336       }
0337       (*this) *= prod.f();
0338       return *this;
0339    }
0340 };
0341 
0342 } // namespace Minuit2
0343 
0344 } // namespace ROOT
0345 
0346 #endif // ROOT_Minuit2_LAVector
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!