Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/root/TMatrixTSym.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

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_TMatrixTSym
0013 #define ROOT_TMatrixTSym
0014 
0015 //////////////////////////////////////////////////////////////////////////
0016 //                                                                      //
0017 // TMatrixTSym                                                          //
0018 //                                                                      //
0019 // Implementation of a symmetric matrix in the linear algebra package   //
0020 //                                                                      //
0021 // Note that in this implementation both matrix element m[i][j] and     //
0022 // m[j][i] are updated and stored in memory. However, when making the   //
0023 // object persistent only the upper right triangle is stored.           //
0024 //                                                                      //
0025 //////////////////////////////////////////////////////////////////////////
0026 
0027 #include "TMatrixTBase.h"
0028 #include "TMatrixTUtils.h"
0029 
0030 #include <cassert>
0031 
0032 template<class Element>class TMatrixT;
0033 template<class Element>class TMatrixTSymLazy;
0034 template<class Element>class TVectorT;
0035 
0036 template<class Element> class TMatrixTSym : public TMatrixTBase<Element> {
0037 
0038 protected:
0039 
0040    Element  fDataStack[TMatrixTBase<Element>::kSizeMax]; //! data container
0041    Element *fElements;                                   //[fNelems] elements themselves
0042 
0043    Element *New_m   (Int_t size);
0044    void     Delete_m(Int_t size,Element*&);
0045    Int_t    Memcpy_m(Element *newp,const Element *oldp,Int_t copySize,
0046                      Int_t newSize,Int_t oldSize);
0047    void     Allocate(Int_t nrows,Int_t ncols,Int_t row_lwb = 0,Int_t col_lwb = 0,Int_t init = 0,
0048                      Int_t /*nr_nonzeros*/ = -1);
0049 
0050 public:
0051 
0052    enum {kWorkMax = 100}; // size of work array
0053    enum EMatrixCreatorsOp1 { kZero,kUnit,kTransposed,kInverted,kAtA };
0054    enum EMatrixCreatorsOp2 { kPlus,kMinus };
0055 
0056    TMatrixTSym() { fElements = nullptr; }
0057    explicit TMatrixTSym(Int_t nrows);
0058    TMatrixTSym(Int_t row_lwb,Int_t row_upb);
0059    TMatrixTSym(Int_t nrows,const Element *data,Option_t *option="");
0060    TMatrixTSym(Int_t row_lwb,Int_t row_upb,const Element *data,Option_t *option="");
0061    TMatrixTSym(const TMatrixTSym<Element> &another);
0062    template <class Element2> TMatrixTSym(const TMatrixTSym<Element2> &another)
0063    {
0064       R__ASSERT(another.IsValid());
0065       Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
0066       *this = another;
0067    }
0068 
0069    TMatrixTSym(EMatrixCreatorsOp1 op,const TMatrixTSym<Element> &prototype);
0070    TMatrixTSym(EMatrixCreatorsOp1 op,const TMatrixT   <Element> &prototype);
0071    TMatrixTSym(const TMatrixTSym<Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSym<Element> &b);
0072    TMatrixTSym(const TMatrixTSymLazy<Element> &lazy_constructor);
0073 
0074    ~TMatrixTSym() override { TMatrixTSym::Clear(); }
0075 
0076    // Elementary constructors
0077    void TMult(const TMatrixT   <Element> &a);
0078    void TMult(const TMatrixTSym<Element> &a);
0079    void Mult (const TMatrixTSym<Element> &a) { TMult(a); }
0080 
0081    void Plus (const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b);
0082    void Minus(const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b);
0083 
0084    inline void SetElement(Int_t rown, Int_t coln, Element val);
0085 
0086    const Element *GetMatrixArray  () const override;
0087          Element *GetMatrixArray  () override;
0088    const Int_t   *GetRowIndexArray() const override { return nullptr; }
0089          Int_t   *GetRowIndexArray() override       { return nullptr; }
0090    const Int_t   *GetColIndexArray() const override { return nullptr; }
0091          Int_t   *GetColIndexArray() override       { return nullptr; }
0092 
0093          TMatrixTBase<Element> &SetRowIndexArray(Int_t * /*data*/) override { MayNotUse("SetRowIndexArray(Int_t *)"); return *this; }
0094          TMatrixTBase<Element> &SetColIndexArray(Int_t * /*data*/) override { MayNotUse("SetColIndexArray(Int_t *)"); return *this; }
0095 
0096    void   Clear      (Option_t * /*option*/ ="") override { if (this->fIsOwner) Delete_m(this->fNelems,fElements);
0097                                                            else fElements = nullptr;
0098                                                            this->fNelems = 0; }
0099    Bool_t IsSymmetric() const override { return kTRUE; }
0100 
0101            TMatrixTSym <Element> &Use           (Int_t row_lwb,Int_t row_upb,Element *data);
0102    const   TMatrixTSym <Element> &Use           (Int_t row_lwb,Int_t row_upb,const Element *data) const
0103                                                   { return (const TMatrixTSym<Element>&)
0104                                                            ((const_cast<TMatrixTSym<Element> *>(this))->Use(row_lwb,row_upb,const_cast<Element *>(data))); }
0105            TMatrixTSym <Element> &Use           (Int_t nrows,Element *data);
0106    const   TMatrixTSym <Element> &Use           (Int_t nrows,const Element *data) const;
0107            TMatrixTSym <Element> &Use           (TMatrixTSym<Element> &a);
0108    const   TMatrixTSym <Element> &Use           (const TMatrixTSym<Element> &a) const;
0109 
0110            TMatrixTSym <Element> &GetSub        (Int_t row_lwb,Int_t row_upb,TMatrixTSym<Element> &target,Option_t *option="S") const;
0111    TMatrixTBase<Element> &GetSub        (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
0112                                                 TMatrixTBase<Element> &target,Option_t *option="S") const override;
0113            TMatrixTSym <Element>  GetSub        (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Option_t *option="S") const;
0114            TMatrixTSym <Element> &SetSub        (Int_t row_lwb,const TMatrixTBase<Element> &source);
0115    TMatrixTBase<Element> &SetSub        (Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source) override;
0116 
0117    TMatrixTBase<Element> &SetMatrixArray(const Element *data, Option_t *option="") override;
0118 
0119    TMatrixTBase<Element> &Shift         (Int_t row_shift,Int_t col_shift) override;
0120    TMatrixTBase<Element> &ResizeTo      (Int_t nrows,Int_t ncols,Int_t /*nr_nonzeros*/ =-1) override;
0121    TMatrixTBase<Element> &ResizeTo      (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t /*nr_nonzeros*/ =-1) override;
0122    inline  TMatrixTBase<Element> &ResizeTo      (const TMatrixTSym<Element> &m) {
0123                                                 return ResizeTo(m.GetRowLwb(),m.GetRowUpb(),m.GetColLwb(),m.GetColUpb()); }
0124 
0125    Double_t      Determinant   () const override;
0126    void          Determinant   (Double_t &d1,Double_t &d2) const override;
0127 
0128            TMatrixTSym<Element>  &Invert        (Double_t *det=nullptr);
0129            TMatrixTSym<Element>  &InvertFast    (Double_t *det=nullptr);
0130            TMatrixTSym<Element>  &Transpose     (const TMatrixTSym<Element> &source);
0131    inline  TMatrixTSym<Element>  &T             () { return this->Transpose(*this); }
0132            TMatrixTSym<Element>  &Rank1Update   (const TVectorT   <Element> &v,Element alpha=1.0);
0133            TMatrixTSym<Element>  &Similarity    (const TMatrixT   <Element> &n);
0134            TMatrixTSym<Element>  &Similarity    (const TMatrixTSym<Element> &n);
0135            Element                Similarity    (const TVectorT   <Element> &v) const;
0136            TMatrixTSym<Element>  &SimilarityT   (const TMatrixT   <Element> &n);
0137 
0138    // Either access a_ij as a(i,j)
0139    inline       Element                    operator()(Int_t rown,Int_t coln) const override;
0140    inline       Element                   &operator()(Int_t rown,Int_t coln) override; ///< Access element a_ij where i=rown and j=coln. \warning Modifying this element by the caller breaks the symmetry of the matrix if a_ji is not modified accordingly.
0141 
0142    // or as a[i][j]
0143    inline const TMatrixTRow_const<Element> operator[](Int_t rown) const { return TMatrixTRow_const<Element>(*this,rown); }
0144    inline       TMatrixTRow      <Element> operator[](Int_t rown)       { return TMatrixTRow      <Element>(*this,rown); }  ///< Access row a_i where i=rown. \note A concatenated call to [coln] allows to access element a_ij where i=rown and j=coln \warning Modifying this row by the caller breaks the symmetry of the matrix if a_j is not modified accordingly.
0145 
0146    TMatrixTSym<Element> &operator= (const TMatrixTSym    <Element> &source);
0147    TMatrixTSym<Element> &operator= (const TMatrixTSymLazy<Element> &source);
0148    template <class Element2> TMatrixTSym<Element> &operator= (const TMatrixTSym<Element2> &source)
0149    {
0150       if (!AreCompatible(*this,source)) {
0151          Error("operator=(const TMatrixTSym2 &)","matrices not compatible");
0152          return *this;
0153       }
0154 
0155       TObject::operator=(source);
0156       const Element2 * const ps = source.GetMatrixArray();
0157             Element  * const pt = this->GetMatrixArray();
0158       for (Int_t i = 0; i < this->fNelems; i++)
0159          pt[i] = ps[i];
0160       this->fTol = source.GetTol();
0161       return *this;
0162    }
0163 
0164    TMatrixTSym<Element> &operator= (Element val);
0165    TMatrixTSym<Element> &operator-=(Element val);
0166    TMatrixTSym<Element> &operator+=(Element val);
0167    TMatrixTSym<Element> &operator*=(Element val);
0168 
0169    TMatrixTSym &operator+=(const TMatrixTSym &source);
0170    TMatrixTSym &operator-=(const TMatrixTSym &source);
0171 
0172    TMatrixTBase<Element> &Apply(const TElementActionT   <Element> &action) override;
0173    TMatrixTBase<Element> &Apply(const TElementPosActionT<Element> &action) override;
0174 
0175    TMatrixTBase<Element> &Randomize  (Element alpha,Element beta,Double_t &seed) override;
0176    virtual TMatrixTSym <Element> &RandomizePD(Element alpha,Element beta,Double_t &seed);
0177 
0178    const TMatrixT<Element> EigenVectors(TVectorT<Element> &eigenValues) const;
0179 
0180    ClassDefOverride(TMatrixTSym,2) // Template of Symmetric Matrix class
0181 };
0182 // When building with -fmodules, it instantiates all pending instantiations,
0183 // instead of delaying them until the end of the translation unit.
0184 // We 'got away with' probably because the use and the definition of the
0185 // explicit specialization do not occur in the same TU.
0186 //
0187 // In case we are building with -fmodules, we need to forward declare the
0188 // specialization in order to compile the dictionary G__Matrix.cxx.
0189 template <> TClass *TMatrixTSym<double>::Class();
0190 
0191 template <class Element> inline const Element               *TMatrixTSym<Element>::GetMatrixArray() const { return fElements; }
0192 template <class Element> inline       Element               *TMatrixTSym<Element>::GetMatrixArray()       { return fElements; }
0193 
0194 template <class Element> inline       TMatrixTSym<Element>  &TMatrixTSym<Element>::Use           (Int_t nrows,Element *data) { return Use(0,nrows-1,data); }
0195 template <class Element> inline const TMatrixTSym<Element>  &TMatrixTSym<Element>::Use           (Int_t nrows,const Element *data) const
0196                                                                                                    { return Use(0,nrows-1,data); }
0197 template <class Element> inline       TMatrixTSym<Element>  &TMatrixTSym<Element>::Use           (TMatrixTSym<Element> &a)
0198                                                                                                  { return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetMatrixArray()); }
0199 template <class Element> inline const TMatrixTSym<Element>  &TMatrixTSym<Element>::Use           (const TMatrixTSym<Element> &a) const
0200                                                                                                  { return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetMatrixArray()); }
0201 
0202 template <class Element> inline       TMatrixTSym<Element>   TMatrixTSym<Element>::GetSub        (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
0203                                                                                                   Option_t *option) const
0204                                                                                                  {
0205                                                                                                    TMatrixTSym<Element> tmp;
0206                                                                                                    this->GetSub(row_lwb,row_upb,col_lwb,col_upb,tmp,option);
0207                                                                                                    return tmp;
0208                                                                                                  }
0209 
0210 template <class Element> inline Element TMatrixTSym<Element>::operator()(Int_t rown,Int_t coln) const
0211 {
0212    R__ASSERT(this->IsValid());
0213    const Int_t arown = rown-this->fRowLwb;
0214    const Int_t acoln = coln-this->fColLwb;
0215    if (arown >= this->fNrows || arown < 0) {
0216       Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
0217       return TMatrixTBase<Element>::NaNValue();
0218    }
0219    if (acoln >= this->fNcols || acoln < 0) {
0220       Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
0221       return TMatrixTBase<Element>::NaNValue();
0222    }
0223    return (fElements[arown*this->fNcols+acoln]);
0224 }
0225 
0226 template <class Element> inline Element &TMatrixTSym<Element>::operator()(Int_t rown,Int_t coln)
0227 {
0228    R__ASSERT(this->IsValid());
0229    const Int_t arown = rown-this->fRowLwb;
0230    const Int_t acoln = coln-this->fColLwb;
0231    if (arown >= this->fNrows || arown < 0) {
0232       Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
0233       return TMatrixTBase<Element>::NaNValue();
0234    }
0235    if (acoln >= this->fNcols || acoln < 0) {
0236       Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
0237       return TMatrixTBase<Element>::NaNValue();
0238    }
0239    return (fElements[arown*this->fNcols+acoln]);
0240 }
0241 
0242 ////////////////////////////////////////////////////////////////////////////////
0243 /// Efficiently sets element (rown,coln) equal to val
0244 /// Index bound checks can be deactivated by defining NDEBUG
0245 
0246 template <class Element>
0247 inline void TMatrixTSym<Element>::SetElement(Int_t rown, Int_t coln, Element val)
0248 {
0249    assert(this->IsValid());
0250    rown = rown - this->fRowLwb;
0251    coln = coln - this->fColLwb;
0252    assert((rown < this->fNrows && rown >= 0) && "SetElement() error: row index outside matrix range");
0253    assert((coln < this->fNcols && coln >= 0) && "SetElement() error: column index outside matrix range");
0254    fElements[rown * this->fNcols + coln] = val;
0255 }
0256 
0257 template <class Element> Bool_t                operator== (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
0258 template <class Element> TMatrixTSym<Element>  operator+  (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
0259 template <class Element> TMatrixTSym<Element>  operator+  (const TMatrixTSym<Element> &source1,      Element                val);
0260 template <class Element> TMatrixTSym<Element>  operator+  (      Element               val    ,const TMatrixTSym<Element>  &source2);
0261 template <class Element> TMatrixTSym<Element>  operator-  (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
0262 template <class Element> TMatrixTSym<Element>  operator-  (const TMatrixTSym<Element> &source1,      Element                val);
0263 template <class Element> TMatrixTSym<Element>  operator-  (      Element               val    ,const TMatrixTSym<Element>  &source2);
0264 template <class Element> TMatrixTSym<Element>  operator*  (const TMatrixTSym<Element> &source,       Element                val    );
0265 template <class Element> TMatrixTSym<Element>  operator*  (      Element               val,    const TMatrixTSym<Element>  &source );
0266 // Preventing warnings with -Weffc++ in GCC since overloading the || and && operators was a design choice.
0267 #if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) >= 40600
0268 #pragma GCC diagnostic push
0269 #pragma GCC diagnostic ignored "-Weffc++"
0270 #endif
0271 template <class Element> TMatrixTSym<Element>  operator&& (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
0272 template <class Element> TMatrixTSym<Element>  operator|| (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
0273 #if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) >= 40600
0274 #pragma GCC diagnostic pop
0275 #endif
0276 template <class Element> TMatrixTSym<Element>  operator>  (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
0277 template <class Element> TMatrixTSym<Element>  operator>= (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
0278 template <class Element> TMatrixTSym<Element>  operator<= (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
0279 template <class Element> TMatrixTSym<Element>  operator<  (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
0280 
0281 template <class Element> TMatrixTSym<Element> &Add        (TMatrixTSym<Element> &target,      Element               scalar,const TMatrixTSym<Element> &source);
0282 template <class Element> TMatrixTSym<Element> &ElementMult(TMatrixTSym<Element> &target,const TMatrixTSym<Element> &source);
0283 template <class Element> TMatrixTSym<Element> &ElementDiv (TMatrixTSym<Element> &target,const TMatrixTSym<Element> &source);
0284 
0285 #endif
[ 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!