Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-17 09:30:48 

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,
0082 // clang++ <v20 (-Wshadow) complains about shadowing TAttMarker.h global enum EMarkerStyle. Let's silence warning:
0083 #if defined(__clang__) && __clang_major__ < 20
0084 #pragma clang diagnostic push
0085 #pragma clang diagnostic ignored "-Wshadow"
0086 #endif
0087       kPlus,
0088 #if defined(__clang__) && __clang_major__ < 20
0089 #pragma clang diagnostic pop
0090 #endif
0091       kMinus };
0092 
0093    TMatrixTSparse() { fElements = nullptr; fRowIndex = nullptr; fColIndex = nullptr; }
0094    TMatrixTSparse(Int_t nrows,Int_t ncols);
0095    TMatrixTSparse(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Int_t nr_nonzeros = 0);
0096    TMatrixTSparse(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
0097                   Int_t *row, Int_t *col,Element *data);
0098    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);
0099    TMatrixTSparse(const TMatrixTSparse<Element> &another);
0100    TMatrixTSparse(const TMatrixT<Element>       &another);
0101 
0102    TMatrixTSparse(EMatrixCreatorsOp1 op,const TMatrixTSparse<Element> &prototype);
0103    TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b);
0104    TMatrixTSparse(const TMatrixTSparse<Element> &a,EMatrixCreatorsOp2 op,const TMatrixT      <Element> &b);
0105    TMatrixTSparse(const TMatrixT      <Element> &a,EMatrixCreatorsOp2 op,const TMatrixTSparse<Element> &b);
0106 
0107    ~TMatrixTSparse() override { TMatrixTSparse::Clear(); }
0108 
0109    const Element *GetMatrixArray  () const override;
0110          Element *GetMatrixArray  () override;
0111    const Int_t    *GetRowIndexArray() const override;
0112          Int_t    *GetRowIndexArray() override;
0113    const Int_t    *GetColIndexArray() const override;
0114          Int_t    *GetColIndexArray() override;
0115 
0116    TMatrixTBase<Element>   &SetRowIndexArray(Int_t *data) override { memmove(fRowIndex,data,(this->fNrows+1)*sizeof(Int_t)); return *this; }
0117    TMatrixTBase<Element>   &SetColIndexArray(Int_t *data) override { memmove(fColIndex,data,this->fNelems*sizeof(Int_t)); return *this; }
0118 
0119            TMatrixTSparse<Element> &SetSparseIndex  (Int_t nelem_new);
0120            TMatrixTSparse<Element> &SetSparseIndex  (const TMatrixTBase<Element> &another);
0121            TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b);
0122            TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixT      <Element> &a,const TMatrixTSparse<Element> &b);
0123            TMatrixTSparse<Element> &SetSparseIndexAB(const TMatrixTSparse<Element> &a,const TMatrixT      <Element> &b)
0124                                               { return SetSparseIndexAB(b,a); }
0125 
0126    void                     GetMatrix2Array (Element *data,Option_t * /*option*/ ="") const override;
0127    TMatrixTBase<Element>   &SetMatrixArray  (const Element *data,Option_t * /*option*/="") override
0128                                                     { memcpy(fElements,data,this->fNelems*sizeof(Element)); return *this; }
0129    virtual TMatrixTBase<Element>   &SetMatrixArray  (Int_t nr_nonzeros,Int_t *irow,Int_t *icol,Element *data);
0130    virtual TMatrixTBase<Element> &
0131    SetMatrixArray(Int_t nr_nonzeros, Int_t nrows, Int_t ncols, Int_t *irow, Int_t *icol, Element *data);
0132    TMatrixTBase<Element>   &InsertRow       (Int_t row,Int_t col,const Element *v,Int_t n=-1) override;
0133    void                     ExtractRow      (Int_t row,Int_t col,      Element *v,Int_t n=-1) const override;
0134 
0135    TMatrixTBase<Element>   &ResizeTo(Int_t nrows,Int_t ncols,Int_t nr_nonzeros=-1) override;
0136    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;
0137    inline  TMatrixTBase<Element>   &ResizeTo(const TMatrixTSparse<Element> &m) {return ResizeTo(m.GetRowLwb(),m.GetRowUpb(),m.GetColLwb(),
0138                                                                                                 m.GetColUpb(),m.GetNoElements()); }
0139 
0140    void Clear(Option_t * /*option*/ ="") override { if (this->fIsOwner) {
0141                                                       if (fElements) { delete [] fElements; fElements = nullptr; }
0142                                                       if (fRowIndex) { delete [] fRowIndex; fRowIndex = nullptr; }
0143                                                       if (fColIndex) { delete [] fColIndex; fColIndex = nullptr; }
0144                                                    }
0145                                                    this->fNelems    = 0;
0146                                                    this->fNrowIndex = 0;
0147                                                  }
0148 
0149            TMatrixTSparse<Element> &Use   (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
0150                                            Int_t *pRowIndex,Int_t *pColIndex,Element *pData);
0151    const   TMatrixTSparse<Element> &Use   (Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Int_t nr_nonzeros,
0152                                            const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const
0153                                             { return (const TMatrixTSparse<Element>&)
0154                                                      ((const_cast<TMatrixTSparse<Element> *>(this))->Use(row_lwb,row_upb,col_lwb,col_upb,nr_nonzeros,
0155                                                                                              const_cast<Int_t *>(pRowIndex),
0156                                                                                              const_cast<Int_t *>(pColIndex),
0157                                                                                              const_cast<Element *>(pData))); }
0158            TMatrixTSparse<Element> &Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
0159                                            Int_t *pRowIndex,Int_t *pColIndex,Element *pData);
0160    const   TMatrixTSparse<Element> &Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
0161                                            const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const;
0162            TMatrixTSparse<Element> &Use   (TMatrixTSparse<Element> &a);
0163    const   TMatrixTSparse<Element> &Use   (const TMatrixTSparse<Element> &a) const;
0164 
0165    TMatrixTBase<Element>   &GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
0166                                             TMatrixTBase<Element> &target,Option_t *option="S") const override;
0167            TMatrixTSparse<Element>  GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,Option_t *option="S") const;
0168    TMatrixTBase<Element>   &SetSub(Int_t row_lwb,Int_t col_lwb,const TMatrixTBase<Element> &source) override;
0169 
0170    Bool_t IsSymmetric() const override { return (*this == TMatrixTSparse<Element>(kTransposed,*this)); }
0171    TMatrixTSparse<Element> &Transpose (const TMatrixTSparse<Element> &source);
0172    inline TMatrixTSparse<Element> &T () { return this->Transpose(*this); }
0173 
0174    inline void Mult(const TMatrixTSparse<Element> &a,const TMatrixTSparse<Element> &b) { AMultB(a,b,0); }
0175 
0176    TMatrixTBase<Element> &Zero       () override;
0177    TMatrixTBase<Element> &UnitMatrix () override;
0178 
0179    Element RowNorm () const override;
0180    Element ColNorm () const override;
0181    Int_t   NonZeros() const override { return this->fNelems; }
0182 
0183    TMatrixTBase<Element> &NormByDiag(const TVectorT<Element> &/*v*/,Option_t * /*option*/) override
0184                                               { MayNotUse("NormByDiag"); return *this; }
0185 
0186    // Either access a_ij as a(i,j)
0187    Element  operator()(Int_t rown,Int_t coln) const override;
0188    Element &operator()(Int_t rown,Int_t coln) override;
0189 
0190    // or as a[i][j]
0191    inline const TMatrixTSparseRow_const<Element> operator[](Int_t rown) const { return TMatrixTSparseRow_const<Element>(*this,rown); }
0192    inline       TMatrixTSparseRow      <Element> operator[](Int_t rown)       { return TMatrixTSparseRow      <Element>(*this,rown); }
0193 
0194    TMatrixTSparse<Element> &operator=(const TMatrixT<Element>       &source);
0195    TMatrixTSparse<Element> &operator=(const TMatrixTSparse<Element> &source);
0196 
0197    TMatrixTSparse<Element> &operator= (Element val);
0198    TMatrixTSparse<Element> &operator-=(Element val);
0199    TMatrixTSparse<Element> &operator+=(Element val);
0200    TMatrixTSparse<Element> &operator*=(Element val);
0201 
0202    TMatrixTSparse<Element> &operator+=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0203                                                                                 if (this == &source) APlusB (tmp,tmp,1);
0204                                                                                 else                 APlusB (tmp,source,1);
0205                                                                                 return *this; }
0206    TMatrixTSparse<Element> &operator+=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0207                                                                                 APlusB(tmp,source,1); return *this; }
0208    TMatrixTSparse<Element> &operator-=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0209                                                                                 if (this == &source) AMinusB (tmp,tmp,1);
0210                                                                                 else                 AMinusB(tmp,source,1);
0211                                                                                 return *this; }
0212    TMatrixTSparse<Element> &operator-=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0213                                                                                 AMinusB(tmp,source,1); return *this; }
0214    TMatrixTSparse<Element> &operator*=(const TMatrixTSparse<Element> &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0215                                                                                 if (this == &source) AMultB (tmp,tmp,1);
0216                                                                                 else                 AMultB (tmp,source,1);
0217                                                                                 return *this; }
0218    TMatrixTSparse<Element> &operator*=(const TMatrixT<Element>       &source) { TMatrixTSparse<Element> tmp(*this); Clear();
0219                                                                                 AMultB(tmp,source,1);
0220                                                                                 return *this; }
0221 
0222    TMatrixTBase  <Element> &Randomize  (Element alpha,Element beta,Double_t &seed) override;
0223    virtual TMatrixTSparse<Element> &RandomizePD(Element alpha,Element beta,Double_t &seed);
0224 
0225    ClassDefOverride(TMatrixTSparse,3) // Template of Sparse Matrix class
0226 };
0227 
0228 // When building with -fmodules, it instantiates all pending instantiations,
0229 // instead of delaying them until the end of the translation unit.
0230 // We 'got away with' probably because the use and the definition of the
0231 // explicit specialization do not occur in the same TU.
0232 //
0233 // In case we are building with -fmodules, we need to forward declare the
0234 // specialization in order to compile the dictionary G__Matrix.cxx.
0235 template <> TClass *TMatrixTSparse<double>::Class();
0236 
0237 template <class Element> inline const Element *TMatrixTSparse<Element>::GetMatrixArray  () const { return fElements; }
0238 template <class Element> inline       Element *TMatrixTSparse<Element>::GetMatrixArray  ()       { return fElements; }
0239 template <class Element> inline const Int_t   *TMatrixTSparse<Element>::GetRowIndexArray() const { return fRowIndex; }
0240 template <class Element> inline       Int_t   *TMatrixTSparse<Element>::GetRowIndexArray()       { return fRowIndex; }
0241 template <class Element> inline const Int_t   *TMatrixTSparse<Element>::GetColIndexArray() const { return fColIndex; }
0242 template <class Element> inline       Int_t   *TMatrixTSparse<Element>::GetColIndexArray()       { return fColIndex; }
0243 
0244 template <class Element>
0245 inline       TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
0246                                                                       Int_t *pRowIndex,Int_t *pColIndex,Element *pData)
0247                                                                         { return Use(0,nrows-1,0,ncols-1,nr_nonzeros,pRowIndex,pColIndex,pData); }
0248 template <class Element>
0249 inline const TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (Int_t nrows,Int_t ncols,Int_t nr_nonzeros,
0250                                                                       const Int_t *pRowIndex,const Int_t *pColIndex,const Element *pData) const
0251                                                                         { return Use(0,nrows-1,0,ncols-1,nr_nonzeros,pRowIndex,pColIndex,pData); }
0252 template <class Element>
0253 inline       TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (TMatrixTSparse<Element> &a)
0254                                                                         { R__ASSERT(a.IsValid());
0255                                                                            return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetColLwb(),a.GetColUpb(),
0256                                                                                       a.GetNoElements(),a.GetRowIndexArray(),
0257                                                                                       a.GetColIndexArray(),a.GetMatrixArray()); }
0258 template <class Element>
0259 inline const TMatrixTSparse<Element> &TMatrixTSparse<Element>::Use   (const TMatrixTSparse<Element> &a) const
0260                                                                         { R__ASSERT(a.IsValid());
0261                                                                            return Use(a.GetRowLwb(),a.GetRowUpb(),a.GetColLwb(),a.GetColUpb(),
0262                                                                                       a.GetNoElements(),a.GetRowIndexArray(),
0263                                                                                       a.GetColIndexArray(),a.GetMatrixArray()); }
0264 
0265 template <class Element>
0266 inline       TMatrixTSparse<Element>  TMatrixTSparse<Element>::GetSub(Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb,
0267                                                                       Option_t *option) const
0268                                                                         {
0269                                                                           TMatrixTSparse<Element> tmp;
0270                                                                           this->GetSub(row_lwb,row_upb,col_lwb,col_upb,tmp,option);
0271                                                                           return tmp;
0272                                                                         }
0273 
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+ (const TMatrixTSparse<Element> &source ,      Element                  val    );
0278 template <class Element> TMatrixTSparse<Element>  operator+ (      Element                  val    ,const TMatrixTSparse<Element> &source );
0279 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
0280 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
0281 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
0282 template <class Element> TMatrixTSparse<Element>  operator- (const TMatrixTSparse<Element> &source ,      Element                  val    );
0283 template <class Element> TMatrixTSparse<Element>  operator- (      Element                  val    ,const TMatrixTSparse<Element> &source );
0284 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source1,const TMatrixTSparse<Element> &source2);
0285 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source1,const TMatrixT<Element>       &source2);
0286 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixT<Element>       &source1,const TMatrixTSparse<Element> &source2);
0287 template <class Element> TMatrixTSparse<Element>  operator* (      Element                  val    ,const TMatrixTSparse<Element> &source );
0288 template <class Element> TMatrixTSparse<Element>  operator* (const TMatrixTSparse<Element> &source,       Element                  val    );
0289 
0290 template <class Element> TMatrixTSparse<Element> &Add        (TMatrixTSparse<Element> &target,      Element                   scalar,
0291                                                               const TMatrixTSparse<Element> &source);
0292 template <class Element> TMatrixTSparse<Element> &ElementMult(TMatrixTSparse<Element> &target,const TMatrixTSparse<Element>  &source);
0293 template <class Element> TMatrixTSparse<Element> &ElementDiv (TMatrixTSparse<Element> &target,const TMatrixTSparse<Element>  &source);
0294 
0295 template <class Element> Bool_t AreCompatible(const TMatrixTSparse<Element> &m1,const TMatrixTSparse<Element> &m2,Int_t verbose=0);
0296 
0297 #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!