Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /include/root/TMatrixTUtils.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // @(#)root/matrix:$Id$
0002 // Authors: Fons Rademakers, Eddy Offermann   Nov 2003
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_TMatrixTUtils
0013 #define ROOT_TMatrixTUtils
0014 
0015 //////////////////////////////////////////////////////////////////////////
0016 //                                                                      //
0017 // Matrix utility classes.                                              //
0018 //                                                                      //
0019 // Templates of utility classes in the Linear Algebra Package.          //
0020 // The following classes are defined here:                              //
0021 //                                                                      //
0022 // Different matrix views without copying data elements :               //
0023 //   TMatrixTRow_const        TMatrixTRow                               //
0024 //   TMatrixTColumn_const     TMatrixTColumn                            //
0025 //   TMatrixTDiag_const       TMatrixTDiag                              //
0026 //   TMatrixTFlat_const       TMatrixTFlat                              //
0027 //   TMatrixTSub_const        TMatrixTSub                               //
0028 //   TMatrixTSparseRow_const  TMatrixTSparseRow                         //
0029 //   TMatrixTSparseDiag_const TMatrixTSparseDiag                        //
0030 //                                                                      //
0031 //   TElementActionT                                                    //
0032 //   TElementPosActionT                                                 //
0033 //                                                                      //
0034 //////////////////////////////////////////////////////////////////////////
0035 
0036 #include "TMatrixTBase.h"
0037 
0038 #include <initializer_list>
0039 
0040 template<class Element> class TVectorT;
0041 template<class Element> class TMatrixT;
0042 template<class Element> class TMatrixTSym;
0043 template<class Element> class TMatrixTSparse;
0044 
0045 //////////////////////////////////////////////////////////////////////////
0046 //                                                                      //
0047 // TElementActionT                                                      //
0048 //                                                                      //
0049 // A class to do a specific operation on every vector or matrix element //
0050 // (regardless of it position) as the object is being traversed.        //
0051 // This is an abstract class. Derived classes need to implement the     //
0052 // action function Operation().                                         //
0053 //                                                                      //
0054 //////////////////////////////////////////////////////////////////////////
0055 
0056 template<class Element> class TElementActionT {
0057 
0058 friend class TMatrixTBase  <Element>;
0059 friend class TMatrixT      <Element>;
0060 friend class TMatrixTSym   <Element>;
0061 friend class TMatrixTSparse<Element>;
0062 friend class TVectorT      <Element>;
0063 
0064 protected:
0065    virtual ~TElementActionT() { }
0066    virtual void Operation(Element &element) const = 0;
0067 
0068 private:
0069    TElementActionT& operator=(const TElementActionT<Element> &) {return *this;}
0070 };
0071 
0072 //////////////////////////////////////////////////////////////////////////
0073 //                                                                      //
0074 // TElementPosActionT                                                   //
0075 //                                                                      //
0076 // A class to do a specific operation on every vector or matrix element //
0077 // as the object is being traversed. This is an abstract class.         //
0078 // Derived classes need to implement the action function Operation().   //
0079 // In the action function the location of the current element is        //
0080 // known (fI=row, fJ=columns).                                          //
0081 //                                                                      //
0082 //////////////////////////////////////////////////////////////////////////
0083 
0084 template<class Element> class TElementPosActionT {
0085 
0086 friend class TMatrixTBase  <Element>;
0087 friend class TMatrixT      <Element>;
0088 friend class TMatrixTSym   <Element>;
0089 friend class TMatrixTSparse<Element>;
0090 friend class TVectorT      <Element>;
0091 
0092 protected:
0093    mutable Int_t fI; // i position of element being passed to Operation()
0094    mutable Int_t fJ; // j position of element being passed to Operation()
0095    virtual ~TElementPosActionT() { }
0096    virtual void Operation(Element &element) const = 0;
0097 
0098 private:
0099    TElementPosActionT<Element>& operator=(const TElementPosActionT<Element> &) {return *this;}
0100 };
0101 
0102 //////////////////////////////////////////////////////////////////////////
0103 //                                                                      //
0104 // TMatrixTRow_const                                                    //
0105 //                                                                      //
0106 // Template class represents a row of a TMatrixT/TMatrixTSym            //
0107 //                                                                      //
0108 //////////////////////////////////////////////////////////////////////////
0109 
0110 template<class Element> class TMatrixTRow_const {
0111 
0112 protected:
0113    const TMatrixTBase<Element> *fMatrix;  //  the matrix I am a row of
0114          Int_t                  fRowInd;  //  effective row index
0115          Int_t                  fInc;     //  if ptr = @a[row,i], then ptr+inc = @a[row,i+1]
0116    const Element               *fPtr;     //  pointer to the a[row,0]
0117 
0118 public:
0119    TMatrixTRow_const() { fMatrix = nullptr; fRowInd = 0; fInc = 0; fPtr = nullptr; }
0120    TMatrixTRow_const(const TMatrixT   <Element> &matrix,Int_t row);
0121    TMatrixTRow_const(const TMatrixTSym<Element> &matrix,Int_t row);
0122   TMatrixTRow_const(const TMatrixTRow_const<Element>& trc):
0123     fMatrix(trc.fMatrix), fRowInd(trc.fRowInd), fInc(trc.fInc), fPtr(trc.fPtr) { }
0124   TMatrixTRow_const<Element>& operator=(const TMatrixTRow_const<Element>& trc) {
0125     if(this != &trc) { fMatrix=trc.fMatrix; fRowInd=trc.fRowInd; fInc=trc.fInc; fPtr=trc.fPtr; } return *this;}
0126    virtual ~TMatrixTRow_const() { }
0127 
0128    inline const TMatrixTBase<Element> *GetMatrix  () const { return fMatrix; }
0129    inline       Int_t                  GetRowIndex() const { return fRowInd; }
0130    inline       Int_t                  GetInc     () const { return fInc; }
0131    inline const Element               *GetPtr     () const { return fPtr; }
0132    inline const Element               &operator   ()(Int_t i) const {
0133       if (!fMatrix) return TMatrixTBase<Element>::NaNValue();
0134       R__ASSERT(fMatrix->IsValid());
0135       const Int_t acoln = i-fMatrix->GetColLwb();
0136       if (acoln < fMatrix->GetNcols() && acoln >= 0)
0137          return fPtr[acoln];
0138       else {
0139          Error("operator()","Request col(%d) outside matrix range of %d - %d",
0140                             i,fMatrix->GetColLwb(),fMatrix->GetColLwb()+fMatrix->GetNcols());
0141          return TMatrixTBase<Element>::NaNValue();
0142       }
0143    }
0144    inline const Element               &operator   [](Int_t i) const { return (*(const TMatrixTRow_const<Element> *)this)(i); }
0145 
0146    ClassDef(TMatrixTRow_const,0)  // Template of General Matrix Row Access class
0147 };
0148 
0149 template<class Element> class TMatrixTRow : public TMatrixTRow_const<Element> {
0150 
0151 public:
0152    TMatrixTRow() {}
0153    TMatrixTRow(TMatrixT   <Element> &matrix,Int_t row);
0154    TMatrixTRow(TMatrixTSym<Element> &matrix,Int_t row);
0155    TMatrixTRow(const TMatrixTRow<Element> &mr);
0156 
0157    inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
0158 
0159    inline const Element &operator()(Int_t i) const {
0160       if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
0161       R__ASSERT(this->fMatrix->IsValid());
0162       const Int_t acoln = i-this->fMatrix->GetColLwb();
0163       if (acoln < this->fMatrix->GetNcols() || acoln >= 0)
0164          return (this->fPtr)[acoln];
0165       else {
0166          Error("operator()","Request col(%d) outside matrix range of %d - %d",
0167                             i,this->fMatrix->GetColLwb(),this->fMatrix->GetColLwb()+this->fMatrix->GetNcols());
0168          return TMatrixTBase<Element>::NaNValue();
0169      }
0170    }
0171    inline       Element &operator()(Int_t i) {
0172       if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
0173       R__ASSERT(this->fMatrix->IsValid());
0174       const Int_t acoln = i-this->fMatrix->GetColLwb();
0175       if (acoln < this->fMatrix->GetNcols() && acoln >= 0)
0176          return (const_cast<Element *>(this->fPtr))[acoln];
0177       else {
0178          Error("operator()","Request col(%d) outside matrix range of %d - %d",
0179                             i,this->fMatrix->GetColLwb(),this->fMatrix->GetColLwb()+this->fMatrix->GetNcols());
0180          //return (const_cast<Element *>(this->fPtr))[0];
0181          return TMatrixTBase<Element>::NaNValue();
0182       }
0183    }
0184    inline const Element &operator[](Int_t i) const { return (*(const TMatrixTRow<Element> *)this)(i); }
0185    inline       Element &operator[](Int_t i)       { return (*(      TMatrixTRow<Element> *)this)(i); }
0186 
0187    void Assign (Element  val);
0188    void operator= (std::initializer_list<Element>  l);
0189    void operator+=(Element val);
0190    void operator*=(Element val);
0191 
0192    void operator=(const TMatrixTRow_const<Element> &r);
0193    TMatrixTRow<Element>& operator=(const TMatrixTRow      <Element> &r) { operator=((TMatrixTRow_const<Element> &)r); return *this;}
0194    void operator=(const TVectorT         <Element> &vec);
0195 
0196    void operator+=(const TMatrixTRow_const<Element> &r);
0197    void operator*=(const TMatrixTRow_const<Element> &r);
0198 
0199    ClassDefOverride(TMatrixTRow,0)  // Template of General Matrix Row Access class
0200 };
0201 
0202 //////////////////////////////////////////////////////////////////////////
0203 //                                                                      //
0204 // TMatrixTColumn_const                                                 //
0205 //                                                                      //
0206 // Template class represents a column of a TMatrixT/TMatrixTSym         //
0207 //                                                                      //
0208 //////////////////////////////////////////////////////////////////////////
0209 
0210 template<class Element> class TMatrixTColumn_const {
0211 
0212 protected:
0213    const TMatrixTBase<Element> *fMatrix;  //  the matrix I am a column of
0214          Int_t                  fColInd;  //  effective column index
0215          Int_t                  fInc;     //  if ptr = @a[i,col], then ptr+inc = @a[i+1,col]
0216    const Element               *fPtr;     //  pointer to the a[0,col] column
0217 
0218 public:
0219    TMatrixTColumn_const() { fMatrix = nullptr; fColInd = 0; fInc = 0; fPtr = nullptr; }
0220    TMatrixTColumn_const(const TMatrixT   <Element> &matrix,Int_t col);
0221    TMatrixTColumn_const(const TMatrixTSym<Element> &matrix,Int_t col);
0222    TMatrixTColumn_const(const TMatrixTColumn_const<Element>& trc):
0223      fMatrix(trc.fMatrix), fColInd(trc.fColInd), fInc(trc.fInc), fPtr(trc.fPtr) { }
0224    TMatrixTColumn_const<Element>& operator=(const TMatrixTColumn_const<Element>& trc) {
0225      if(this != &trc) { fMatrix=trc.fMatrix; fColInd=trc.fColInd; fInc=trc.fInc; fPtr=trc.fPtr; } return *this;}
0226    virtual ~TMatrixTColumn_const() { }
0227 
0228    inline const TMatrixTBase <Element> *GetMatrix  () const { return fMatrix; }
0229    inline       Int_t                   GetColIndex() const { return fColInd; }
0230    inline       Int_t                   GetInc     () const { return fInc; }
0231    inline const Element                *GetPtr     () const { return fPtr; }
0232    inline const Element                &operator   ()(Int_t i) const {
0233       if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
0234       R__ASSERT(fMatrix->IsValid());
0235       const Int_t arown = i-fMatrix->GetRowLwb();
0236       if (arown < fMatrix->GetNrows() && arown >= 0)
0237          return fPtr[arown*fInc];
0238       else {
0239          Error("operator()","Request row(%d) outside matrix range of %d - %d",
0240                             i,fMatrix->GetRowLwb(),fMatrix->GetRowLwb()+fMatrix->GetNrows());
0241          return TMatrixTBase<Element>::NaNValue();
0242       }
0243    }
0244    inline const Element                &operator [](Int_t i) const { return (*(const TMatrixTColumn_const<Element> *)this)(i); }
0245 
0246    ClassDef(TMatrixTColumn_const,0)  // Template of General Matrix Column Access class
0247 };
0248 
0249 template<class Element> class TMatrixTColumn : public TMatrixTColumn_const<Element> {
0250 
0251 public:
0252    TMatrixTColumn() {}
0253    TMatrixTColumn(TMatrixT   <Element>&matrix,Int_t col);
0254    TMatrixTColumn(TMatrixTSym<Element>&matrix,Int_t col);
0255    TMatrixTColumn(const TMatrixTColumn <Element>&mc);
0256 
0257    inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
0258 
0259    inline const Element &operator()(Int_t i) const {
0260       if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
0261       R__ASSERT(this->fMatrix->IsValid());
0262       const Int_t arown = i-this->fMatrix->GetRowLwb();
0263       if (arown < this->fMatrix->GetNrows() && arown >= 0)
0264          return (this->fPtr)[arown*this->fInc];
0265       else {
0266          Error("operator()","Request row(%d) outside matrix range of %d - %d",
0267                             i,this->fMatrix->GetRowLwb(),this->fMatrix->GetRowLwb()+this->fMatrix->GetNrows());
0268          return TMatrixTBase<Element>::NaNValue();
0269       }
0270    }
0271    inline       Element &operator()(Int_t i) {
0272       if (!this->fMatrix) return TMatrixTBase<Element>::NaNValue();
0273       R__ASSERT(this->fMatrix->IsValid());
0274       const Int_t arown = i-this->fMatrix->GetRowLwb();
0275 
0276       if (arown < this->fMatrix->GetNrows() && arown >= 0)
0277          return (const_cast<Element *>(this->fPtr))[arown*this->fInc];
0278       else {
0279          Error("operator()","Request row(%d) outside matrix range of %d - %d",
0280                             i,this->fMatrix->GetRowLwb(),this->fMatrix->GetRowLwb()+this->fMatrix->GetNrows());
0281          return TMatrixTBase<Element>::NaNValue();
0282       }
0283    }
0284    inline const Element &operator[](Int_t i) const { return (*(const TMatrixTColumn<Element> *)this)(i); }
0285    inline       Element &operator[](Int_t i)       { return (*(      TMatrixTColumn<Element> *)this)(i); }
0286 
0287    void Assign (Element val);
0288    // keep it for backward compatibility (but it has been removed for TMatrixTRow)
0289    void operator= (Element val) { return Assign(val); }
0290    void operator= (std::initializer_list<Element>  l);
0291    void operator+=(Element val);
0292    void operator*=(Element val);
0293 
0294    void operator=(const TMatrixTColumn_const<Element> &c);
0295    TMatrixTColumn<Element>& operator=(const TMatrixTColumn <Element> &c) { operator=((TMatrixTColumn_const<Element> &)c); return *this;}
0296    void operator=(const TVectorT            <Element> &vec);
0297 
0298    void operator+=(const TMatrixTColumn_const<Element> &c);
0299    void operator*=(const TMatrixTColumn_const<Element> &c);
0300 
0301    ClassDefOverride(TMatrixTColumn,0)  // Template of General Matrix Column Access class
0302 };
0303 
0304 //////////////////////////////////////////////////////////////////////////
0305 //                                                                      //
0306 // TMatrixTDiag_const                                                   //
0307 //                                                                      //
0308 // Template class represents the diagonal of a TMatrixT/TMatrixTSym     //
0309 //                                                                      //
0310 //////////////////////////////////////////////////////////////////////////
0311 
0312 template<class Element> class TMatrixTDiag_const {
0313 
0314 protected:
0315    const TMatrixTBase<Element> *fMatrix;  //  the matrix I am the diagonal of
0316          Int_t                  fInc;     //  if ptr=@a[i,i], then ptr+inc = @a[i+1,i+1]
0317          Int_t                  fNdiag;   //  number of diag elems, min(nrows,ncols)
0318    const Element               *fPtr;     //  pointer to the a[0,0]
0319 
0320 public:
0321    TMatrixTDiag_const() { fMatrix = nullptr; fInc = 0; fNdiag = 0; fPtr = nullptr; }
0322    TMatrixTDiag_const(const TMatrixT   <Element> &matrix);
0323    TMatrixTDiag_const(const TMatrixTSym<Element> &matrix);
0324    TMatrixTDiag_const(const TMatrixTDiag_const<Element>& trc):
0325     fMatrix(trc.fMatrix), fInc(trc.fInc), fNdiag(trc.fNdiag), fPtr(trc.fPtr) { }
0326    TMatrixTDiag_const<Element>& operator=(const TMatrixTDiag_const<Element>& trc) {
0327       if(this != &trc) { fMatrix=trc.fMatrix; fInc=trc.fInc; fNdiag=trc.fNdiag; fPtr=trc.fPtr; } return *this;}
0328    virtual ~TMatrixTDiag_const() { }
0329 
0330    inline const TMatrixTBase<Element> *GetMatrix() const { return fMatrix; }
0331    inline const Element               *GetPtr   () const { return fPtr; }
0332    inline       Int_t                  GetInc   () const { return fInc; }
0333    inline const Element               &operator ()(Int_t i) const {
0334       R__ASSERT(fMatrix->IsValid());
0335       if (i < fNdiag && i >= 0)
0336          return fPtr[i*fInc];
0337       else {
0338          Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,fNdiag);
0339          return TMatrixTBase<Element>::NaNValue();
0340       }
0341    }
0342    inline const Element               &operator [](Int_t i) const { return (*(const TMatrixTDiag_const<Element> *)this)(i); }
0343 
0344    Int_t GetNdiags() const { return fNdiag; }
0345 
0346    ClassDef(TMatrixTDiag_const,0)  // Template of General Matrix Diagonal Access class
0347 };
0348 
0349 template<class Element> class TMatrixTDiag : public TMatrixTDiag_const<Element> {
0350 
0351 public:
0352    TMatrixTDiag() {}
0353    TMatrixTDiag(TMatrixT   <Element>&matrix);
0354    TMatrixTDiag(TMatrixTSym<Element>&matrix);
0355    TMatrixTDiag(const TMatrixTDiag<Element> &md);
0356 
0357    inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
0358 
0359    inline const Element &operator()(Int_t i) const {
0360       R__ASSERT(this->fMatrix->IsValid());
0361       if (i < this->fNdiag && i >= 0)
0362          return (this->fPtr)[i*this->fInc];
0363       else {
0364          Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,this->fNdiag);
0365          return TMatrixTBase<Element>::NaNValue();
0366       }
0367    }
0368    inline       Element &operator()(Int_t i) {
0369       R__ASSERT(this->fMatrix->IsValid());
0370       if (i < this->fNdiag && i >= 0)
0371          return (const_cast<Element *>(this->fPtr))[i*this->fInc];
0372       else {
0373          Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,this->fNdiag);
0374          return (const_cast<Element *>(this->fPtr))[0];
0375       }
0376    }
0377    inline const Element &operator[](Int_t i) const { return (*(const TMatrixTDiag<Element> *)this)(i); }
0378    inline       Element &operator[](Int_t i)       { return (*(      TMatrixTDiag *)this)(i); }
0379 
0380    void operator= (Element val);
0381    void operator+=(Element val);
0382    void operator*=(Element val);
0383 
0384    void operator=(const TMatrixTDiag_const<Element> &d);
0385    TMatrixTDiag<Element>& operator=(const TMatrixTDiag <Element> &d) { operator=((TMatrixTDiag_const<Element> &)d); return *this;}
0386    void operator=(const TVectorT          <Element> &vec);
0387 
0388    void operator+=(const TMatrixTDiag_const<Element> &d);
0389    void operator*=(const TMatrixTDiag_const<Element> &d);
0390 
0391    ClassDefOverride(TMatrixTDiag,0)  // Template of General Matrix Diagonal Access class
0392 };
0393 
0394 //////////////////////////////////////////////////////////////////////////
0395 //                                                                      //
0396 // TMatrixTFlat_const                                                   //
0397 //                                                                      //
0398 // Template class represents a flat TMatrixT/TMatrixTSym                //
0399 //                                                                      //
0400 //////////////////////////////////////////////////////////////////////////
0401 
0402 template<class Element> class TMatrixTFlat_const {
0403 
0404 protected:
0405    const TMatrixTBase<Element> *fMatrix;  //  the matrix I am the diagonal of
0406          Int_t                  fNelems;  //
0407    const Element               *fPtr;     //  pointer to the a[0,0]
0408 
0409 public:
0410    TMatrixTFlat_const() { fMatrix = nullptr; fNelems = 0; fPtr = nullptr; }
0411    TMatrixTFlat_const(const TMatrixT   <Element> &matrix);
0412    TMatrixTFlat_const(const TMatrixTSym<Element> &matrix);
0413    TMatrixTFlat_const(const TMatrixTFlat_const<Element>& trc):
0414      fMatrix(trc.fMatrix), fNelems(trc.fNelems), fPtr(trc.fPtr) { }
0415    TMatrixTFlat_const<Element>& operator=(const TMatrixTFlat_const<Element>& trc) {
0416       if(this != &trc) { fMatrix=trc.fMatrix; fNelems=trc.fNelems; fPtr=trc.fPtr; } return *this;}
0417    virtual ~TMatrixTFlat_const() { }
0418 
0419    inline const TMatrixTBase<Element> *GetMatrix() const { return fMatrix; }
0420    inline const Element               *GetPtr   () const { return fPtr; }
0421    inline const Element               &operator ()(Int_t i) const {
0422       R__ASSERT(fMatrix->IsValid());
0423       if (i < fNelems && i >= 0)
0424          return fPtr[i];
0425       else {
0426          Error("operator()","Request element(%d) outside matrix range of 0 - %d",i,fNelems);
0427          return TMatrixTBase<Element>::NaNValue();
0428       }
0429    }
0430    inline const Element               &operator [](Int_t i) const { return (*(const TMatrixTFlat_const<Element> *)this)(i); }
0431 
0432    ClassDef(TMatrixTFlat_const,0)  // Template of General Matrix Flat Representation class
0433 };
0434 
0435 template<class Element> class TMatrixTFlat : public TMatrixTFlat_const<Element> {
0436 
0437 public:
0438    TMatrixTFlat() {}
0439    TMatrixTFlat(TMatrixT   <Element> &matrix);
0440    TMatrixTFlat(TMatrixTSym<Element> &matrix);
0441    TMatrixTFlat(const TMatrixTFlat<Element> &mf);
0442 
0443    inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
0444 
0445    inline const Element &operator()(Int_t i) const {
0446       R__ASSERT(this->fMatrix->IsValid());
0447       if (i < this->fNelems && i >= 0)
0448          return (this->fPtr)[i];
0449       else {
0450          Error("operator()","Request element(%d) outside matrix range of 0 - %d",i,this->fNelems);
0451          return TMatrixTBase<Element>::NaNValue();
0452       }
0453    }
0454    inline       Element &operator()(Int_t i) {
0455       R__ASSERT(this->fMatrix->IsValid());
0456       if (i < this->fNelems && i >= 0)
0457          return (const_cast<Element *>(this->fPtr))[i];
0458       else {
0459          Error("operator()","Request element(%d) outside matrix range of 0 - %d",i,this->fNelems);
0460          return TMatrixTBase<Element>::NaNValue();
0461       }
0462    }
0463    inline const Element &operator[](Int_t i) const { return (*(const TMatrixTFlat<Element> *)this)(i); }
0464    inline       Element &operator[](Int_t i)       { return (*(      TMatrixTFlat<Element> *)this)(i); }
0465 
0466    void operator= (Element val);
0467    void operator+=(Element val);
0468    void operator*=(Element val);
0469 
0470    void operator=(const TMatrixTFlat_const<Element> &f);
0471    TMatrixTFlat<Element>& operator=(const TMatrixTFlat <Element> &f) { operator=((TMatrixTFlat_const<Element> &)f); return *this;}
0472    void operator=(const TVectorT          <Element> &vec);
0473 
0474    void operator+=(const TMatrixTFlat_const<Element> &f);
0475    void operator*=(const TMatrixTFlat_const<Element> &f);
0476 
0477    ClassDefOverride(TMatrixTFlat,0)  // Template of General Matrix Flat Representation class
0478 };
0479 
0480 //////////////////////////////////////////////////////////////////////////
0481 //                                                                      //
0482 // TMatrixTSub_const                                                    //
0483 //                                                                      //
0484 // Template class represents a sub matrix of TMatrixT/TMatrixTSym       //
0485 //                                                                      //
0486 //////////////////////////////////////////////////////////////////////////
0487 
0488 template<class Element> class TMatrixTSub_const {
0489 
0490 protected:
0491    const TMatrixTBase<Element> *fMatrix;    //  the matrix I am a submatrix of
0492          Int_t                  fRowOff;    //
0493          Int_t                  fColOff;    //
0494          Int_t                  fNrowsSub;  //
0495          Int_t                  fNcolsSub;  //
0496 
0497 public:
0498    TMatrixTSub_const() { fRowOff = fColOff = fNrowsSub = fNcolsSub = 0; fMatrix = nullptr; }
0499    TMatrixTSub_const(const TMatrixT   <Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
0500    TMatrixTSub_const(const TMatrixTSym<Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
0501    virtual ~TMatrixTSub_const() { }
0502 
0503    inline const TMatrixTBase<Element> *GetMatrix() const { return fMatrix; }
0504    inline       Int_t                  GetRowOff() const { return fRowOff; }
0505    inline       Int_t                  GetColOff() const { return fColOff; }
0506    inline       Int_t                  GetNrows () const { return fNrowsSub; }
0507    inline       Int_t                  GetNcols () const { return fNcolsSub; }
0508    inline const Element               &operator ()(Int_t rown,Int_t coln) const {
0509       R__ASSERT(fMatrix->IsValid());
0510 
0511       const Element *ptr = fMatrix->GetMatrixArray();
0512       if (rown >= fNrowsSub || rown < 0) {
0513          Error("operator()","Request row(%d) outside matrix range of 0 - %d",rown,fNrowsSub);
0514          return TMatrixTBase<Element>::NaNValue();
0515       }
0516       if (coln >= fNcolsSub || coln < 0) {
0517          Error("operator()","Request column(%d) outside matrix range of 0 - %d",coln,fNcolsSub);
0518          return TMatrixTBase<Element>::NaNValue();
0519       }
0520       const Int_t index = (rown+fRowOff)*fMatrix->GetNcols()+coln+fColOff;
0521       return ptr[index];
0522    }
0523 
0524    ClassDef(TMatrixTSub_const,0)  // Template of Sub Matrix Access class
0525 };
0526 
0527 template<class Element> class TMatrixTSub : public TMatrixTSub_const<Element> {
0528 
0529 public:
0530 
0531    enum {kWorkMax = 100};
0532 
0533    TMatrixTSub() {}
0534    TMatrixTSub(TMatrixT   <Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
0535    TMatrixTSub(TMatrixTSym<Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
0536    TMatrixTSub(const TMatrixTSub<Element> &ms);
0537 
0538    inline Element &operator()(Int_t rown,Int_t coln) {
0539       R__ASSERT(this->fMatrix->IsValid());
0540 
0541       const Element *ptr = this->fMatrix->GetMatrixArray();
0542       if (rown >= this->fNrowsSub || rown < 0) {
0543          Error("operator()","Request row(%d) outside matrix range of 0 - %d",rown,this->fNrowsSub);
0544          return TMatrixTBase<Element>::NaNValue();
0545       }
0546       if (coln >= this->fNcolsSub || coln < 0) {
0547          Error("operator()","Request column(%d) outside matrix range of 0 - %d",coln,this->fNcolsSub);
0548          return TMatrixTBase<Element>::NaNValue();
0549       }
0550       const Int_t index = (rown+this->fRowOff)*this->fMatrix->GetNcols()+coln+this->fColOff;
0551       return (const_cast<Element *>(ptr))[index];
0552    }
0553 
0554    void Rank1Update(const TVectorT<Element> &vec,Element alpha=1.0);
0555 
0556    void operator= (Element val);
0557    void operator+=(Element val);
0558    void operator*=(Element val);
0559 
0560    void operator=(const TMatrixTSub_const<Element> &s);
0561    TMatrixTSub<Element>& operator=(const TMatrixTSub <Element> &s) { operator=((TMatrixTSub_const<Element> &)s); return *this;}
0562    void operator=(const TMatrixTBase     <Element> &m);
0563 
0564    void operator+=(const TMatrixTSub_const<Element> &s);
0565    void operator*=(const TMatrixTSub_const<Element> &s);
0566    void operator+=(const TMatrixTBase     <Element> &m);
0567    void operator*=(const TMatrixT         <Element> &m);
0568    void operator*=(const TMatrixTSym      <Element> &m);
0569 
0570    ClassDefOverride(TMatrixTSub,0)  // Template of Sub Matrix Access class
0571 };
0572 
0573 //////////////////////////////////////////////////////////////////////////
0574 //                                                                      //
0575 // TMatrixTSparseRow_const                                              //
0576 //                                                                      //
0577 // Template class represents a row of TMatrixTSparse                    //
0578 //                                                                      //
0579 //////////////////////////////////////////////////////////////////////////
0580 
0581 template<class Element> class TMatrixTSparseRow_const {
0582 
0583 protected:
0584    const TMatrixTSparse<Element> *fMatrix;  // the matrix I am a row of
0585          Int_t                  fRowInd;  // effective row index
0586          Int_t                  fNindex;  // index range
0587    const Int_t                 *fColPtr;  // column index pointer
0588    const Element               *fDataPtr; // data pointer
0589 
0590 public:
0591    TMatrixTSparseRow_const() { fMatrix = nullptr; fRowInd = 0; fNindex = 0; fColPtr = nullptr; fDataPtr = nullptr; }
0592    TMatrixTSparseRow_const(const TMatrixTSparse<Element> &matrix,Int_t row);
0593    TMatrixTSparseRow_const(const TMatrixTSparseRow_const<Element>& trc):
0594      fMatrix(trc.fMatrix), fRowInd(trc.fRowInd), fNindex(trc.fNindex), fColPtr(trc.fColPtr), fDataPtr(trc.fDataPtr) { }
0595    TMatrixTSparseRow_const<Element>& operator=(const TMatrixTSparseRow_const<Element>& trc) {
0596      if(this != &trc) { fMatrix=trc.fMatrix; fRowInd=trc.fRowInd; fNindex=trc.fNindex; fColPtr=trc.fColPtr; fDataPtr=trc.fDataPtr; }  return *this;}
0597    virtual ~TMatrixTSparseRow_const() { }
0598 
0599    inline const TMatrixTBase<Element> *GetMatrix  () const { return fMatrix; }
0600    inline const Element               *GetDataPtr () const { return fDataPtr; }
0601    inline const Int_t                 *GetColPtr  () const { return fColPtr; }
0602    inline       Int_t                  GetRowIndex() const { return fRowInd; }
0603    inline       Int_t                  GetNindex  () const { return fNindex; }
0604 
0605           Element operator()(Int_t i) const;
0606    inline Element operator[](Int_t i) const { return (*(const TMatrixTSparseRow_const<Element> *)this)(i); }
0607 
0608    ClassDef(TMatrixTSparseRow_const,0)  // Template of Sparse Matrix Row Access class
0609 };
0610 
0611 template<class Element> class TMatrixTSparseRow : public TMatrixTSparseRow_const<Element> {
0612 
0613 public:
0614    TMatrixTSparseRow() {}
0615    TMatrixTSparseRow(TMatrixTSparse<Element> &matrix,Int_t row);
0616    TMatrixTSparseRow(const TMatrixTSparseRow<Element> &mr);
0617 
0618    inline Element *GetDataPtr() const { return const_cast<Element *>(this->fDataPtr); }
0619 
0620           Element  operator()(Int_t i) const;
0621           Element &operator()(Int_t i);
0622    inline Element  operator[](Int_t i) const { return (*(const TMatrixTSparseRow<Element> *)this)(i); }
0623    inline Element &operator[](Int_t i)       { return (*(TMatrixTSparseRow<Element> *)this)(i); }
0624 
0625    void operator= (Element val);
0626    void operator+=(Element val);
0627    void operator*=(Element val);
0628 
0629    void operator=(const TMatrixTSparseRow_const<Element> &r);
0630    TMatrixTSparseRow<Element>& operator=(const TMatrixTSparseRow <Element> &r) { operator=((TMatrixTSparseRow_const<Element> &)r); return *this;}
0631    void operator=(const TVectorT               <Element> &vec);
0632 
0633    void operator+=(const TMatrixTSparseRow_const<Element> &r);
0634    void operator*=(const TMatrixTSparseRow_const<Element> &r);
0635 
0636    ClassDefOverride(TMatrixTSparseRow,0)  // Template of Sparse Matrix Row Access class
0637 };
0638 
0639 //////////////////////////////////////////////////////////////////////////
0640 //                                                                      //
0641 // TMatrixTSparseDiag_const                                             //
0642 //                                                                      //
0643 // Template class represents the diagonal of TMatrixTSparse             //
0644 //                                                                      //
0645 //////////////////////////////////////////////////////////////////////////
0646 
0647 template<class Element> class TMatrixTSparseDiag_const {
0648 
0649 protected:
0650    const TMatrixTSparse<Element> *fMatrix;  //  the matrix I am the diagonal of
0651          Int_t                  fNdiag;   //  number of diag elems, min(nrows,ncols)
0652    const Element               *fDataPtr; //  data pointer
0653 
0654 public:
0655    TMatrixTSparseDiag_const() { fMatrix = nullptr; fNdiag = 0; fDataPtr = nullptr; }
0656    TMatrixTSparseDiag_const(const TMatrixTSparse<Element> &matrix);
0657    TMatrixTSparseDiag_const(const TMatrixTSparseDiag_const<Element>& trc):
0658      fMatrix(trc.fMatrix), fNdiag(trc.fNdiag), fDataPtr(trc.fDataPtr) { }
0659    TMatrixTSparseDiag_const<Element>& operator=(const TMatrixTSparseDiag_const<Element>& trc) {
0660       if(this !=  &trc) { fMatrix=trc.fMatrix; fNdiag=trc.fNdiag; fDataPtr=trc.fDataPtr; } return *this;}
0661    virtual ~TMatrixTSparseDiag_const() { }
0662 
0663    inline const TMatrixTBase<Element> *GetMatrix () const { return fMatrix; }
0664    inline const Element               *GetDataPtr() const { return fDataPtr; }
0665    inline       Int_t                  GetNdiags () const { return fNdiag; }
0666 
0667           Element operator ()(Int_t i) const;
0668    inline Element operator [](Int_t i) const { return (*(const TMatrixTSparseRow_const<Element> *)this)(i); }
0669 
0670    ClassDef(TMatrixTSparseDiag_const,0)  // Template of Sparse Matrix Diagonal Access class
0671 };
0672 
0673 template<class Element> class TMatrixTSparseDiag : public TMatrixTSparseDiag_const<Element> {
0674 
0675 public:
0676    TMatrixTSparseDiag() {}
0677    TMatrixTSparseDiag(TMatrixTSparse<Element> &matrix);
0678    TMatrixTSparseDiag(const TMatrixTSparseDiag<Element> &md);
0679 
0680    inline Element *GetDataPtr() const { return const_cast<Element *>(this->fDataPtr); }
0681 
0682                 Element  operator()(Int_t i) const;
0683                 Element &operator()(Int_t i);
0684    inline       Element  operator[](Int_t i) const { return (*(const TMatrixTSparseDiag<Element> *)this)(i); }
0685    inline       Element &operator[](Int_t i)       { return (*(TMatrixTSparseDiag<Element> *)this)(i); }
0686 
0687    void operator= (Element val);
0688    void operator+=(Element val);
0689    void operator*=(Element val);
0690 
0691    void operator=(const TMatrixTSparseDiag_const<Element> &d);
0692    TMatrixTSparseDiag<Element>& operator=(const TMatrixTSparseDiag <Element> &d) { operator=((TMatrixTSparseDiag_const<Element> &)d); return *this;}
0693    void operator=(const TVectorT                <Element> &vec);
0694 
0695    void operator+=(const TMatrixTSparseDiag_const<Element> &d);
0696    void operator*=(const TMatrixTSparseDiag_const<Element> &d);
0697 
0698    ClassDefOverride(TMatrixTSparseDiag,0)  // Template of Sparse Matrix Diagonal Access class
0699 };
0700 
0701 Double_t Drand(Double_t &ix);
0702 #endif