Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-11-04 10:28:01 

0001 // @(#)root/matrix:$Id$
0002 // Authors: Fons Rademakers, Eddy Offermann   Feb 2004
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_TMatrixTSparse
0013 #define ROOT_TMatrixTSparse
0014 
0015 #include "TMatrixTBase.h"
0016 #include "TMatrixTUtils.h"
0017 
0018 #include <cstring>
0019 
0020 #ifdef CBLAS
0021 #include <vecLib/vBLAS.h>
0022 //#include <cblas.h>
0023 #endif
0024 
0025 //////////////////////////////////////////////////////////////////////////
0026 //                                                                      //
0027 // TMatrixTSparse                                                       //
0028 //                                                                      //
0029 // Template class of a general sparse matrix in the Harwell-Boeing      //
0030 // format                                                               //
0031 //                                                                      //
0032 //////////////////////////////////////////////////////////////////////////
0033 
0034 template<class Element> class TMatrixT;
0035 
0036 template<class Element> class TMatrixTSparse : public TMatrixTBase<Element> {
0037 
0038 protected:
0039 
0040    Int_t   *fRowIndex;  //[fNrowIndex] row index
0041    Int_t   *fColIndex;  //[fNelems]    column index
0042    Element *fElements;  //[fNelems]
0043 
0044    void Allocate(Int_t nrows,Int_t ncols,Int_t row_lwb = 0,Int_t col_lwb = 0,
0045                  Int_t init = 0,Int_t nr_nonzeros = 0);
0046 
0047   // Elementary constructors
0048    void AMultB (const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=0) {
0049       conservative_sparse_sparse_product_impl(b, a, constr);
0050    }
0051    void AMultB (const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0) {
0052                 const TMatrixTSparse<Element> bsp = b;
0053                 const TMatrixTSparse<Element> bt(TMatrixTSparse::kTransposed,bsp); AMultBt(a,bt,constr); }
0054    void AMultB (const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0) {
0055                 const TMatrixTSparse<Element> bt(TMatrixTSparse::kTransposed,b); AMultBt(a,bt,constr); }
0056 
0057    void AMultBt(const TMatrixTSparse<Element> &a, const TMatrixTSparse<Element> &b, Int_t constr = 0)
0058    {
0059       const TMatrixTSparse<Element> bt(TMatrixTSparse::kTransposed, b);
0060       conservative_sparse_sparse_product_impl(bt, a, constr);
0061    }
0062    void AMultBt(const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0);
0063    void AMultBt(const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
0064 
0065    void APlusB (const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
0066    void APlusB (const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0);
0067    void APlusB (const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0) { APlusB(b,a,constr); }
0068 
0069    void AMinusB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
0070    void AMinusB(const TMatrixTSparse<Element> &a,const TMatrixT<Element>       &b,Int_t constr=0);
0071    void AMinusB(const TMatrixT<Element>       &a,const TMatrixTSparse<Element> &b,Int_t constr=0);
0072 
0073    void conservative_sparse_sparse_product_impl(const TMatrixTSparse<Element> &lhs, const TMatrixTSparse<Element> &rhs,
0074                                                 Int_t constr = 0);
0075 
0076    Int_t ReduceSparseMatrix(Int_t nr, Int_t *row, Int_t *col, Element *data);
0077 
0078 public:
0079 
0080    enum EMatrixCreatorsOp1 { kZero,kUnit,kTransposed,kAtA };
0081    enum EMatrixCreatorsOp2 { kMult,kMultTranspose,kPlus,kMinus };
0082 
0083    TMatrixTSparse() { fElements = nullptr; fRowIndex = nullptr; fColIndex = nullptr; }
0084    TMatrixTSparse(Int_t nrows,Int_t ncols);
0085    TMatrixTSparse(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Int_t nr_nonzeros = 0);
0086    TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
0087                   Int_t *row, Int_t *col,Element *data);
0088    TMatrixTSparse(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Int_t *rowptr, Int_t *col, Element *data);
0089    TMatrixTSparse(const TMatrixTSparse<Element> &another);
0090    TMatrixTSparse(const TMatrixT<Element>       &another);
0091 
0092    TMatrixTSparse(EMatrixCreatorsOp1 op,const TMatrixTSparse<Element> &prototype);
0093    TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b);
0094    TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixT      <Element> &b);
0095    TMatrixTSparse(const TMatrixT      <Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b);
0096 
0097    ~TMatrixTSparse() override { TMatrixTSparse::Clear(); }
0098 
0099    const Element *GetMatrixArray  () const override;
0100          Element *GetMatrixArray  () override;
0101    const Int_t    *GetRowIndexArray() const override;
0102          Int_t    *GetRowIndexArray() override;
0103    const Int_t    *GetColIndexArray() const override;
0104          Int_t    *GetColIndexArray() override;
0105 
0106    TMatrixTBase<Element>   &SetRowIndexArray(Int_t *data) override { memmove(fRowIndex,data,(this->fNrows+1)*sizeof(Int_t)); return *this; }
0107    TMatrixTBase<Element>   &SetColIndexArray(Int_t *data) override { memmove(fColIndex,data,this->fNelems*sizeof(Int_t)); return *this; }
0108 
0109            TMatrixTSparse<Element> &SetSparseIndex  (Int_t nelem_new);
0110            TMatrixTSparse<Element> &SetSparseIndex  (const TMatrixTBase<Element> &another);
0111            TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b);
0112            TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixT      <Element> &a,const TMatrixTSparse<Element> &b);
0113            TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixTSparse<Element> &a,const TMatrixT      <Element> &b)
0114                                               { return SetSparseIndexAB(b,a); }
0115 
0116    void                     GetMatrix2Array (Element *data,Option_t * /*option*/ ="") const override;
0117    TMatrixTBase<Element>   &SetMatrixArray  (const Element *data,Option_t * /*option*/="") override
0118                                                     { memcpy(fElements,data,this->fNelems*sizeof(Element)); return *this; }
0119    virtual TMatrixTBase<Element>   &SetMatrixArray  (Int_t nr_nonzeros,Int_t *irow,Int_t *icol,Element *data);
0120    virtual TMatrixTBase<Element> &
0121    SetMatrixArray(Int_t nr_nonzeros, Int_t nrows, Int_t ncols, Int_t *irow, Int_t *icol, Element *data);
0122    TMatrixTBase<Element>   &InsertRow       (Int_t row,Int_t col,const Element *v,Int_t n=-1) override;
0123    void                     ExtractRow      (Int_t row,Int_t col,      Element *v,Int_t n=-1) const override;
0124 
0125    TMatrixTBase<Element>   &ResizeTo(Int_t nrows,Int_t ncols,Int_t nr_nonzeros=-1) override;
0126    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;
0127    inline  TMatrixTBase<Element>   &ResizeTo(const TMatrixTSparse<Element> &m) {return ResizeTo(m.GetRowLwb(),m.GetRowUpb(),m.GetColLwb(),
0128                                                                                                 m.GetColUpb(),m.GetNoElements()); }
0129 
0130    void Clear(Option_t * /*option*/ ="") override { if (this->fIsOwner) {
0131                                                       if (fElements) { delete [] fElements; fElements = nullptr; }
0132                                                       if (fRowIndex) { delete [] fRowIndex; fRowIndex = nullptr; }
0133                                                       if (fColIndex) { delete [] fColIndex; fColIndex = nullptr; }
0134                                                    }
0135                                                    this->fNelems    = 0;
0136                                                    this->fNrowIndex = 0;
0137                                                  }
0138 
0139            TMatrixTSparse<Element> &Use   (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
0140                                            Int_t *pRowIndex,Int_t *pColIndex,Element *pData);
0141    const   TMatrixTSparse<Element> &Use   (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
0142                                            const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const
0143                                             { return (const TMatrixTSparse<Element>&)
0144                                                      ((const_cast<TMatrixTSparse<Element> *>(this))->Use(row_lwb,row_upb,col_lwb,col_upb,nr_nonzeros,
0145                                                                                              const_cast<Int_t *>(pRowIndex),
0146                                                                                              const_cast<Int_t *>(pColIndex),
0147                                                                                              const_cast<Element *>(pData))); }
0148            TMatrixTSparse<Element> &Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
0149                                            Int_t *pRowIndex,Int_t *pColIndex,Element *pData);
0150    const   TMatrixTSparse<Element> &Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
0151                                            const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const;
0152            TMatrixTSparse<Element> &Use   (TMatrixTSparse<Element> &a);
0153    const   TMatrixTSparse<Element> &Use   (const TMatrixTSparse<Element> &a) const;
0154 
0155    TMatrixTBase<Element>   &GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
0156                                             TMatrixTBase<Element> &target,Option_t *option="S") const override;
0157            TMatrixTSparse<Element>  GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Option_t *option="S") const;
0158    TMatrixTBase<Element>   &SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source) override;
0159 
0160    Bool_t IsSymmetric() const override { return (*this == TMatrixTSparse<Element>(kTransposed,*this)); }
0161    TMatrixTSparse<Element> &Transpose (const TMatrixTSparse<Element> &source);
0162    inline TMatrixTSparse<Element> &T () { return this->Transpose(*this); }
0163 
0164    inline void Mult(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b) { AMultB(a,b,0); }
0165 
0166    TMatrixTBase<Element> &Zero       () override;
0167    TMatrixTBase<Element> &UnitMatrix () override;
0168 
0169    Element RowNorm () const override;
0170    Element ColNorm () const override;
0171    Int_t   NonZeros() const override { return this->fNelems; }
0172 
0173    TMatrixTBase<Element> &NormByDiag(const TVectorT<Element> &/*v*/,Option_t * /*option*/) override
0174                                               { MayNotUse("NormByDiag"); return *this; }
0175 
0176    // Either access a_ij as a(i,j)
0177    Element  operator()(Int_t rown,Int_t coln) const override;
0178    Element &operator()(Int_t rown,Int_t coln) override;
0179 
0180    // or as a[i][j]
0181    inline const TMatrixTSparseRow_const<Element> operator[](Int_t rown) const { return TMatrixTSparseRow_const<Element>(*this,rown); }
0182    inline       TMatrixTSparseRow      <Element> operator[](Int_t rown)       { return TMatrixTSparseRow      <Element>(*this,rown); }
0183 
0184    TMatrixTSparse<Element> &operator=(const TMatrixT<Element>       &source);
0185    TMatrixTSparse<Element> &operator=(const TMatrixTSparse<Element> &source);
0186 
0187    TMatrixTSparse<Element> &operator= (Element val);
0188    TMatrixTSparse<Element> &operator-=(Element val);
0189    TMatrixTSparse<Element> &operator+=(Element val);
0190    TMatrixTSparse<Element> &operator*=(Element val);
0191 
0192    TMatrixTSparse<Element> &operator+=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0193                                                                                 if (this == &source) APlusB (tmp,tmp,1);
0194                                                                                 else                 APlusB (tmp,source,1);
0195                                                                                 return *this; }
0196    TMatrixTSparse<Element> &operator+=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0197                                                                                 APlusB(tmp,source,1); return *this; }
0198    TMatrixTSparse<Element> &operator-=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0199                                                                                 if (this == &source) AMinusB (tmp,tmp,1);
0200                                                                                 else                 AMinusB(tmp,source,1);
0201                                                                                 return *this; }
0202    TMatrixTSparse<Element> &operator-=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0203                                                                                 AMinusB(tmp,source,1); return *this; }
0204    TMatrixTSparse<Element> &operator*=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0205                                                                                 if (this == &source) AMultB (tmp,tmp,1);
0206                                                                                 else                 AMultB (tmp,source,1);
0207                                                                                 return *this; }
0208    TMatrixTSparse<Element> &operator*=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0209                                                                                 AMultB(tmp,source,1);
0210                                                                                 return *this; }
0211 
0212    TMatrixTBase  <Element> &Randomize  (Element alpha,Element beta,Double_t &seed) override;
0213    virtual TMatrixTSparse<Element> &RandomizePD(Element alpha,Element beta,Double_t &seed);
0214 
0215    ClassDefOverride(TMatrixTSparse,3) // Template of Sparse Matrix class
0216 };
0217 
0218 // When building with -fmodules, it instantiates all pending instantiations,
0219 // instead of delaying them until the end of the translation unit.
0220 // We 'got away with' probably because the use and the definition of the
0221 // explicit specialization do not occur in the same TU.
0222 //
0223 // In case we are building with -fmodules, we need to forward declare the
0224 // specialization in order to compile the dictionary G__Matrix.cxx.
0225 template <> TClass *TMatrixTSparse<double>::Class();
0226 
0227 template <class Element> inline const Element *TMatrixTSparse<Element>::GetMatrixArray  () const { return fElements; }
0228 template <class Element> inline       Element *TMatrixTSparse<Element>::GetMatrixArray  ()       { return fElements; }
0229 template <class Element> inline const Int_t   *TMatrixTSparse<Element>::GetRowIndexArray() const { return fRowIndex; }
0230 template <class Element> inline       Int_t   *TMatrixTSparse<Element>::GetRowIndexArray()       { return fRowIndex; }
0231 template <class Element> inline const Int_t   *TMatrixTSparse<Element>::GetColIndexArray() const { return fColIndex; }
0232 template <class Element> inline       Int_t   *TMatrixTSparse<Element>::GetColIndexArray()       { return fColIndex; }
0233 
0234 template <class Element>
0235 inline       TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
0236                                                                       Int_t *pRowIndex,Int_t *pColIndex,Element *pData)
0237                                                                         { return Use(0,nrows-1,0,ncols-1,nr_nonzeros,pRowIndex,pColIndex,pData); }
0238 template <class Element>
0239 inline const TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
0240                                                                       const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const
0241                                                                         { return Use(0,nrows-1,0,ncols-1,nr_nonzeros,pRowIndex,pColIndex,pData); }
0242 template <class Element>
0243 inline       TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (TMatrixTSparse<Element> &a)
0244                                                                         { R__ASSERT(a.IsValid());
0245                                                                            return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetColLwb(),a.GetColUpb(),
0246                                                                                       a.GetNoElements(),a.GetRowIndexArray(),
0247                                                                                       a.GetColIndexArray(),a.GetMatrixArray()); }
0248 template <class Element>
0249 inline const TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (const TMatrixTSparse<Element> &a) const
0250                                                                         { R__ASSERT(a.IsValid());
0251                                                                            return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetColLwb(),a.GetColUpb(),
0252                                                                                       a.GetNoElements(),a.GetRowIndexArray(),
0253                                                                                       a.GetColIndexArray(),a.GetMatrixArray()); }
0254 
0255 template <class Element>
0256 inline       TMatrixTSparse<Element>  TMatrixTSparse<Element>::GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
0257                                                                       Option_t *option) const
0258                                                                         {
0259                                                                           TMatrixTSparse<Element> tmp;
0260                                                                           this->GetSub(row_lwb,row_upb,col_lwb,col_upb,tmp,option);
0261                                                                           return tmp;
0262                                                                         }
0263 
0264 template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
0265 template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
0266 template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
0267 template <class Element> TMatrixTSparse<Element>  operator+ (const TMatrixTSparse<Element> &source ,      Element                  val    );
0268 template <class Element> TMatrixTSparse<Element>  operator+ (      Element                  val    ,const TMatrixTSparse<Element> &source );
0269 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
0270 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
0271 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
0272 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source ,      Element                  val    );
0273 template <class Element> TMatrixTSparse<Element>  operator- (      Element                  val    ,const TMatrixTSparse<Element> &source );
0274 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
0275 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
0276 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
0277 template <class Element> TMatrixTSparse<Element>  operator* (      Element                  val    ,const TMatrixTSparse<Element> &source );
0278 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source,       Element                  val    );
0279 
0280 template <class Element> TMatrixTSparse<Element> &Add        (TMatrixTSparse<Element> &target,      Element                   scalar,
0281                                                               const TMatrixTSparse<Element> &source);
0282 template <class Element> TMatrixTSparse<Element> &ElementMult(TMatrixTSparse<Element> &target,const TMatrixTSparse<Element>  &source);
0283 template <class Element> TMatrixTSparse<Element> &ElementDiv (TMatrixTSparse<Element> &target,const TMatrixTSparse<Element>  &source);
0284 
0285 template <class Element> Bool_t AreCompatible(const TMatrixTSparse<Element> &m1,const TMatrixTSparse<Element> &m2,Int_t verbose=0);
0286 
0287 #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!