Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 10:12:32

0001 // @(#)root/matrix:$Id$
0002 // Authors: Fons Rademakers, Eddy Offermann   Nov 2003
0003 
0004 /*************************************************************************
0005  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011 
0012 #ifndef ROOT_TVectorT
0013 #define ROOT_TVectorT
0014 
0015 //////////////////////////////////////////////////////////////////////////
0016 //                                                                      //
0017 // TVectorT                                                             //
0018 //                                                                      //
0019 // Template class of Vectors in the linear algebra package              //
0020 //                                                                      //
0021 //////////////////////////////////////////////////////////////////////////
0022 
0023 #include "TMatrixT.h"
0024 #include "TMatrixTSym.h"
0025 #include "TMatrixTSparse.h"
0026 
0027 template<class Element> class TVectorT : public TObject {
0028 
0029 protected:
0030    Int_t    fNrows{0};             // number of rows
0031    Int_t    fRowLwb{0};            // lower bound of the row index
0032    Element *fElements{nullptr};    //[fNrows] elements themselves
0033 
0034    enum {kSizeMax = 5};            // size data container on stack, see New_m(),Delete_m()
0035    enum {kWorkMax = 100};          // size of work array's in several routines
0036 
0037    Element  fDataStack[kSizeMax];  //! data container
0038    Bool_t   fIsOwner{kTRUE};       //!default kTRUE, when Use array kFALSE
0039 
0040    Element* New_m   (Int_t size);
0041    void     Delete_m(Int_t size,Element*&);
0042    Int_t    Memcpy_m(Element *newp,const Element *oldp,Int_t copySize,
0043                      Int_t newSize,Int_t oldSize);
0044 
0045    void     Allocate(Int_t nrows,Int_t row_lwb = 0,Int_t init = 0);
0046 
0047    enum EVectorStatusBits {
0048      kStatus = BIT(14) // set if vector object is valid
0049    };
0050 
0051 public:
0052 
0053    TVectorT() : fNrows(0), fRowLwb(0), fElements(nullptr), fDataStack (), fIsOwner(kTRUE) { }
0054    explicit TVectorT(Int_t n);
0055    TVectorT(Int_t lwb,Int_t upb);
0056    TVectorT(Int_t n,const Element *elements);
0057    TVectorT(Int_t lwb,Int_t upb,const Element *elements);
0058    TVectorT(const TVectorT            <Element> &another);
0059    TVectorT(const TMatrixTRow_const   <Element> &mr);
0060    TVectorT(const TMatrixTColumn_const<Element> &mc);
0061    TVectorT(const TMatrixTDiag_const  <Element> &md);
0062    template <class Element2> TVectorT(const TVectorT<Element2> &another)
0063    {
0064       R__ASSERT(another.IsValid());
0065       Allocate(another.GetUpb()-another.GetLwb()+1,another.GetLwb());
0066       *this = another;
0067    }
0068    TVectorT(Int_t lwb,Int_t upb,Double_t iv1, ...);
0069    ~TVectorT() override { TVectorT::Clear(); }
0070 
0071    inline          Int_t     GetLwb       () const { return fRowLwb; }
0072    inline          Int_t     GetUpb       () const { return fNrows+fRowLwb-1; }
0073    inline          Int_t     GetNrows     () const { return fNrows; }
0074    inline          Int_t     GetNoElements() const { return fNrows; }
0075 
0076    inline          Element  *GetMatrixArray  ()       { return fElements; }
0077    inline const    Element  *GetMatrixArray  () const { return fElements; }
0078 
0079    inline void     Invalidate ()       { SetBit(kStatus); }
0080    inline void     MakeValid  ()       { ResetBit(kStatus); }
0081    inline Bool_t   IsValid    () const { return !TestBit(kStatus); }
0082    inline Bool_t   IsOwner    () const { return fIsOwner; }
0083    inline void     SetElements(const Element *elements) { R__ASSERT(IsValid());
0084                                                           memcpy(fElements,elements,fNrows*sizeof(Element)); }
0085    inline TVectorT<Element> &Shift     (Int_t row_shift)            { fRowLwb += row_shift; return *this; }
0086           TVectorT<Element> &ResizeTo  (Int_t lwb,Int_t upb);
0087    inline TVectorT<Element> &ResizeTo  (Int_t n)                    { return ResizeTo(0,n-1); }
0088    inline TVectorT<Element> &ResizeTo  (const TVectorT<Element> &v) { return ResizeTo(v.GetLwb(),v.GetUpb()); }
0089 
0090           TVectorT<Element> &Use       (Int_t lwb,Int_t upb,Element *data);
0091    const  TVectorT<Element> &Use       (Int_t lwb,Int_t upb,const Element *data) const
0092           { return (const TVectorT<Element>&)(const_cast<TVectorT<Element> *>(this))->Use(lwb,upb,const_cast<Element *>(data)); }
0093           TVectorT<Element> &Use       (Int_t n,Element *data);
0094    const  TVectorT<Element> &Use       (Int_t n,const Element *data) const ;
0095           TVectorT<Element> &Use       (TVectorT<Element> &v);
0096    const  TVectorT<Element> &Use       (const TVectorT<Element> &v) const ;
0097 
0098           TVectorT<Element> &GetSub    (Int_t row_lwb,Int_t row_upb,TVectorT<Element> &target,Option_t *option="S") const;
0099           TVectorT<Element>  GetSub    (Int_t row_lwb,Int_t row_upb,Option_t *option="S") const;
0100           TVectorT<Element> &SetSub    (Int_t row_lwb,const TVectorT<Element> &source);
0101 
0102    TVectorT<Element> &Zero();
0103    TVectorT<Element> &Abs ();
0104    TVectorT<Element> &Sqr ();
0105    TVectorT<Element> &Sqrt();
0106    TVectorT<Element> &Invert();
0107    TVectorT<Element> &SelectNonZeros(const TVectorT<Element> &select);
0108 
0109    Element Norm1   () const;
0110    Element Norm2Sqr() const;
0111    Element NormInf () const;
0112    Int_t   NonZeros() const;
0113    Element Sum     () const;
0114    Element Min     () const;
0115    Element Max     () const;
0116 
0117    inline const Element &operator()(Int_t index) const;
0118    inline       Element &operator()(Int_t index);
0119    inline const Element &operator[](Int_t index) const { return (*this)(index); }
0120    inline       Element &operator[](Int_t index)       { return (*this)(index); }
0121 
0122    TVectorT<Element> &operator= (const TVectorT                <Element> &source);
0123    TVectorT<Element> &operator= (const TMatrixTRow_const       <Element> &mr);
0124    TVectorT<Element> &operator= (const TMatrixTColumn_const    <Element> &mc);
0125    TVectorT<Element> &operator= (const TMatrixTDiag_const      <Element> &md);
0126    TVectorT<Element> &operator= (const TMatrixTSparseRow_const <Element> &md);
0127    TVectorT<Element> &operator= (const TMatrixTSparseDiag_const<Element> &md);
0128    template <class Element2> TVectorT<Element> &operator= (const TVectorT<Element2> &source)
0129    {
0130       if (!AreCompatible(*this,source)) {
0131          Error("operator=(const TVectorT2 &)","vectors not compatible");
0132          return *this;
0133       }
0134 
0135      TObject::operator=(source);
0136      const Element2 * const ps = source.GetMatrixArray();
0137            Element  * const pt = GetMatrixArray();
0138      for (Int_t i = 0; i < this->fNrows; i++)
0139         pt[i] = ps[i];
0140      return *this;
0141    }
0142 
0143    TVectorT<Element> &operator= (Element val);
0144    TVectorT<Element> &operator+=(Element val);
0145    TVectorT<Element> &operator-=(Element val);
0146    TVectorT<Element> &operator*=(Element val);
0147 
0148    TVectorT<Element> &operator+=(const TVectorT      <Element> &source);
0149    TVectorT<Element> &operator-=(const TVectorT      <Element> &source);
0150    TVectorT<Element> &operator*=(const TMatrixT      <Element> &a);
0151    TVectorT<Element> &operator*=(const TMatrixTSym   <Element> &a);
0152    TVectorT<Element> &operator*=(const TMatrixTSparse<Element> &a);
0153 
0154    Bool_t operator==(Element val) const;
0155    Bool_t operator!=(Element val) const;
0156    Bool_t operator< (Element val) const;
0157    Bool_t operator<=(Element val) const;
0158    Bool_t operator> (Element val) const;
0159    Bool_t operator>=(Element val) const;
0160 
0161    Bool_t MatchesNonZeroPattern(const TVectorT<Element> &select);
0162    Bool_t SomePositive         (const TVectorT<Element> &select);
0163    void   AddSomeConstant      (Element val,const TVectorT<Element> &select);
0164 
0165    void   Randomize            (Element alpha,Element beta,Double_t &seed);
0166 
0167    TVectorT<Element> &Apply(const TElementActionT   <Element> &action);
0168    TVectorT<Element> &Apply(const TElementPosActionT<Element> &action);
0169 
0170    void Add(const TVectorT<Element> &v);
0171    void Add(const TVectorT<Element> &v1, const TVectorT<Element> &v2);
0172    void Clear(Option_t * /*option*/ = "") override
0173    {
0174       if (fIsOwner)
0175          Delete_m(fNrows, fElements);
0176       else
0177          fElements = nullptr;
0178       fNrows = 0;
0179    }
0180    void Draw(Option_t *option = "") override;       // *MENU*
0181    void Print(Option_t *option = "") const override;  // *MENU*
0182 
0183    ClassDefOverride(TVectorT,4)  // Template of Vector class
0184 };
0185 
0186 // When building with -fmodules, it instantiates all pending instantiations,
0187 // instead of delaying them until the end of the translation unit.
0188 // We 'got away with' probably because the use and the definition of the
0189 // explicit specialization do not occur in the same TU.
0190 //
0191 // In case we are building with -fmodules, we need to forward declare the
0192 // specialization in order to compile the dictionary G__Matrix.cxx.
0193 template <> TClass *TVectorT<double>::Class();
0194 
0195 template<class Element> inline       TVectorT<Element> &TVectorT<Element>::Use     (Int_t n,Element *data) { return Use(0,n-1,data); }
0196 template<class Element> inline const TVectorT<Element> &TVectorT<Element>::Use     (Int_t n,const Element *data) const { return Use(0,n-1,data); }
0197 template<class Element> inline       TVectorT<Element> &TVectorT<Element>::Use     (TVectorT &v)
0198                                                                                    {
0199                                                                                      R__ASSERT(v.IsValid());
0200                                                                                      return Use(v.GetLwb(),v.GetUpb(),v.GetMatrixArray());
0201                                                                                    }
0202 template<class Element> inline const TVectorT<Element> &TVectorT<Element>::Use     (const TVectorT &v) const
0203                                                                                    {
0204                                                                                      R__ASSERT(v.IsValid());
0205                                                                                      return Use(v.GetLwb(),v.GetUpb(),v.GetMatrixArray());
0206                                                                                    }
0207 template<class Element> inline       TVectorT<Element>  TVectorT<Element>::GetSub  (Int_t row_lwb,Int_t row_upb,Option_t *option) const
0208                                                                                    {
0209                                                                                      TVectorT tmp;
0210                                                                                      this->GetSub(row_lwb,row_upb,tmp,option);
0211                                                                                      return tmp;
0212                                                                                    }
0213 
0214 template<class Element> inline const Element &TVectorT<Element>::operator()(Int_t ind) const
0215 {
0216    // Access a vector element.
0217 
0218    R__ASSERT(IsValid());
0219    const Int_t aind = ind-fRowLwb;
0220    if (aind >= fNrows || aind < 0) {
0221       Error("operator()","Request index(%d) outside vector range of %d - %d",ind,fRowLwb,fRowLwb+fNrows);
0222       return TMatrixTBase<Element>::NaNValue();
0223    }
0224 
0225    return fElements[aind];
0226 }
0227 template<class Element> inline Element &TVectorT<Element>::operator()(Int_t ind)
0228 {
0229    // Access a vector element.
0230 
0231    R__ASSERT(IsValid());
0232    const Int_t aind = ind-fRowLwb;
0233    if (aind >= fNrows || aind < 0) {
0234       Error("operator()","Request index(%d) outside vector range of %d - %d",ind,fRowLwb,fRowLwb+fNrows);
0235       return TMatrixTBase<Element>::NaNValue();
0236    }
0237 
0238    return fElements[aind];
0239 }
0240 inline namespace TMatrixTAutoloadOps {
0241 
0242 template<class Element> Bool_t              operator==  (const TVectorT      <Element>  &source1,const TVectorT <Element>  &source2);
0243 template<class Element> TVectorT<Element>   operator+   (const TVectorT      <Element>  &source1,const TVectorT <Element>  &source2);
0244 template<class Element> TVectorT<Element>   operator-   (const TVectorT      <Element>  &source1,const TVectorT <Element>  &source2);
0245 template<class Element> Element             operator*   (const TVectorT      <Element>  &source1,const TVectorT <Element>  &source2);
0246 template<class Element> TVectorT<Element>   operator*   (const TMatrixT      <Element>  &a,      const TVectorT <Element>  &source);
0247 template<class Element> TVectorT<Element>   operator*   (const TMatrixTSym   <Element>  &a,      const TVectorT <Element>  &source);
0248 template<class Element> TVectorT<Element>   operator*   (const TMatrixTSparse<Element>  &a,      const TVectorT <Element>  &source);
0249 template<class Element> TVectorT<Element>   operator*   (      Element                   val,    const TVectorT <Element>  &source);
0250 template<class Element>
0251 inline
0252 TVectorT<Element> operator*   (const TVectorT <Element>  &source, Element val) { return val * source; }
0253 
0254 template<class Element> Element             Dot         (const TVectorT      <Element>  &source1,const TVectorT <Element>  &source2);
0255 template <class Element1,class Element2>
0256                         TMatrixT<Element1>  OuterProduct(const TVectorT      <Element1> &v1,     const TVectorT <Element2> &v2);
0257 template <class Element1,class Element2,class Element3>
0258                         TMatrixT<Element1> &OuterProduct(      TMatrixT      <Element1> &target, const TVectorT <Element2> &v1,     const TVectorT      <Element3> &v2);
0259 template <class Element1,class Element2,class Element3>
0260                         Element1            Mult        (const TVectorT      <Element1> &v1,     const TMatrixT <Element2> &m,      const TVectorT      <Element3> &v2);
0261 
0262 template<class Element> TVectorT<Element>  &Add         (      TVectorT      <Element>  &target,       Element              scalar, const TVectorT      <Element>  &source);
0263 template<class Element> TVectorT<Element>  &Add         (      TVectorT      <Element>  &target,       Element              scalar, const TMatrixT      <Element>  &a,
0264                                                          const TVectorT<Element> &source);
0265 template<class Element> TVectorT<Element>  &Add         (      TVectorT      <Element>  &target,       Element              scalar, const TMatrixTSym   <Element>  &a,
0266                                                          const TVectorT<Element> &source);
0267 template<class Element> TVectorT<Element>  &Add         (      TVectorT      <Element>  &target,       Element              scalar, const TMatrixTSparse<Element>  &a,
0268                                                          const TVectorT<Element> &source);
0269 template<class Element> TVectorT<Element>  &AddElemMult (      TVectorT      <Element>  &target,       Element              scalar, const TVectorT      <Element>  &source1,
0270                                                          const TVectorT      <Element>  &source2);
0271 template<class Element> TVectorT<Element>  &AddElemMult (      TVectorT      <Element>  &target,       Element              scalar, const TVectorT      <Element>  &source1,
0272                                                          const TVectorT      <Element>  &source2,const TVectorT <Element>  &select);
0273 template<class Element> TVectorT<Element>  &AddElemDiv  (      TVectorT      <Element>  &target,       Element              scalar, const TVectorT      <Element>  &source1,
0274                                                          const TVectorT      <Element>  &source2);
0275 template<class Element> TVectorT<Element>  &AddElemDiv  (      TVectorT      <Element>  &target,       Element              scalar, const TVectorT      <Element>  &source1,
0276                                                          const TVectorT      <Element>  &source2,const TVectorT <Element>  &select);
0277 template<class Element> TVectorT<Element>  &ElementMult (      TVectorT      <Element>  &target, const TVectorT <Element>  &source);
0278 template<class Element> TVectorT<Element>  &ElementMult (      TVectorT      <Element>  &target, const TVectorT <Element>  &source, const TVectorT      <Element>  &select);
0279 template<class Element> TVectorT<Element>  &ElementDiv  (      TVectorT      <Element>  &target, const TVectorT <Element>  &source);
0280 template<class Element> TVectorT<Element>  &ElementDiv  (      TVectorT      <Element>  &target, const TVectorT <Element>  &source, const TVectorT      <Element>  &select);
0281 
0282 template<class Element1,class Element2> Bool_t AreCompatible(const TVectorT<Element1> &v1,const TVectorT<Element2> &v2,Int_t verbose=0);
0283 // Check matrix and vector for compatibility in multiply:  M * v and v * M
0284 template<class Element1,class Element2> Bool_t AreCompatible(const TMatrixT<Element1> &m, const TVectorT<Element2> &v, Int_t verbose=0);
0285 template<class Element1,class Element2> Bool_t AreCompatible(const TVectorT<Element1> &v, const TMatrixT<Element2> &m, Int_t verbose=0);
0286 
0287 template<class Element> void   Compare              (const TVectorT <Element>  &source1,const TVectorT <Element>  &source2);
0288 template<class Element> Bool_t VerifyVectorValue    (const TVectorT <Element>  &m,            Element val,Int_t verbose, Element maxDevAllow);
0289 template<class Element> Bool_t VerifyVectorValue    (const TVectorT <Element>  &m,            Element val,Int_t verbose)
0290                                                      { return VerifyVectorValue(m,val,verbose,Element(0.0)); }
0291 template<class Element> Bool_t VerifyVectorValue    (const TVectorT <Element>  &m,            Element val)
0292                                                      { return VerifyVectorValue(m,val,1,Element(0.0)); }
0293 template<class Element> Bool_t VerifyVectorIdentity (const TVectorT <Element>  &m1,const TVectorT <Element> &m2, Int_t verbose, Element maxDevAllow);
0294 template<class Element> Bool_t VerifyVectorIdentity (const TVectorT <Element>  &m1,const TVectorT <Element> &m2, Int_t verbose)
0295                                                      { return VerifyVectorIdentity(m1,m2,verbose,Element(0.0)); }
0296 template<class Element> Bool_t VerifyVectorIdentity (const TVectorT <Element>  &m1,const TVectorT <Element> &m2)
0297                                                      { return VerifyVectorIdentity(m1,m2,1,Element(0.0)); }
0298 } // inline namespace TMatrixTAutoloadOps
0299 #endif