ROOT_528-00b_version: include/TMatrixTSym.h Source File

00001 // @(#)root/matrix:$Id: TMatrixTSym.h 34744 2010-08-07 06:16:36Z brun $
00002 // Authors: Fons Rademakers, Eddy Offermann   Nov 2003
00003 
00004 /*************************************************************************
00005  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
00006  * All rights reserved.                                                  *
00007  *                                                                       *
00008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
00009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
00010  *************************************************************************/
00011 
00012 #ifndef ROOT_TMatrixTSym
00013 #define ROOT_TMatrixTSym
00014 
00015 //////////////////////////////////////////////////////////////////////////
00016 //                                                                      //
00017 // TMatrixTSym                                                          //
00018 //                                                                      //
00019 // Implementation of a symmetric matrix in the linear algebra package   //
00020 //                                                                      //
00021 // Note that in this implementation both matrix element m[i][j] and     //
00022 // m[j][i] are updated and stored in memory . However, when making the  //
00023 // object persistent only the upper right triangle is stored .          //
00024 //                                                                      //
00025 //////////////////////////////////////////////////////////////////////////
00026 
00027 #ifndef ROOT_TMatrixTBase
00028 #include "TMatrixTBase.h"
00029 #endif
00030 #ifndef ROOT_TMatrixTUtils
00031 #include "TMatrixTUtils.h"
00032 #endif
00033 
00034 template<class Element>class TMatrixT;
00035 template<class Element>class TMatrixTSymLazy;
00036 template<class Element>class TVectorT;
00037 
00038 template<class Element> class TMatrixTSym : public TMatrixTBase<Element> {
00039 
00040 protected:
00041 
00042    Element  fDataStack[TMatrixTBase<Element>::kSizeMax]; //! data container
00043    Element *fElements;                                   //[fNelems] elements themselves
00044 
00045    Element *New_m   (Int_t size);
00046    void     Delete_m(Int_t size,Element*&);
00047    Int_t    Memcpy_m(Element *newp,const Element *oldp,Int_t copySize,
00048                      Int_t newSize,Int_t oldSize);
00049    void     Allocate(Int_t nrows,Int_t ncols,Int_t row_lwb = 0,Int_t col_lwb = 0,Int_t init = 0,
00050                      Int_t /*nr_nonzeros*/ = -1);
00051 
00052 public:
00053 
00054    enum {kWorkMax = 100}; // size of work array
00055    enum EMatrixCreatorsOp1 { kZero,kUnit,kTransposed,kInverted,kAtA };
00056    enum EMatrixCreatorsOp2 { kPlus,kMinus };
00057 
00058    TMatrixTSym() { fElements = 0; }
00059    explicit TMatrixTSym(Int_t nrows);
00060    TMatrixTSym(Int_t row_lwb,Int_t row_upb);
00061    TMatrixTSym(Int_t nrows,const Element *data,Option_t *option="");
00062    TMatrixTSym(Int_t row_lwb,Int_t row_upb,const Element *data,Option_t *option="");
00063    TMatrixTSym(const TMatrixTSym<Element> &another);
00064    template <class Element2> TMatrixTSym(const TMatrixTSym<Element2> &another)
00065    {
00066       R__ASSERT(another.IsValid());
00067       Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
00068       *this = another;
00069    }
00070 
00071    TMatrixTSym(EMatrixCreatorsOp1 op,const TMatrixTSym<Element> &prototype);
00072    TMatrixTSym(EMatrixCreatorsOp1 op,const TMatrixT   <Element> &prototype);
00073    TMatrixTSym(const TMatrixTSym<Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSym<Element> &b);
00074    TMatrixTSym(const TMatrixTSymLazy<Element> &lazy_constructor);
00075 
00076    virtual ~TMatrixTSym() { Clear(); }
00077 
00078    // Elementary constructors
00079    void TMult(const TMatrixT   <Element> &a);
00080    void TMult(const TMatrixTSym<Element> &a);
00081    void Mult (const TMatrixTSym<Element> &a) { TMult(a); }
00082 
00083    void Plus (const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b);
00084    void Minus(const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b);
00085 
00086    virtual const Element *GetMatrixArray  () const;
00087    virtual       Element *GetMatrixArray  ();
00088    virtual const Int_t   *GetRowIndexArray() const { return 0; }
00089    virtual       Int_t   *GetRowIndexArray()       { return 0; }
00090    virtual const Int_t   *GetColIndexArray() const { return 0; }
00091    virtual       Int_t   *GetColIndexArray()       { return 0; }
00092 
00093    virtual       TMatrixTBase<Element> &SetRowIndexArray(Int_t * /*data*/) { MayNotUse("SetRowIndexArray(Int_t *)"); return *this; }
00094    virtual       TMatrixTBase<Element> &SetColIndexArray(Int_t * /*data*/) { MayNotUse("SetColIndexArray(Int_t *)"); return *this; }
00095 
00096    virtual void   Clear      (Option_t * /*option*/ ="") { if (this->fIsOwner) Delete_m(this->fNelems,fElements);
00097                                                            else fElements = 0; this->fNelems = 0; }
00098    virtual Bool_t IsSymmetric() const { return kTRUE; }
00099 
00100            TMatrixTSym <Element> &Use           (Int_t row_lwb,Int_t row_upb,Element *data);
00101    const   TMatrixTSym <Element> &Use           (Int_t row_lwb,Int_t row_upb,const Element *data) const
00102                                                   { return (const TMatrixTSym<Element>&)
00103                                                            ((const_cast<TMatrixTSym<Element> *>(this))->Use(row_lwb,row_upb,const_cast<Element *>(data))); }
00104            TMatrixTSym <Element> &Use           (Int_t nrows,Element *data);
00105    const   TMatrixTSym <Element> &Use           (Int_t nrows,const Element *data) const;
00106            TMatrixTSym <Element> &Use           (TMatrixTSym<Element> &a);
00107    const   TMatrixTSym <Element> &Use           (const TMatrixTSym<Element> &a) const;
00108 
00109            TMatrixTSym <Element> &GetSub        (Int_t row_lwb,Int_t row_upb,TMatrixTSym<Element> &target,Option_t *option="S") const;
00110    virtual TMatrixTBase<Element> &GetSub        (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
00111                                                 TMatrixTBase<Element> &target,Option_t *option="S") const;
00112            TMatrixTSym <Element>  GetSub        (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Option_t *option="S") const;
00113            TMatrixTSym <Element> &SetSub        (Int_t row_lwb,const TMatrixTBase<Element> &source);
00114    virtual TMatrixTBase<Element> &SetSub        (Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source);
00115 
00116    virtual TMatrixTBase<Element> &SetMatrixArray(const Element *data, Option_t *option="");
00117 
00118    virtual TMatrixTBase<Element> &Shift         (Int_t row_shift,Int_t col_shift);
00119    virtual TMatrixTBase<Element> &ResizeTo      (Int_t nrows,Int_t ncols,Int_t /*nr_nonzeros*/ =-1);
00120    virtual TMatrixTBase<Element> &ResizeTo      (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t /*nr_nonzeros*/ =-1);
00121    inline  TMatrixTBase<Element> &ResizeTo      (const TMatrixTSym<Element> &m) {
00122                                                 return ResizeTo(m.GetRowLwb(),m.GetRowUpb(),m.GetColLwb(),m.GetColUpb()); }
00123 
00124    virtual Double_t      Determinant   () const;
00125    virtual void          Determinant   (Double_t &d1,Double_t &d2) const;
00126 
00127            TMatrixTSym<Element>  &Invert        (Double_t *det=0);
00128            TMatrixTSym<Element>  &InvertFast    (Double_t *det=0);
00129            TMatrixTSym<Element>  &Transpose     (const TMatrixTSym<Element> &source);
00130    inline  TMatrixTSym<Element>  &T             () { return this->Transpose(*this); }
00131            TMatrixTSym<Element>  &Rank1Update   (const TVectorT   <Element> &v,Element alpha=1.0);
00132            TMatrixTSym<Element>  &Similarity    (const TMatrixT   <Element> &n);
00133            TMatrixTSym<Element>  &Similarity    (const TMatrixTSym<Element> &n);
00134            Element                Similarity    (const TVectorT   <Element> &v) const;
00135            TMatrixTSym<Element>  &SimilarityT   (const TMatrixT   <Element> &n);
00136 
00137    // Either access a_ij as a(i,j)
00138    inline       Element                    operator()(Int_t rown,Int_t coln) const;
00139    inline       Element                   &operator()(Int_t rown,Int_t coln);
00140 
00141    // or as a[i][j]
00142    inline const TMatrixTRow_const<Element> operator[](Int_t rown) const { return TMatrixTRow_const<Element>(*this,rown); }
00143    inline       TMatrixTRow      <Element> operator[](Int_t rown)       { return TMatrixTRow      <Element>(*this,rown); }
00144 
00145    TMatrixTSym<Element> &operator= (const TMatrixTSym    <Element> &source);
00146    TMatrixTSym<Element> &operator= (const TMatrixTSymLazy<Element> &source);
00147    template <class Element2> TMatrixTSym<Element> &operator= (const TMatrixTSym<Element2> &source)
00148    {
00149       if (!AreCompatible(*this,source)) {
00150          Error("operator=(const TMatrixTSym2 &)","matrices not compatible");
00151          return *this;
00152       }
00153 
00154       TObject::operator=(source);
00155       const Element2 * const ps = source.GetMatrixArray();
00156             Element  * const pt = this->GetMatrixArray();
00157       for (Int_t i = 0; i < this->fNelems; i++)
00158          pt[i] = ps[i];
00159       this->fTol = source.GetTol();
00160       return *this;
00161    }
00162 
00163    TMatrixTSym<Element> &operator= (Element val);
00164    TMatrixTSym<Element> &operator-=(Element val);
00165    TMatrixTSym<Element> &operator+=(Element val);
00166    TMatrixTSym<Element> &operator*=(Element val);
00167 
00168    TMatrixTSym &operator+=(const TMatrixTSym &source);
00169    TMatrixTSym &operator-=(const TMatrixTSym &source);
00170 
00171    TMatrixTBase<Element> &Apply(const TElementActionT   <Element> &action);
00172    TMatrixTBase<Element> &Apply(const TElementPosActionT<Element> &action);
00173 
00174    virtual TMatrixTBase<Element> &Randomize  (Element alpha,Element beta,Double_t &seed);
00175    virtual TMatrixTSym <Element> &RandomizePD(Element alpha,Element beta,Double_t &seed);
00176 
00177    const TMatrixT<Element> EigenVectors(TVectorT<Element> &eigenValues) const;
00178 
00179    ClassDef(TMatrixTSym,2) // Template of Symmetric Matrix class
00180 };
00181 
00182 template <class Element> inline const Element               *TMatrixTSym<Element>::GetMatrixArray() const { return fElements; }
00183 template <class Element> inline       Element               *TMatrixTSym<Element>::GetMatrixArray()       { return fElements; }
00184 
00185 template <class Element> inline       TMatrixTSym<Element>  &TMatrixTSym<Element>::Use           (Int_t nrows,Element *data) { return Use(0,nrows-1,data); }
00186 template <class Element> inline const TMatrixTSym<Element>  &TMatrixTSym<Element>::Use           (Int_t nrows,const Element *data) const
00187                                                                                                    { return Use(0,nrows-1,data); }
00188 template <class Element> inline       TMatrixTSym<Element>  &TMatrixTSym<Element>::Use           (TMatrixTSym<Element> &a)
00189                                                                                                  { return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetMatrixArray()); }
00190 template <class Element> inline const TMatrixTSym<Element>  &TMatrixTSym<Element>::Use           (const TMatrixTSym<Element> &a) const
00191                                                                                                  { return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetMatrixArray()); }
00192 
00193 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,
00194                                                                                                   Option_t *option) const
00195                                                                                                  {
00196                                                                                                    TMatrixTSym<Element> tmp;
00197                                                                                                    this->GetSub(row_lwb,row_upb,col_lwb,col_upb,tmp,option);
00198                                                                                                    return tmp;
00199                                                                                                  }
00200 
00201 template <class Element> inline Element TMatrixTSym<Element>::operator()(Int_t rown,Int_t coln) const
00202 {
00203    R__ASSERT(this->IsValid());
00204    const Int_t arown = rown-this->fRowLwb;
00205    const Int_t acoln = coln-this->fColLwb;
00206    if (arown >= this->fNrows || arown < 0) {
00207       Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
00208       return 0.0;
00209    }
00210    if (acoln >= this->fNcols || acoln < 0) {
00211       Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
00212       return 0.0;
00213    }
00214    return (fElements[arown*this->fNcols+acoln]);
00215 }
00216 
00217 template <class Element> inline Element &TMatrixTSym<Element>::operator()(Int_t rown,Int_t coln)
00218 {
00219    R__ASSERT(this->IsValid());
00220    const Int_t arown = rown-this->fRowLwb;
00221    const Int_t acoln = coln-this->fColLwb;
00222    if (arown >= this->fNrows || arown < 0) {
00223       Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
00224       return fElements[0];
00225    }
00226    if (acoln >= this->fNcols || acoln < 0) {
00227       Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
00228       return fElements[0];
00229    }
00230    return (fElements[arown*this->fNcols+acoln]);
00231 }
00232 
00233 template <class Element> Bool_t                operator== (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
00234 template <class Element> TMatrixTSym<Element>  operator+  (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
00235 template <class Element> TMatrixTSym<Element>  operator+  (const TMatrixTSym<Element> &source1,      Element                val);
00236 template <class Element> TMatrixTSym<Element>  operator+  (      Element               val    ,const TMatrixTSym<Element>  &source2);
00237 template <class Element> TMatrixTSym<Element>  operator-  (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
00238 template <class Element> TMatrixTSym<Element>  operator-  (const TMatrixTSym<Element> &source1,      Element                val);
00239 template <class Element> TMatrixTSym<Element>  operator-  (      Element               val    ,const TMatrixTSym<Element>  &source2);
00240 template <class Element> TMatrixTSym<Element>  operator*  (const TMatrixTSym<Element> &source,       Element                val    );
00241 template <class Element> TMatrixTSym<Element>  operator*  (      Element               val,    const TMatrixTSym<Element>  &source );
00242 // Preventing warnings with -Weffc++ in GCC since overloading the || and && operators was a design choice.
00243 #if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) >= 40600
00244 #pragma GCC diagnostic push
00245 #pragma GCC diagnostic ignored "-Weffc++"
00246 #endif
00247 template <class Element> TMatrixTSym<Element>  operator&& (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
00248 template <class Element> TMatrixTSym<Element>  operator|| (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
00249 #if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) >= 40600
00250 #pragma GCC diagnostic pop
00251 #endif
00252 template <class Element> TMatrixTSym<Element>  operator>  (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
00253 template <class Element> TMatrixTSym<Element>  operator>= (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
00254 template <class Element> TMatrixTSym<Element>  operator<= (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
00255 template <class Element> TMatrixTSym<Element>  operator<  (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element>  &source2);
00256 
00257 template <class Element> TMatrixTSym<Element> &Add        (TMatrixTSym<Element> &target,      Element               scalar,const TMatrixTSym<Element> &source);
00258 template <class Element> TMatrixTSym<Element> &ElementMult(TMatrixTSym<Element> &target,const TMatrixTSym<Element> &source);
00259 template <class Element> TMatrixTSym<Element> &ElementDiv (TMatrixTSym<Element> &target,const TMatrixTSym<Element> &source);
00260 
00261 #endif