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

00001 // @(#)root/matrix:$Id: TMatrixTSparse.h 32616 2010-03-15 16:56:42Z rdm $
00002 // Authors: Fons Rademakers, Eddy Offermann   Feb 2004
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_TMatrixTSparse
00013 #define ROOT_TMatrixTSparse
00014 
00015 #ifndef ROOT_TMatrixTBase
00016 #include "TMatrixTBase.h"
00017 #endif
00018 #ifndef ROOT_TMatrixTUtils
00019 #include "TMatrixTUtils.h"
00020 #endif
00021 
00022 
00023 #ifdef CBLAS
00024 #include <vecLib/vBLAS.h>
00025 //#include <cblas.h>
00026 #endif
00027 
00028 //////////////////////////////////////////////////////////////////////////
00029 //                                                                      //
00030 // TMatrixTSparse                                                       //
00031 //                                                                      //
00032 // Template class of a general sparse matrix in the Harwell-Boeing      //
00033 // format                                                               //
00034 //                                                                      //
00035 //////////////////////////////////////////////////////////////////////////
00036 
00037 template<class Element> class TMatrixT;
00038 
00039 template<class Element> class TMatrixTSparse : public TMatrixTBase<Element> {
00040 
00041 protected:
00042 
00043    Int_t   *fRowIndex;  //[fNrowIndex] row index
00044    Int_t   *fColIndex;  //[fNelems]    column index
00045    Element *fElements;  //[fNelems]
00046 
00047    void Allocate(Int_t nrows,Int_t ncols,Int_t row_lwb = 0,Int_t col_lwb = 0,
00048                  Int_t init = 0,Int_t nr_nonzeros = 0);
00049 
00050   // Elementary constructors
00051    void AMultB (const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=0) {
00052                 const TMatrixTSparse<Element> bt(TMatrixTSparse::kTransposed,b); AMultBt(a,bt,constr); }
00053    void AMultB (const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0) {
00054                 const TMatrixTSparse<Element> bsp = b;
00055                 const TMatrixTSparse<Element> bt(TMatrixTSparse::kTransposed,bsp); AMultBt(a,bt,constr); }
00056    void AMultB (const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0) {
00057                 const TMatrixTSparse<Element> bt(TMatrixTSparse::kTransposed,b); AMultBt(a,bt,constr); }
00058 
00059    void AMultBt(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
00060    void AMultBt(const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0);
00061    void AMultBt(const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
00062 
00063    void APlusB (const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
00064    void APlusB (const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0);
00065    void APlusB (const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0) { APlusB(b,a,constr); }
00066 
00067    void AMinusB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
00068    void AMinusB(const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0);
00069    void AMinusB(const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
00070 
00071 public:
00072 
00073    enum EMatrixCreatorsOp1 { kZero,kUnit,kTransposed,kAtA };
00074    enum EMatrixCreatorsOp2 { kMult,kMultTranspose,kPlus,kMinus };
00075 
00076    TMatrixTSparse() { fElements = 0; fRowIndex = 0; fColIndex = 0; }
00077    TMatrixTSparse(Int_t nrows,Int_t ncols);
00078    TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
00079    TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
00080                   Int_t *row, Int_t *col,Element *data);
00081    TMatrixTSparse(const TMatrixTSparse<Element> &another);
00082    TMatrixTSparse(const TMatrixT<Element>       &another);
00083 
00084    TMatrixTSparse(EMatrixCreatorsOp1 op,const TMatrixTSparse<Element> &prototype);
00085    TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b);
00086    TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixT      <Element> &b);
00087    TMatrixTSparse(const TMatrixT      <Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b);
00088 
00089    virtual ~TMatrixTSparse() { Clear(); }
00090 
00091    virtual const Element *GetMatrixArray  () const;
00092    virtual       Element *GetMatrixArray  ();
00093    virtual const Int_t    *GetRowIndexArray() const;
00094    virtual       Int_t    *GetRowIndexArray();
00095    virtual const Int_t    *GetColIndexArray() const;
00096    virtual       Int_t    *GetColIndexArray();
00097 
00098    virtual TMatrixTBase<Element>   &SetRowIndexArray(Int_t *data) { memmove(fRowIndex,data,(this->fNrows+1)*sizeof(Int_t)); return *this; }
00099    virtual TMatrixTBase<Element>   &SetColIndexArray(Int_t *data) { memmove(fColIndex,data,this->fNelems*sizeof(Int_t)); return *this; }
00100 
00101            TMatrixTSparse<Element> &SetSparseIndex  (Int_t nelem_new);
00102            TMatrixTSparse<Element> &SetSparseIndex  (const TMatrixTBase<Element> &another);
00103            TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b);
00104            TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixT      <Element> &a,const TMatrixTSparse<Element> &b);
00105            TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixTSparse<Element> &a,const TMatrixT      <Element> &b)
00106                                               { return SetSparseIndexAB(b,a); }
00107 
00108    virtual void                     GetMatrix2Array (Element *data,Option_t * /*option*/ ="") const;
00109    virtual TMatrixTBase<Element>   &SetMatrixArray  (const Element *data,Option_t * /*option*/="")
00110                                                     { memcpy(fElements,data,this->fNelems*sizeof(Element)); return *this; }
00111    virtual TMatrixTBase<Element>   &SetMatrixArray  (Int_t nr_nonzeros,Int_t *irow,Int_t *icol,Element *data);
00112    virtual TMatrixTBase<Element>   &InsertRow       (Int_t row,Int_t col,const Element *v,Int_t n=-1);
00113    virtual void                     ExtractRow      (Int_t row,Int_t col,      Element *v,Int_t n=-1) const;
00114 
00115    virtual TMatrixTBase<Element>   &ResizeTo(Int_t nrows,Int_t ncols,Int_t nr_nonzeros=-1);
00116    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);
00117    inline  TMatrixTBase<Element>   &ResizeTo(const TMatrixTSparse<Element> &m) {return ResizeTo(m.GetRowLwb(),m.GetRowUpb(),m.GetColLwb(),
00118                                                                                                 m.GetColUpb(),m.GetNoElements()); }
00119 
00120    virtual void Clear(Option_t * /*option*/ ="") { if (this->fIsOwner) {
00121                                                      if (fElements) delete [] fElements; fElements = 0;
00122                                                      if (fRowIndex) delete [] fRowIndex; fRowIndex = 0;
00123                                                      if (fColIndex) delete [] fColIndex; fColIndex = 0;
00124                                                    }
00125                                                    this->fNelems    = 0;
00126                                                    this->fNrowIndex = 0;
00127                                                  }
00128 
00129            TMatrixTSparse<Element> &Use   (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
00130                                            Int_t *pRowIndex,Int_t *pColIndex,Element *pData);
00131    const   TMatrixTSparse<Element> &Use   (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
00132                                            const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const
00133                                             { return (const TMatrixTSparse<Element>&)
00134                                                      ((const_cast<TMatrixTSparse<Element> *>(this))->Use(row_lwb,row_upb,col_lwb,col_upb,nr_nonzeros,
00135                                                                                              const_cast<Int_t *>(pRowIndex),
00136                                                                                              const_cast<Int_t *>(pColIndex),
00137                                                                                              const_cast<Element *>(pData))); }
00138            TMatrixTSparse<Element> &Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
00139                                            Int_t *pRowIndex,Int_t *pColIndex,Element *pData);
00140    const   TMatrixTSparse<Element> &Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
00141                                            const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const;
00142            TMatrixTSparse<Element> &Use   (TMatrixTSparse<Element> &a);
00143    const   TMatrixTSparse<Element> &Use   (const TMatrixTSparse<Element> &a) const;
00144 
00145    virtual TMatrixTBase<Element>   &GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
00146                                             TMatrixTBase<Element> &target,Option_t *option="S") const;
00147            TMatrixTSparse<Element>  GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Option_t *option="S") const;
00148    virtual TMatrixTBase<Element>   &SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source);
00149 
00150    virtual Bool_t IsSymmetric() const { return (*this == TMatrixTSparse<Element>(kTransposed,*this)); }
00151    TMatrixTSparse<Element> &Transpose (const TMatrixTSparse<Element> &source);
00152    inline TMatrixTSparse<Element> &T () { return this->Transpose(*this); }
00153 
00154    inline void Mult(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b) { AMultB(a,b,0); }
00155 
00156    virtual TMatrixTBase<Element> &Zero       ();
00157    virtual TMatrixTBase<Element> &UnitMatrix ();
00158 
00159    virtual Element RowNorm () const;
00160    virtual Element ColNorm () const;
00161    virtual Int_t   NonZeros() const { return this->fNelems; }
00162 
00163    virtual TMatrixTBase<Element> &NormByDiag(const TVectorT<Element> &/*v*/,Option_t * /*option*/)
00164                                               { MayNotUse("NormByDiag"); return *this; }
00165 
00166    // Either access a_ij as a(i,j)
00167    Element  operator()(Int_t rown,Int_t coln) const;
00168    Element &operator()(Int_t rown,Int_t coln);
00169 
00170    // or as a[i][j]
00171    inline const TMatrixTSparseRow_const<Element> operator[](Int_t rown) const { return TMatrixTSparseRow_const<Element>(*this,rown); }
00172    inline       TMatrixTSparseRow      <Element> operator[](Int_t rown)       { return TMatrixTSparseRow      <Element>(*this,rown); }
00173 
00174    TMatrixTSparse<Element> &operator=(const TMatrixT<Element>       &source);
00175    TMatrixTSparse<Element> &operator=(const TMatrixTSparse<Element> &source);
00176 
00177    TMatrixTSparse<Element> &operator= (Element val);
00178    TMatrixTSparse<Element> &operator-=(Element val);
00179    TMatrixTSparse<Element> &operator+=(Element val);
00180    TMatrixTSparse<Element> &operator*=(Element val);
00181 
00182    TMatrixTSparse<Element> &operator+=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
00183                                                                                 if (this == &source) APlusB (tmp,tmp,1);
00184                                                                                 else                 APlusB (tmp,source,1); return *this; }
00185    TMatrixTSparse<Element> &operator+=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this); Clear();
00186                                                                                 APlusB(tmp,source,1); return *this; }
00187    TMatrixTSparse<Element> &operator-=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
00188                                                                                 if (this == &source) AMinusB (tmp,tmp,1);
00189                                                                                 else                 AMinusB(tmp,source,1); return *this; }
00190    TMatrixTSparse<Element> &operator-=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this); Clear();
00191                                                                                 AMinusB(tmp,source,1); return *this; }
00192    TMatrixTSparse<Element> &operator*=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
00193                                                                                 if (this == &source) AMultB (tmp,tmp,1);
00194                                                                                 else                 AMultB (tmp,source,1); return *this; }
00195    TMatrixTSparse<Element> &operator*=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this); Clear();
00196                                                                                 AMultB(tmp,source,1); return *this; }
00197 
00198    virtual TMatrixTBase  <Element> &Randomize  (Element alpha,Element beta,Double_t &seed);
00199    virtual TMatrixTSparse<Element> &RandomizePD(Element alpha,Element beta,Double_t &seed);
00200 
00201    ClassDef(TMatrixTSparse,3) // Template of Sparse Matrix class
00202 };
00203 
00204 template <class Element> inline const Element *TMatrixTSparse<Element>::GetMatrixArray  () const { return fElements; }
00205 template <class Element> inline       Element *TMatrixTSparse<Element>::GetMatrixArray  ()       { return fElements; }
00206 template <class Element> inline const Int_t   *TMatrixTSparse<Element>::GetRowIndexArray() const { return fRowIndex; }
00207 template <class Element> inline       Int_t   *TMatrixTSparse<Element>::GetRowIndexArray()       { return fRowIndex; }
00208 template <class Element> inline const Int_t   *TMatrixTSparse<Element>::GetColIndexArray() const { return fColIndex; }
00209 template <class Element> inline       Int_t   *TMatrixTSparse<Element>::GetColIndexArray()       { return fColIndex; }
00210 
00211 template <class Element>
00212 inline       TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
00213                                                                       Int_t *pRowIndex,Int_t *pColIndex,Element *pData)
00214                                                                         { return Use(0,nrows-1,0,ncols-1,nr_nonzeros,pRowIndex,pColIndex,pData); }
00215 template <class Element>
00216 inline const TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
00217                                                                       const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const
00218                                                                         { return Use(0,nrows-1,0,ncols-1,nr_nonzeros,pRowIndex,pColIndex,pData); }
00219 template <class Element>
00220 inline       TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (TMatrixTSparse<Element> &a)
00221                                                                         { R__ASSERT(a.IsValid());
00222                                                                            return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetColLwb(),a.GetColUpb(),
00223                                                                                       a.GetNoElements(),a.GetRowIndexArray(),
00224                                                                                       a.GetColIndexArray(),a.GetMatrixArray()); }
00225 template <class Element>
00226 inline const TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (const TMatrixTSparse<Element> &a) const
00227                                                                         { R__ASSERT(a.IsValid());
00228                                                                            return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetColLwb(),a.GetColUpb(),
00229                                                                                       a.GetNoElements(),a.GetRowIndexArray(),
00230                                                                                       a.GetColIndexArray(),a.GetMatrixArray()); }
00231 
00232 template <class Element>
00233 inline       TMatrixTSparse<Element>  TMatrixTSparse<Element>::GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
00234                                                                       Option_t *option) const
00235                                                                         {
00236                                                                           TMatrixTSparse<Element> tmp;
00237                                                                           this->GetSub(row_lwb,row_upb,col_lwb,col_upb,tmp,option);
00238                                                                           return tmp;
00239                                                                         }
00240 
00241 template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
00242 template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
00243 template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
00244 template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixTSparse<Element> &source ,      Element                  val    );
00245 template <class Element> TMatrixTSparse<Element>  operator+ (      Element                  val    ,const TMatrixTSparse<Element> &source );
00246 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
00247 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
00248 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
00249 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source ,      Element                  val    );
00250 template <class Element> TMatrixTSparse<Element>  operator- (      Element                  val    ,const TMatrixTSparse<Element> &source );
00251 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
00252 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
00253 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
00254 template <class Element> TMatrixTSparse<Element>  operator* (      Element                  val    ,const TMatrixTSparse<Element> &source );
00255 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source,       Element                  val    );
00256 
00257 template <class Element> TMatrixTSparse<Element> &Add        (TMatrixTSparse<Element> &target,      Element                   scalar,
00258                                                               const TMatrixTSparse<Element> &source);
00259 template <class Element> TMatrixTSparse<Element> &ElementMult(TMatrixTSparse<Element> &target,const TMatrixTSparse<Element>  &source);
00260 template <class Element> TMatrixTSparse<Element> &ElementDiv (TMatrixTSparse<Element> &target,const TMatrixTSparse<Element>  &source);
00261 
00262 template <class Element> Bool_t AreCompatible(const TMatrixTSparse<Element> &m1,const TMatrixTSparse<Element> &m2,Int_t verbose=0);
00263 
00264 #endif