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

00001 // @(#)root/matrix:$Id: TMatrixT.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_TMatrixT
00013 #define ROOT_TMatrixT
00014 
00015 //////////////////////////////////////////////////////////////////////////
00016 //                                                                      //
00017 // TMatrixT                                                             //
00018 //                                                                      //
00019 // Template class of a general matrix in the linear algebra package     //
00020 //                                                                      //
00021 //////////////////////////////////////////////////////////////////////////
00022 
00023 #ifndef ROOT_TMatrixTBase
00024 #include "TMatrixTBase.h"
00025 #endif
00026 #ifndef ROOT_TMatrixTUtils
00027 #include "TMatrixTUtils.h"
00028 #endif
00029 
00030 #ifdef CBLAS
00031 #include <vecLib/vBLAS.h>
00032 //#include <cblas.h>
00033 #endif
00034 
00035 template<class Element> class TMatrixTSym;
00036 template<class Element> class TMatrixTSparse;
00037 template<class Element> class TMatrixTLazy;
00038 
00039 template<class Element> class TMatrixT : public TMatrixTBase<Element> {
00040 
00041 protected:
00042 
00043    Element  fDataStack[TMatrixTBase<Element>::kSizeMax]; //! data container
00044    Element *fElements;                                   //[fNelems] elements themselves
00045 
00046    Element *New_m   (Int_t size);
00047    void     Delete_m(Int_t size,Element*&);
00048    Int_t    Memcpy_m(Element *newp,const Element *oldp,Int_t copySize,
00049                       Int_t newSize,Int_t oldSize);
00050    void     Allocate(Int_t nrows,Int_t ncols,Int_t row_lwb = 0,Int_t col_lwb = 0,Int_t init = 0,
00051                      Int_t /*nr_nonzeros*/ = -1);
00052 
00053 public:
00054 
00055    enum {kWorkMax = 100};
00056    enum EMatrixCreatorsOp1 { kZero,kUnit,kTransposed,kInverted,kAtA };
00057    enum EMatrixCreatorsOp2 { kMult,kTransposeMult,kInvMult,kMultTranspose,kPlus,kMinus };
00058 
00059    TMatrixT(): fElements(0) { }
00060    TMatrixT(Int_t nrows,Int_t ncols);
00061    TMatrixT(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
00062    TMatrixT(Int_t nrows,Int_t ncols,const Element *data,Option_t *option="");
00063    TMatrixT(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,const Element *data,Option_t *option="");
00064    TMatrixT(const TMatrixT      <Element> &another);
00065    TMatrixT(const TMatrixTSym   <Element> &another);
00066    TMatrixT(const TMatrixTSparse<Element> &another);
00067    template <class Element2> TMatrixT(const TMatrixT<Element2> &another): fElements(0)
00068    {
00069       R__ASSERT(another.IsValid());
00070       Allocate(another.GetNrows(),another.GetNcols(),another.GetRowLwb(),another.GetColLwb());
00071       *this = another;
00072    }
00073 
00074    TMatrixT(EMatrixCreatorsOp1 op,const TMatrixT<Element> &prototype);
00075    TMatrixT(const TMatrixT    <Element> &a,EMatrixCreatorsOp2 op,const TMatrixT   <Element> &b);
00076    TMatrixT(const TMatrixT    <Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSym<Element> &b);
00077    TMatrixT(const TMatrixTSym <Element> &a,EMatrixCreatorsOp2 op,const TMatrixT   <Element> &b);
00078    TMatrixT(const TMatrixTSym <Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSym<Element> &b);
00079    TMatrixT(const TMatrixTLazy<Element> &lazy_constructor);
00080 
00081    virtual ~TMatrixT() { Clear(); }
00082 
00083    // Elementary constructors
00084 
00085    void Plus (const TMatrixT   <Element> &a,const TMatrixT   <Element> &b);
00086    void Plus (const TMatrixT   <Element> &a,const TMatrixTSym<Element> &b);
00087    void Plus (const TMatrixTSym<Element> &a,const TMatrixT   <Element> &b) { Plus(b,a); }
00088 
00089    void Minus(const TMatrixT   <Element> &a,const TMatrixT   <Element> &b);
00090    void Minus(const TMatrixT   <Element> &a,const TMatrixTSym<Element> &b);
00091    void Minus(const TMatrixTSym<Element> &a,const TMatrixT   <Element> &b) { Minus(b,a); }
00092 
00093    void Mult (const TMatrixT   <Element> &a,const TMatrixT   <Element> &b);
00094    void Mult (const TMatrixT   <Element> &a,const TMatrixTSym<Element> &b);
00095    void Mult (const TMatrixTSym<Element> &a,const TMatrixT   <Element> &b);
00096    void Mult (const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b);
00097 
00098    void TMult(const TMatrixT   <Element> &a,const TMatrixT   <Element> &b);
00099    void TMult(const TMatrixT   <Element> &a,const TMatrixTSym<Element> &b);
00100    void TMult(const TMatrixTSym<Element> &a,const TMatrixT   <Element> &b) { Mult(a,b); }
00101    void TMult(const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b) { Mult(a,b); }
00102 
00103    void MultT(const TMatrixT   <Element> &a,const TMatrixT   <Element> &b);
00104    void MultT(const TMatrixT   <Element> &a,const TMatrixTSym<Element> &b) { Mult(a,b); }
00105    void MultT(const TMatrixTSym<Element> &a,const TMatrixT   <Element> &b);
00106    void MultT(const TMatrixTSym<Element> &a,const TMatrixTSym<Element> &b) { Mult(a,b); }
00107 
00108    virtual const Element *GetMatrixArray  () const;
00109    virtual       Element *GetMatrixArray  ();
00110    virtual const Int_t   *GetRowIndexArray() const { return 0; }
00111    virtual       Int_t   *GetRowIndexArray()       { return 0; }
00112    virtual const Int_t   *GetColIndexArray() const { return 0; }
00113    virtual       Int_t   *GetColIndexArray()       { return 0; }
00114 
00115    virtual       TMatrixTBase<Element> &SetRowIndexArray(Int_t * /*data*/) { MayNotUse("SetRowIndexArray(Int_t *)"); return *this; }
00116    virtual       TMatrixTBase<Element> &SetColIndexArray(Int_t * /*data*/) { MayNotUse("SetColIndexArray(Int_t *)"); return *this; }
00117 
00118    virtual void Clear(Option_t * /*option*/ ="") { if (this->fIsOwner) Delete_m(this->fNelems,fElements);
00119                                                    else fElements = 0;  this->fNelems = 0; }
00120 
00121            TMatrixT    <Element> &Use     (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Element *data);
00122    const   TMatrixT    <Element> &Use     (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,const Element *data) const
00123                                             { return (const TMatrixT<Element>&)
00124                                                      ((const_cast<TMatrixT<Element> *>(this))->Use(row_lwb,row_upb,col_lwb,col_upb, const_cast<Element *>(data))); }
00125            TMatrixT    <Element> &Use     (Int_t nrows,Int_t ncols,Element *data);
00126    const   TMatrixT    <Element> &Use     (Int_t nrows,Int_t ncols,const Element *data) const;
00127            TMatrixT    <Element> &Use     (TMatrixT<Element> &a);
00128    const   TMatrixT    <Element> &Use     (const TMatrixT<Element> &a) const;
00129 
00130    virtual TMatrixTBase<Element> &GetSub  (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
00131                                            TMatrixTBase<Element> &target,Option_t *option="S") const;
00132            TMatrixT    <Element>  GetSub  (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Option_t *option="S") const;
00133    virtual TMatrixTBase<Element> &SetSub  (Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source);
00134 
00135    virtual TMatrixTBase<Element> &ResizeTo(Int_t nrows,Int_t ncols,Int_t /*nr_nonzeros*/ =-1);
00136    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);
00137    inline  TMatrixTBase<Element> &ResizeTo(const TMatrixT<Element> &m) {
00138                                             return ResizeTo(m.GetRowLwb(),m.GetRowUpb(),m.GetColLwb(),m.GetColUpb());
00139                                  }
00140 
00141    virtual Double_t Determinant  () const;
00142    virtual void     Determinant  (Double_t &d1,Double_t &d2) const;
00143 
00144            TMatrixT<Element> &Invert      (Double_t *det=0);
00145            TMatrixT<Element> &InvertFast  (Double_t *det=0);
00146            TMatrixT<Element> &Transpose   (const TMatrixT<Element> &source);
00147    inline  TMatrixT<Element> &T           () { return this->Transpose(*this); }
00148            TMatrixT<Element> &Rank1Update (const TVectorT<Element> &v,Element alpha=1.0);
00149            TMatrixT<Element> &Rank1Update (const TVectorT<Element> &v1,const TVectorT<Element> &v2,Element alpha=1.0);
00150            Element            Similarity  (const TVectorT<Element> &v) const;
00151 
00152    TMatrixT<Element> &NormByColumn(const TVectorT<Element> &v,Option_t *option="D");
00153    TMatrixT<Element> &NormByRow   (const TVectorT<Element> &v,Option_t *option="D");
00154 
00155    // Either access a_ij as a(i,j)
00156    inline       Element                     operator()(Int_t rown,Int_t coln) const;
00157    inline       Element                    &operator()(Int_t rown,Int_t coln);
00158 
00159    // or as a[i][j]
00160    inline const TMatrixTRow_const<Element>  operator[](Int_t rown) const { return TMatrixTRow_const<Element>(*this,rown); }
00161    inline       TMatrixTRow      <Element>  operator[](Int_t rown)       { return TMatrixTRow      <Element>(*this,rown); }
00162 
00163    TMatrixT<Element> &operator= (const TMatrixT      <Element> &source);
00164    TMatrixT<Element> &operator= (const TMatrixTSym   <Element> &source);
00165    TMatrixT<Element> &operator= (const TMatrixTSparse<Element> &source);
00166    TMatrixT<Element> &operator= (const TMatrixTLazy  <Element> &source);
00167    template <class Element2> TMatrixT<Element> &operator= (const TMatrixT<Element2> &source)
00168    {
00169       if (!AreCompatible(*this,source)) {
00170          Error("operator=(const TMatrixT2 &)","matrices not compatible");
00171          return *this;
00172       }
00173 
00174      TObject::operator=(source);
00175      const Element2 * const ps = source.GetMatrixArray();
00176            Element  * const pt = this->GetMatrixArray();
00177      for (Int_t i = 0; i < this->fNelems; i++)
00178         pt[i] = ps[i];
00179      this->fTol = source.GetTol();
00180      return *this;
00181    }
00182 
00183    TMatrixT<Element> &operator= (Element val);
00184    TMatrixT<Element> &operator-=(Element val);
00185    TMatrixT<Element> &operator+=(Element val);
00186    TMatrixT<Element> &operator*=(Element val);
00187 
00188    TMatrixT<Element> &operator+=(const TMatrixT   <Element> &source);
00189    TMatrixT<Element> &operator+=(const TMatrixTSym<Element> &source);
00190    TMatrixT<Element> &operator-=(const TMatrixT   <Element> &source);
00191    TMatrixT<Element> &operator-=(const TMatrixTSym<Element> &source);
00192 
00193    TMatrixT<Element> &operator*=(const TMatrixT            <Element> &source);
00194    TMatrixT<Element> &operator*=(const TMatrixTSym         <Element> &source);
00195    TMatrixT<Element> &operator*=(const TMatrixTDiag_const  <Element> &diag);
00196    TMatrixT<Element> &operator/=(const TMatrixTDiag_const  <Element> &diag);
00197    TMatrixT<Element> &operator*=(const TMatrixTRow_const   <Element> &row);
00198    TMatrixT<Element> &operator/=(const TMatrixTRow_const   <Element> &row);
00199    TMatrixT<Element> &operator*=(const TMatrixTColumn_const<Element> &col);
00200    TMatrixT<Element> &operator/=(const TMatrixTColumn_const<Element> &col);
00201 
00202    const TMatrixT<Element> EigenVectors(TVectorT<Element> &eigenValues) const;
00203 
00204    ClassDef(TMatrixT,4) // Template of General Matrix class
00205 };
00206 
00207 template <class Element> inline const Element           *TMatrixT<Element>::GetMatrixArray() const { return fElements; }
00208 template <class Element> inline       Element           *TMatrixT<Element>::GetMatrixArray()       { return fElements; }
00209 
00210 template <class Element> inline       TMatrixT<Element> &TMatrixT<Element>::Use           (Int_t nrows,Int_t ncols,Element *data)
00211                                                                                           { return Use(0,nrows-1,0,ncols-1,data); }
00212 template <class Element> inline const TMatrixT<Element> &TMatrixT<Element>::Use           (Int_t nrows,Int_t ncols,const Element *data) const
00213                                                                                           { return Use(0,nrows-1,0,ncols-1,data); }
00214 template <class Element> inline       TMatrixT<Element> &TMatrixT<Element>::Use           (TMatrixT &a)
00215                                                                                           {
00216                                                                                             R__ASSERT(a.IsValid());
00217                                                                                             return Use(a.GetRowLwb(),a.GetRowUpb(),
00218                                                                                                        a.GetColLwb(),a.GetColUpb(),a.GetMatrixArray());
00219                                                                                           }
00220 template <class Element> inline const TMatrixT<Element> &TMatrixT<Element>::Use           (const TMatrixT &a) const
00221                                                                                           {
00222                                                                                             R__ASSERT(a.IsValid());
00223                                                                                             return Use(a.GetRowLwb(),a.GetRowUpb(),
00224                                                                                                        a.GetColLwb(),a.GetColUpb(),a.GetMatrixArray());
00225                                                                                           }
00226 
00227 template <class Element> inline       TMatrixT<Element>  TMatrixT<Element>::GetSub        (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
00228                                                                                            Option_t *option) const
00229                                                                                           {
00230                                                                                             TMatrixT tmp;
00231                                                                                             this->GetSub(row_lwb,row_upb,col_lwb,col_upb,tmp,option);
00232                                                                                             return tmp;
00233                                                                                           }
00234 
00235 template <class Element> inline Element TMatrixT<Element>::operator()(Int_t rown,Int_t coln) const
00236 {
00237    R__ASSERT(this->IsValid());
00238    const Int_t arown = rown-this->fRowLwb;
00239    const Int_t acoln = coln-this->fColLwb;
00240    if (arown >= this->fNrows || arown < 0) {
00241       Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
00242       return fElements[0];
00243    }
00244    if (acoln >= this->fNcols || acoln < 0) {
00245       Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
00246       return fElements[0];
00247    }
00248    return (fElements[arown*this->fNcols+acoln]);
00249 }
00250 
00251 template <class Element> inline Element &TMatrixT<Element>::operator()(Int_t rown,Int_t coln)
00252 {
00253    R__ASSERT(this->IsValid());
00254    const Int_t arown = rown-this->fRowLwb;
00255    const Int_t acoln = coln-this->fColLwb;
00256    if (arown >= this->fNrows || arown < 0) {
00257       Error("operator()","Request row(%d) outside matrix range of %d - %d",rown,this->fRowLwb,this->fRowLwb+this->fNrows);
00258       return fElements[0];
00259    }
00260    if (acoln >= this->fNcols || acoln < 0) {
00261       Error("operator()","Request column(%d) outside matrix range of %d - %d",coln,this->fColLwb,this->fColLwb+this->fNcols);
00262       return fElements[0];
00263    }
00264    return (fElements[arown*this->fNcols+acoln]);
00265 }
00266 
00267 template <class Element> TMatrixT<Element>  operator+  (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
00268 template <class Element> TMatrixT<Element>  operator+  (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
00269 template <class Element> TMatrixT<Element>  operator+  (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
00270 template <class Element> TMatrixT<Element>  operator+  (const TMatrixT   <Element> &source ,      Element               val    );
00271 template <class Element> TMatrixT<Element>  operator+  (      Element               val    ,const TMatrixT   <Element> &source );
00272 template <class Element> TMatrixT<Element>  operator-  (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
00273 template <class Element> TMatrixT<Element>  operator-  (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
00274 template <class Element> TMatrixT<Element>  operator-  (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
00275 template <class Element> TMatrixT<Element>  operator-  (const TMatrixT   <Element> &source ,      Element               val    );
00276 template <class Element> TMatrixT<Element>  operator-  (      Element               val    ,const TMatrixT   <Element> &source );
00277 template <class Element> TMatrixT<Element>  operator*  (      Element               val    ,const TMatrixT   <Element> &source );
00278 template <class Element> TMatrixT<Element>  operator*  (const TMatrixT   <Element> &source ,      Element               val    );
00279 template <class Element> TMatrixT<Element>  operator*  (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
00280 template <class Element> TMatrixT<Element>  operator*  (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
00281 template <class Element> TMatrixT<Element>  operator*  (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
00282 template <class Element> TMatrixT<Element>  operator*  (const TMatrixTSym<Element> &source1,const TMatrixTSym<Element> &source2);
00283 // Preventing warnings with -Weffc++ in GCC since overloading the || and && operators was a design choice.
00284 #if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) >= 40600
00285 #pragma GCC diagnostic push
00286 #pragma GCC diagnostic ignored "-Weffc++"
00287 #endif
00288 template <class Element> TMatrixT<Element>  operator&& (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
00289 template <class Element> TMatrixT<Element>  operator&& (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
00290 template <class Element> TMatrixT<Element>  operator&& (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
00291 template <class Element> TMatrixT<Element>  operator|| (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
00292 template <class Element> TMatrixT<Element>  operator|| (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
00293 template <class Element> TMatrixT<Element>  operator|| (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
00294 #if (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) >= 40600
00295 #pragma GCC diagnostic pop
00296 #endif
00297 template <class Element> TMatrixT<Element>  operator>  (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
00298 template <class Element> TMatrixT<Element>  operator>  (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
00299 template <class Element> TMatrixT<Element>  operator>  (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
00300 template <class Element> TMatrixT<Element>  operator>= (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
00301 template <class Element> TMatrixT<Element>  operator>= (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
00302 template <class Element> TMatrixT<Element>  operator>= (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
00303 template <class Element> TMatrixT<Element>  operator<= (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
00304 template <class Element> TMatrixT<Element>  operator<= (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
00305 template <class Element> TMatrixT<Element>  operator<= (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
00306 template <class Element> TMatrixT<Element>  operator<  (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
00307 template <class Element> TMatrixT<Element>  operator<  (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
00308 template <class Element> TMatrixT<Element>  operator<  (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
00309 template <class Element> TMatrixT<Element>  operator!= (const TMatrixT   <Element> &source1,const TMatrixT   <Element> &source2);
00310 template <class Element> TMatrixT<Element>  operator!= (const TMatrixT   <Element> &source1,const TMatrixTSym<Element> &source2);
00311 template <class Element> TMatrixT<Element>  operator!= (const TMatrixTSym<Element> &source1,const TMatrixT   <Element> &source2);
00312 
00313 template <class Element> TMatrixT<Element> &Add        (TMatrixT<Element> &target,      Element               scalar,const TMatrixT   <Element> &source);
00314 template <class Element> TMatrixT<Element> &Add        (TMatrixT<Element> &target,      Element               scalar,const TMatrixTSym<Element> &source);
00315 template <class Element> TMatrixT<Element> &ElementMult(TMatrixT<Element> &target,const TMatrixT   <Element> &source);
00316 template <class Element> TMatrixT<Element> &ElementMult(TMatrixT<Element> &target,const TMatrixTSym<Element> &source);
00317 template <class Element> TMatrixT<Element> &ElementDiv (TMatrixT<Element> &target,const TMatrixT   <Element> &source);
00318 template <class Element> TMatrixT<Element> &ElementDiv (TMatrixT<Element> &target,const TMatrixTSym<Element> &source);
00319 
00320 template <class Element> void AMultB (const Element * const ap,Int_t na,Int_t ncolsa,
00321                                       const Element * const bp,Int_t nb,Int_t ncolsb,Element *cp);
00322 template <class Element> void AtMultB(const Element * const ap,Int_t ncolsa,
00323                                       const Element * const bp,Int_t nb,Int_t ncolsb,Element *cp);
00324 template <class Element> void AMultBt(const Element * const ap,Int_t na,Int_t ncolsa,
00325                                       const Element * const bp,Int_t nb,Int_t ncolsb,Element *cp);
00326 
00327 #endif