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

00001 // @(#)root/matrix:$Id: TMatrixTUtils.h 36047 2010-10-04 06:43:15Z brun $
00002 // Authors: Fons Rademakers, Eddy Offermann   Nov 2003
00003 
00004 /*************************************************************************
00005  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers.               *
00006  * All rights reserved.                                                  *
00007  *                                                                       *
00008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
00009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
00010  *************************************************************************/
00011 
00012 #ifndef ROOT_TMatrixTUtils
00013 #define ROOT_TMatrixTUtils
00014 
00015 //////////////////////////////////////////////////////////////////////////
00016 //                                                                      //
00017 // Matrix utility classes.                                              //
00018 //                                                                      //
00019 // Templates of utility classes in the Linear Algebra Package.          //
00020 // The following classes are defined here:                              //
00021 //                                                                      //
00022 // Different matrix views without copying data elements :               //
00023 //   TMatrixTRow_const        TMatrixTRow                               //
00024 //   TMatrixTColumn_const     TMatrixTColumn                            //
00025 //   TMatrixTDiag_const       TMatrixTDiag                              //
00026 //   TMatrixTFlat_const       TMatrixTFlat                              //
00027 //   TMatrixTSub_const        TMatrixTSub                               //
00028 //   TMatrixTSparseRow_const  TMatrixTSparseRow                         //
00029 //   TMatrixTSparseDiag_const TMatrixTSparseDiag                        //
00030 //                                                                      //
00031 //   TElementActionT                                                    //
00032 //   TElementPosActionT                                                 //
00033 //                                                                      //
00034 //////////////////////////////////////////////////////////////////////////
00035 
00036 #ifndef ROOT_TMatrixTBase
00037 #include "TMatrixTBase.h"
00038 #endif
00039 
00040 template<class Element> class TVectorT;
00041 template<class Element> class TMatrixT;
00042 template<class Element> class TMatrixTSym;
00043 template<class Element> class TMatrixTSparse;
00044 
00045 //////////////////////////////////////////////////////////////////////////
00046 //                                                                      //
00047 // TElementActionT                                                      //
00048 //                                                                      //
00049 // A class to do a specific operation on every vector or matrix element //
00050 // (regardless of it position) as the object is being traversed.        //
00051 // This is an abstract class. Derived classes need to implement the     //
00052 // action function Operation().                                         //
00053 //                                                                      //
00054 //////////////////////////////////////////////////////////////////////////
00055 
00056 template<class Element> class TElementActionT {
00057 
00058 #ifndef __CINT__
00059 friend class TMatrixTBase  <Element>;
00060 friend class TMatrixT      <Element>;
00061 friend class TMatrixTSym   <Element>;
00062 friend class TMatrixTSparse<Element>;
00063 friend class TVectorT      <Element>;
00064 #endif
00065 
00066 protected:
00067    virtual ~TElementActionT() { }
00068    virtual void Operation(Element &element) const = 0;
00069 
00070 private:
00071    TElementActionT& operator=(const TElementActionT<Element> &) {return *this;}
00072 };
00073 
00074 //////////////////////////////////////////////////////////////////////////
00075 //                                                                      //
00076 // TElementPosActionT                                                   //
00077 //                                                                      //
00078 // A class to do a specific operation on every vector or matrix element //
00079 // as the object is being traversed. This is an abstract class.         //
00080 // Derived classes need to implement the action function Operation().   //
00081 // In the action function the location of the current element is        //
00082 // known (fI=row, fJ=columns).                                          //
00083 //                                                                      //
00084 //////////////////////////////////////////////////////////////////////////
00085 
00086 template<class Element> class TElementPosActionT {
00087 
00088 #ifndef __CINT__
00089 friend class TMatrixTBase  <Element>;
00090 friend class TMatrixT      <Element>;
00091 friend class TMatrixTSym   <Element>;
00092 friend class TMatrixTSparse<Element>;
00093 friend class TVectorT      <Element>;
00094 #endif
00095 
00096 protected:
00097    mutable Int_t fI; // i position of element being passed to Operation()
00098    mutable Int_t fJ; // j position of element being passed to Operation()
00099    virtual ~TElementPosActionT() { }
00100    virtual void Operation(Element &element) const = 0;
00101 
00102 private:
00103    TElementPosActionT<Element>& operator=(const TElementPosActionT<Element> &) {return *this;}
00104 };
00105 
00106 //////////////////////////////////////////////////////////////////////////
00107 //                                                                      //
00108 // TMatrixTRow_const                                                    //
00109 //                                                                      //
00110 // Template class represents a row of a TMatrixT/TMatrixTSym            //
00111 //                                                                      //
00112 //////////////////////////////////////////////////////////////////////////
00113 
00114 template<class Element> class TMatrixTRow_const {
00115 
00116 protected:
00117    const TMatrixTBase<Element> *fMatrix;  //  the matrix I am a row of
00118          Int_t                  fRowInd;  //  effective row index
00119          Int_t                  fInc;     //  if ptr = @a[row,i], then ptr+inc = @a[row,i+1]
00120    const Element               *fPtr;     //  pointer to the a[row,0]
00121 
00122 public:
00123    TMatrixTRow_const() { fMatrix = 0; fRowInd = 0; fInc = 0; fPtr = 0; }
00124    TMatrixTRow_const(const TMatrixT   <Element> &matrix,Int_t row);
00125    TMatrixTRow_const(const TMatrixTSym<Element> &matrix,Int_t row);
00126   TMatrixTRow_const(const TMatrixTRow_const<Element>& trc):
00127     fMatrix(trc.fMatrix), fRowInd(trc.fRowInd), fInc(trc.fInc), fPtr(trc.fPtr) { }
00128   TMatrixTRow_const<Element>& operator=(const TMatrixTRow_const<Element>& trc) {
00129     fMatrix=trc.fMatrix; fRowInd=trc.fRowInd; fInc=trc.fInc; fPtr=trc.fPtr; return *this;}
00130    virtual ~TMatrixTRow_const() { }
00131 
00132    inline const TMatrixTBase<Element> *GetMatrix  () const { return fMatrix; }
00133    inline       Int_t                  GetRowIndex() const { return fRowInd; }
00134    inline       Int_t                  GetInc     () const { return fInc; }
00135    inline const Element               *GetPtr     () const { return fPtr; }
00136    inline const Element               &operator   ()(Int_t i) const {
00137       R__ASSERT(fMatrix->IsValid());
00138       const Int_t acoln = i-fMatrix->GetColLwb();
00139       if (acoln < fMatrix->GetNcols() && acoln >= 0)
00140          return fPtr[acoln];
00141       else {
00142          Error("operator()","Request col(%d) outside matrix range of %d - %d",
00143                             i,fMatrix->GetColLwb(),fMatrix->GetColLwb()+fMatrix->GetNcols());
00144          return fPtr[0];
00145       }
00146    }
00147    inline const Element               &operator   [](Int_t i) const { return (*(const TMatrixTRow_const<Element> *)this)(i); }
00148 
00149    ClassDef(TMatrixTRow_const,0)  // Template of General Matrix Row Access class
00150 };
00151 
00152 template<class Element> class TMatrixTRow : public TMatrixTRow_const<Element> {
00153 
00154 public:
00155    TMatrixTRow() {}
00156    TMatrixTRow(TMatrixT   <Element> &matrix,Int_t row);
00157    TMatrixTRow(TMatrixTSym<Element> &matrix,Int_t row);
00158    TMatrixTRow(const TMatrixTRow<Element> &mr);
00159 
00160    inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
00161 
00162    inline const Element &operator()(Int_t i) const {
00163       R__ASSERT(this->fMatrix->IsValid());
00164       const Int_t acoln = i-this->fMatrix->GetColLwb();
00165       if (acoln < this->fMatrix->GetNcols() || acoln >= 0)
00166          return (this->fPtr)[acoln];
00167       else {
00168          Error("operator()","Request col(%d) outside matrix range of %d - %d",
00169                             i,this->fMatrix->GetColLwb(),this->fMatrix->GetColLwb()+this->fMatrix->GetNcols());
00170          return (this->fPtr)[0];
00171      }
00172    }
00173    inline       Element &operator()(Int_t i) {
00174       R__ASSERT(this->fMatrix->IsValid());
00175       const Int_t acoln = i-this->fMatrix->GetColLwb();
00176       if (acoln < this->fMatrix->GetNcols() && acoln >= 0)
00177          return (const_cast<Element *>(this->fPtr))[acoln];
00178       else {
00179          Error("operator()","Request col(%d) outside matrix range of %d - %d",
00180                             i,this->fMatrix->GetColLwb(),this->fMatrix->GetColLwb()+this->fMatrix->GetNcols());
00181          return (const_cast<Element *>(this->fPtr))[0];
00182       }
00183    }
00184    inline const Element &operator[](Int_t i) const { return (*(const TMatrixTRow<Element> *)this)(i); }
00185    inline       Element &operator[](Int_t i)       { return (*(      TMatrixTRow<Element> *)this)(i); }
00186 
00187    void operator= (Element val);
00188    void operator+=(Element val);
00189    void operator*=(Element val);
00190 
00191    void operator=(const TMatrixTRow_const<Element> &r);
00192    TMatrixTRow<Element>& operator=(const TMatrixTRow      <Element> &r) { operator=((TMatrixTRow_const<Element> &)r); return *this;}
00193    void operator=(const TVectorT         <Element> &vec);
00194 
00195    void operator+=(const TMatrixTRow_const<Element> &r);
00196    void operator*=(const TMatrixTRow_const<Element> &r);
00197 
00198    ClassDef(TMatrixTRow,0)  // Template of General Matrix Row Access class
00199 };
00200 
00201 //////////////////////////////////////////////////////////////////////////
00202 //                                                                      //
00203 // TMatrixTColumn_const                                                 //
00204 //                                                                      //
00205 // Template class represents a column of a TMatrixT/TMatrixTSym         //
00206 //                                                                      //
00207 //////////////////////////////////////////////////////////////////////////
00208 
00209 template<class Element> class TMatrixTColumn_const {
00210 
00211 protected:
00212    const TMatrixTBase<Element> *fMatrix;  //  the matrix I am a column of
00213          Int_t                  fColInd;  //  effective column index
00214          Int_t                  fInc;     //  if ptr = @a[i,col], then ptr+inc = @a[i+1,col]
00215    const Element               *fPtr;     //  pointer to the a[0,col] column
00216 
00217 public:
00218    TMatrixTColumn_const() { fMatrix = 0; fColInd = 0; fInc = 0; fPtr = 0; }
00219    TMatrixTColumn_const(const TMatrixT   <Element> &matrix,Int_t col);
00220    TMatrixTColumn_const(const TMatrixTSym<Element> &matrix,Int_t col);
00221    TMatrixTColumn_const(const TMatrixTColumn_const<Element>& trc):
00222      fMatrix(trc.fMatrix), fColInd(trc.fColInd), fInc(trc.fInc), fPtr(trc.fPtr) { }
00223    TMatrixTColumn_const<Element>& operator=(const TMatrixTColumn_const<Element>& trc) {
00224      fMatrix=trc.fMatrix; fColInd=trc.fColInd; fInc=trc.fInc; fPtr=trc.fPtr; return *this;}
00225    virtual ~TMatrixTColumn_const() { }
00226 
00227    inline const TMatrixTBase <Element> *GetMatrix  () const { return fMatrix; }
00228    inline       Int_t                   GetColIndex() const { return fColInd; }
00229    inline       Int_t                   GetInc     () const { return fInc; }
00230    inline const Element                *GetPtr     () const { return fPtr; }
00231    inline const Element                &operator   ()(Int_t i) const {
00232       R__ASSERT(fMatrix->IsValid());
00233       const Int_t arown = i-fMatrix->GetRowLwb();
00234       if (arown < fMatrix->GetNrows() && arown >= 0)
00235          return fPtr[arown*fInc];
00236       else {
00237          Error("operator()","Request row(%d) outside matrix range of %d - %d",
00238                             i,fMatrix->GetRowLwb(),fMatrix->GetRowLwb()+fMatrix->GetNrows());
00239          return fPtr[0];
00240       }
00241    }
00242    inline const Element                &operator [](Int_t i) const { return (*(const TMatrixTColumn_const<Element> *)this)(i); }
00243 
00244    ClassDef(TMatrixTColumn_const,0)  // Template of General Matrix Column Access class
00245 };
00246 
00247 template<class Element> class TMatrixTColumn : public TMatrixTColumn_const<Element> {
00248 
00249 public:
00250    TMatrixTColumn() {}
00251    TMatrixTColumn(TMatrixT   <Element>&matrix,Int_t col);
00252    TMatrixTColumn(TMatrixTSym<Element>&matrix,Int_t col);
00253    TMatrixTColumn(const TMatrixTColumn <Element>&mc);
00254 
00255    inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
00256  
00257    inline const Element &operator()(Int_t i) const {
00258       R__ASSERT(this->fMatrix->IsValid());
00259       const Int_t arown = i-this->fMatrix->GetRowLwb();
00260       if (arown < this->fMatrix->GetNrows() && arown >= 0)
00261          return (this->fPtr)[arown*this->fInc];
00262       else {
00263          Error("operator()","Request row(%d) outside matrix range of %d - %d",
00264                             i,this->fMatrix->GetRowLwb(),this->fMatrix->GetRowLwb()+this->fMatrix->GetNrows());
00265          return (this->fPtr)[0];
00266       }
00267    }
00268    inline       Element &operator()(Int_t i) {
00269       R__ASSERT(this->fMatrix->IsValid());
00270       const Int_t arown = i-this->fMatrix->GetRowLwb();
00271 
00272       if (arown < this->fMatrix->GetNrows() && arown >= 0)
00273          return (const_cast<Element *>(this->fPtr))[arown*this->fInc];
00274       else {
00275          Error("operator()","Request row(%d) outside matrix range of %d - %d",
00276                             i,this->fMatrix->GetRowLwb(),this->fMatrix->GetRowLwb()+this->fMatrix->GetNrows());
00277          return (const_cast<Element *>(this->fPtr))[0];
00278       }
00279    }
00280    inline const Element &operator[](Int_t i) const { return (*(const TMatrixTColumn<Element> *)this)(i); }
00281    inline       Element &operator[](Int_t i)       { return (*(      TMatrixTColumn<Element> *)this)(i); }
00282 
00283    void operator= (Element val);
00284    void operator+=(Element val);
00285    void operator*=(Element val);
00286 
00287    void operator=(const TMatrixTColumn_const<Element> &c);
00288    TMatrixTColumn<Element>& operator=(const TMatrixTColumn <Element> &c) { operator=((TMatrixTColumn_const<Element> &)c); return *this;}
00289    void operator=(const TVectorT            <Element> &vec);
00290 
00291    void operator+=(const TMatrixTColumn_const<Element> &c);
00292    void operator*=(const TMatrixTColumn_const<Element> &c);
00293 
00294    ClassDef(TMatrixTColumn,0)  // Template of General Matrix Column Access class
00295 };
00296 
00297 //////////////////////////////////////////////////////////////////////////
00298 //                                                                      //
00299 // TMatrixTDiag_const                                                   //
00300 //                                                                      //
00301 // Template class represents the diagonal of a TMatrixT/TMatrixTSym     //
00302 //                                                                      //
00303 //////////////////////////////////////////////////////////////////////////
00304 
00305 template<class Element> class TMatrixTDiag_const {
00306 
00307 protected:
00308    const TMatrixTBase<Element> *fMatrix;  //  the matrix I am the diagonal of
00309          Int_t                  fInc;     //  if ptr=@a[i,i], then ptr+inc = @a[i+1,i+1]
00310          Int_t                  fNdiag;   //  number of diag elems, min(nrows,ncols)
00311    const Element               *fPtr;     //  pointer to the a[0,0]
00312 
00313 public:
00314    TMatrixTDiag_const() { fMatrix = 0; fInc = 0; fNdiag = 0; fPtr = 0; }
00315    TMatrixTDiag_const(const TMatrixT   <Element> &matrix);
00316    TMatrixTDiag_const(const TMatrixTSym<Element> &matrix);
00317    TMatrixTDiag_const(const TMatrixTDiag_const<Element>& trc):
00318     fMatrix(trc.fMatrix), fInc(trc.fInc), fNdiag(trc.fNdiag), fPtr(trc.fPtr) { }
00319    TMatrixTDiag_const<Element>& operator=(const TMatrixTDiag_const<Element>& trc) {
00320      fMatrix=trc.fMatrix; fInc=trc.fInc; fNdiag=trc.fNdiag; fPtr=trc.fPtr; return *this;}
00321    virtual ~TMatrixTDiag_const() { }
00322 
00323    inline const TMatrixTBase<Element> *GetMatrix() const { return fMatrix; }
00324    inline const Element               *GetPtr   () const { return fPtr; }
00325    inline       Int_t                  GetInc   () const { return fInc; }
00326    inline const Element               &operator ()(Int_t i) const {
00327       R__ASSERT(fMatrix->IsValid());
00328       if (i < fNdiag && i >= 0)
00329          return fPtr[i*fInc];
00330       else {
00331          Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,fNdiag);
00332          return fPtr[0];
00333       }
00334    }
00335    inline const Element               &operator [](Int_t i) const { return (*(const TMatrixTDiag_const<Element> *)this)(i); }
00336 
00337    Int_t GetNdiags() const { return fNdiag; }
00338 
00339    ClassDef(TMatrixTDiag_const,0)  // Template of General Matrix Diagonal Access class
00340 };
00341 
00342 template<class Element> class TMatrixTDiag : public TMatrixTDiag_const<Element> {
00343 
00344 public:
00345    TMatrixTDiag() {}
00346    TMatrixTDiag(TMatrixT   <Element>&matrix);
00347    TMatrixTDiag(TMatrixTSym<Element>&matrix);
00348    TMatrixTDiag(const TMatrixTDiag<Element> &md);
00349 
00350    inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
00351 
00352    inline const Element &operator()(Int_t i) const {
00353       R__ASSERT(this->fMatrix->IsValid());
00354       if (i < this->fNdiag && i >= 0)
00355          return (this->fPtr)[i*this->fInc];
00356       else {
00357          Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,this->fNdiag);
00358          return (this->fPtr)[0];
00359       }
00360    }
00361    inline       Element &operator()(Int_t i) {
00362       R__ASSERT(this->fMatrix->IsValid());
00363       if (i < this->fNdiag && i >= 0)
00364          return (const_cast<Element *>(this->fPtr))[i*this->fInc];
00365       else {
00366          Error("operator()","Request diagonal(%d) outside matrix range of 0 - %d",i,this->fNdiag);
00367          return (const_cast<Element *>(this->fPtr))[0];
00368       }
00369    }
00370    inline const Element &operator[](Int_t i) const { return (*(const TMatrixTDiag<Element> *)this)(i); }
00371    inline       Element &operator[](Int_t i)       { return (*(      TMatrixTDiag *)this)(i); }
00372 
00373    void operator= (Element val);
00374    void operator+=(Element val);
00375    void operator*=(Element val);
00376 
00377    void operator=(const TMatrixTDiag_const<Element> &d);
00378    TMatrixTDiag<Element>& operator=(const TMatrixTDiag <Element> &d) { operator=((TMatrixTDiag_const<Element> &)d); return *this;}
00379    void operator=(const TVectorT          <Element> &vec);
00380 
00381    void operator+=(const TMatrixTDiag_const<Element> &d);
00382    void operator*=(const TMatrixTDiag_const<Element> &d);
00383 
00384    ClassDef(TMatrixTDiag,0)  // Template of General Matrix Diagonal Access class
00385 };
00386 
00387 //////////////////////////////////////////////////////////////////////////
00388 //                                                                      //
00389 // TMatrixTFlat_const                                                   //
00390 //                                                                      //
00391 // Template class represents a flat TMatrixT/TMatrixTSym                //
00392 //                                                                      //
00393 //////////////////////////////////////////////////////////////////////////
00394 
00395 template<class Element> class TMatrixTFlat_const {
00396 
00397 protected:
00398    const TMatrixTBase<Element> *fMatrix;  //  the matrix I am the diagonal of
00399          Int_t                  fNelems;  //
00400    const Element               *fPtr;     //  pointer to the a[0,0]
00401 
00402 public:
00403    TMatrixTFlat_const() { fMatrix = 0; fNelems = 0; fPtr = 0; }
00404    TMatrixTFlat_const(const TMatrixT   <Element> &matrix);
00405    TMatrixTFlat_const(const TMatrixTSym<Element> &matrix);
00406    TMatrixTFlat_const(const TMatrixTFlat_const<Element>& trc):
00407      fMatrix(trc.fMatrix), fNelems(trc.fNelems), fPtr(trc.fPtr) { }
00408    TMatrixTFlat_const<Element>& operator=(const TMatrixTFlat_const<Element>& trc) {
00409     fMatrix=trc.fMatrix; fNelems=trc.fNelems; fPtr=trc.fPtr; return *this;}
00410    virtual ~TMatrixTFlat_const() { }
00411 
00412    inline const TMatrixTBase<Element> *GetMatrix() const { return fMatrix; }
00413    inline const Element               *GetPtr   () const { return fPtr; }
00414    inline const Element               &operator ()(Int_t i) const {
00415       R__ASSERT(fMatrix->IsValid());
00416       if (i < fNelems && i >= 0)
00417          return fPtr[i];
00418       else {
00419          Error("operator()","Request element(%d) outside matrix range of 0 - %d",i,fNelems);
00420          return fPtr[0];
00421       }
00422    }
00423    inline const Element               &operator [](Int_t i) const { return (*(const TMatrixTFlat_const<Element> *)this)(i); }
00424 
00425    ClassDef(TMatrixTFlat_const,0)  // Template of General Matrix Flat Representation class
00426 };
00427 
00428 template<class Element> class TMatrixTFlat : public TMatrixTFlat_const<Element> {
00429 
00430 public:
00431    TMatrixTFlat() {}
00432    TMatrixTFlat(TMatrixT   <Element> &matrix);
00433    TMatrixTFlat(TMatrixTSym<Element> &matrix);
00434    TMatrixTFlat(const TMatrixTFlat<Element> &mf);
00435 
00436    inline Element *GetPtr() const { return const_cast<Element *>(this->fPtr); }
00437 
00438    inline const Element &operator()(Int_t i) const {
00439       R__ASSERT(this->fMatrix->IsValid());
00440       if (i < this->fNelems && i >= 0)
00441          return (this->fPtr)[i];
00442       else {
00443          Error("operator()","Request element(%d) outside matrix range of 0 - %d",i,this->fNelems);
00444          return (this->fPtr)[0];
00445       }
00446    }
00447    inline       Element &operator()(Int_t i) {
00448       R__ASSERT(this->fMatrix->IsValid());
00449       if (i < this->fNelems && i >= 0)
00450          return (const_cast<Element *>(this->fPtr))[i];
00451       else {
00452          Error("operator()","Request element(%d) outside matrix range of 0 - %d",i,this->fNelems);
00453          return (const_cast<Element *>(this->fPtr))[0];
00454       }
00455    }
00456    inline const Element &operator[](Int_t i) const { return (*(const TMatrixTFlat<Element> *)this)(i); }
00457    inline       Element &operator[](Int_t i)       { return (*(      TMatrixTFlat<Element> *)this)(i); }
00458 
00459    void operator= (Element val);
00460    void operator+=(Element val);
00461    void operator*=(Element val);
00462 
00463    void operator=(const TMatrixTFlat_const<Element> &f);
00464    TMatrixTFlat<Element>& operator=(const TMatrixTFlat <Element> &f) { operator=((TMatrixTFlat_const<Element> &)f); return *this;}
00465    void operator=(const TVectorT          <Element> &vec);
00466 
00467    void operator+=(const TMatrixTFlat_const<Element> &f);
00468    void operator*=(const TMatrixTFlat_const<Element> &f);
00469 
00470    ClassDef(TMatrixTFlat,0)  // Template of General Matrix Flat Representation class
00471 };
00472 
00473 //////////////////////////////////////////////////////////////////////////
00474 //                                                                      //
00475 // TMatrixTSub_const                                                    //
00476 //                                                                      //
00477 // Template class represents a sub matrix of TMatrixT/TMatrixTSym       //
00478 //                                                                      //
00479 //////////////////////////////////////////////////////////////////////////
00480 
00481 template<class Element> class TMatrixTSub_const {
00482 
00483 protected:
00484    const TMatrixTBase<Element> *fMatrix;    //  the matrix I am a submatrix of
00485          Int_t                  fRowOff;    //
00486          Int_t                  fColOff;    //
00487          Int_t                  fNrowsSub;  //
00488          Int_t                  fNcolsSub;  //
00489 
00490 public:
00491    TMatrixTSub_const() { fRowOff = fColOff = fNrowsSub = fNcolsSub = 0; fMatrix = 0; }
00492    TMatrixTSub_const(const TMatrixT   <Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
00493    TMatrixTSub_const(const TMatrixTSym<Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
00494    virtual ~TMatrixTSub_const() { }
00495 
00496    inline const TMatrixTBase<Element> *GetMatrix() const { return fMatrix; }
00497    inline       Int_t                  GetRowOff() const { return fRowOff; }
00498    inline       Int_t                  GetColOff() const { return fColOff; }
00499    inline       Int_t                  GetNrows () const { return fNrowsSub; }
00500    inline       Int_t                  GetNcols () const { return fNcolsSub; }
00501    inline const Element               &operator ()(Int_t rown,Int_t coln) const {
00502       R__ASSERT(fMatrix->IsValid());
00503 
00504       const Element *ptr = fMatrix->GetMatrixArray();
00505       if (rown >= fNrowsSub || rown < 0) {
00506          Error("operator()","Request row(%d) outside matrix range of 0 - %d",rown,fNrowsSub);
00507          return ptr[0];
00508       }
00509       if (coln >= fNcolsSub || coln < 0) {
00510          Error("operator()","Request column(%d) outside matrix range of 0 - %d",coln,fNcolsSub);
00511          return ptr[0];
00512       }
00513       const Int_t index = (rown+fRowOff)*fMatrix->GetNcols()+coln+fColOff;
00514       return ptr[index];
00515    }
00516 
00517    ClassDef(TMatrixTSub_const,0)  // Template of Sub Matrix Access class
00518 };
00519 
00520 template<class Element> class TMatrixTSub : public TMatrixTSub_const<Element> {
00521 
00522 public:
00523 
00524    enum {kWorkMax = 100};
00525 
00526    TMatrixTSub() {}
00527    TMatrixTSub(TMatrixT   <Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
00528    TMatrixTSub(TMatrixTSym<Element> &matrix,Int_t row_lwb,Int_t row_upb,Int_t col_lwb,Int_t col_upb);
00529    TMatrixTSub(const TMatrixTSub<Element> &ms);
00530 
00531    inline Element &operator()(Int_t rown,Int_t coln) {
00532       R__ASSERT(this->fMatrix->IsValid());
00533 
00534       const Element *ptr = this->fMatrix->GetMatrixArray();
00535       if (rown >= this->fNrowsSub || rown < 0) {
00536          Error("operator()","Request row(%d) outside matrix range of 0 - %d",rown,this->fNrowsSub);
00537          return (const_cast<Element *>(ptr))[0];
00538       }
00539       if (coln >= this->fNcolsSub || coln < 0) {
00540          Error("operator()","Request column(%d) outside matrix range of 0 - %d",coln,this->fNcolsSub);
00541          return (const_cast<Element *>(ptr))[0];
00542       }
00543       const Int_t index = (rown+this->fRowOff)*this->fMatrix->GetNcols()+coln+this->fColOff;
00544       return (const_cast<Element *>(ptr))[index];
00545    }
00546 
00547    void Rank1Update(const TVectorT<Element> &vec,Element alpha=1.0);
00548 
00549    void operator= (Element val);
00550    void operator+=(Element val);
00551    void operator*=(Element val);
00552 
00553    void operator=(const TMatrixTSub_const<Element> &s);
00554    TMatrixTSub<Element>& operator=(const TMatrixTSub <Element> &s) { operator=((TMatrixTSub_const<Element> &)s); return *this;}
00555    void operator=(const TMatrixTBase     <Element> &m);
00556 
00557    void operator+=(const TMatrixTSub_const<Element> &s);
00558    void operator*=(const TMatrixTSub_const<Element> &s);
00559    void operator+=(const TMatrixTBase     <Element> &m);
00560    void operator*=(const TMatrixT         <Element> &m);
00561    void operator*=(const TMatrixTSym      <Element> &m);
00562 
00563    ClassDef(TMatrixTSub,0)  // Template of Sub Matrix Access class
00564 };
00565 
00566 //////////////////////////////////////////////////////////////////////////
00567 //                                                                      //
00568 // TMatrixTSparseRow_const                                              //
00569 //                                                                      //
00570 // Template class represents a row of TMatrixTSparse                    //
00571 //                                                                      //
00572 //////////////////////////////////////////////////////////////////////////
00573 
00574 template<class Element> class TMatrixTSparseRow_const {
00575 
00576 protected:
00577    const TMatrixTBase<Element> *fMatrix;  // the matrix I am a row of
00578          Int_t                  fRowInd;  // effective row index
00579          Int_t                  fNindex;  // index range
00580    const Int_t                 *fColPtr;  // column index pointer
00581    const Element               *fDataPtr; // data pointer
00582 
00583 public:
00584    TMatrixTSparseRow_const() { fMatrix = 0; fRowInd = 0; fNindex = 0; fColPtr = 0; fDataPtr = 0; }
00585    TMatrixTSparseRow_const(const TMatrixTSparse<Element> &matrix,Int_t row);
00586    TMatrixTSparseRow_const(const TMatrixTSparseRow_const<Element>& trc):
00587      fMatrix(trc.fMatrix), fRowInd(trc.fRowInd), fNindex(trc.fNindex), fColPtr(trc.fColPtr), fDataPtr(trc.fDataPtr) { }
00588    TMatrixTSparseRow_const<Element>& operator=(const TMatrixTSparseRow_const<Element>& trc) {
00589      fMatrix=trc.fMatrix; fRowInd=trc.fRowInd; fNindex=trc.fNindex; fColPtr=trc.fColPtr; fDataPtr=trc.fDataPtr; return *this;}
00590    virtual ~TMatrixTSparseRow_const() { }
00591 
00592    inline const TMatrixTBase<Element> *GetMatrix  () const { return fMatrix; }
00593    inline const Element               *GetDataPtr () const { return fDataPtr; }
00594    inline const Int_t                 *GetColPtr  () const { return fColPtr; }
00595    inline       Int_t                  GetRowIndex() const { return fRowInd; }
00596    inline       Int_t                  GetNindex  () const { return fNindex; }
00597 
00598           Element operator()(Int_t i) const;
00599    inline Element operator[](Int_t i) const { return (*(const TMatrixTSparseRow_const<Element> *)this)(i); }
00600 
00601    ClassDef(TMatrixTSparseRow_const,0)  // Template of Sparse Matrix Row Access class
00602 };
00603 
00604 template<class Element> class TMatrixTSparseRow : public TMatrixTSparseRow_const<Element> {
00605 
00606 public:
00607    TMatrixTSparseRow() {}
00608    TMatrixTSparseRow(TMatrixTSparse<Element> &matrix,Int_t row);
00609    TMatrixTSparseRow(const TMatrixTSparseRow<Element> &mr);
00610 
00611    inline Element *GetDataPtr() const { return const_cast<Element *>(this->fDataPtr); }
00612 
00613           Element  operator()(Int_t i) const;
00614           Element &operator()(Int_t i);
00615    inline Element  operator[](Int_t i) const { return (*(const TMatrixTSparseRow<Element> *)this)(i); }
00616    inline Element &operator[](Int_t i)       { return (*(TMatrixTSparseRow<Element> *)this)(i); }
00617 
00618    void operator= (Element val);
00619    void operator+=(Element val);
00620    void operator*=(Element val);
00621 
00622    void operator=(const TMatrixTSparseRow_const<Element> &r);
00623    TMatrixTSparseRow<Element>& operator=(const TMatrixTSparseRow <Element> &r) { operator=((TMatrixTSparseRow_const<Element> &)r); return *this;}
00624    void operator=(const TVectorT               <Element> &vec);
00625 
00626    void operator+=(const TMatrixTSparseRow_const<Element> &r);
00627    void operator*=(const TMatrixTSparseRow_const<Element> &r);
00628 
00629    ClassDef(TMatrixTSparseRow,0)  // Template of Sparse Matrix Row Access class
00630 };
00631 
00632 //////////////////////////////////////////////////////////////////////////
00633 //                                                                      //
00634 // TMatrixTSparseDiag_const                                             //
00635 //                                                                      //
00636 // Template class represents the diagonal of TMatrixTSparse             //
00637 //                                                                      //
00638 //////////////////////////////////////////////////////////////////////////
00639 
00640 template<class Element> class TMatrixTSparseDiag_const {
00641 
00642 protected:
00643    const TMatrixTBase<Element> *fMatrix;  //  the matrix I am the diagonal of
00644          Int_t                  fNdiag;   //  number of diag elems, min(nrows,ncols)
00645    const Element               *fDataPtr; //  data pointer
00646 
00647 public:
00648    TMatrixTSparseDiag_const() { fMatrix = 0; fNdiag = 0; fDataPtr = 0; }
00649    TMatrixTSparseDiag_const(const TMatrixTSparse<Element> &matrix);
00650    TMatrixTSparseDiag_const(const TMatrixTSparseDiag_const<Element>& trc):
00651      fMatrix(trc.fMatrix), fNdiag(trc.fNdiag), fDataPtr(trc.fDataPtr) { }
00652    TMatrixTSparseDiag_const<Element>& operator=(const TMatrixTSparseDiag_const<Element>& trc) {
00653      fMatrix=trc.fMatrix; fNdiag=trc.fNdiag; fDataPtr=trc.fDataPtr; return *this;}
00654    virtual ~TMatrixTSparseDiag_const() { }
00655 
00656    inline const TMatrixTBase<Element> *GetMatrix () const { return fMatrix; }
00657    inline const Element               *GetDataPtr() const { return fDataPtr; }
00658    inline       Int_t                  GetNdiags () const { return fNdiag; }
00659 
00660           Element operator ()(Int_t i) const;
00661    inline Element operator [](Int_t i) const { return (*(const TMatrixTSparseRow_const<Element> *)this)(i); }
00662 
00663    ClassDef(TMatrixTSparseDiag_const,0)  // Template of Sparse Matrix Diagonal Access class
00664 };
00665 
00666 template<class Element> class TMatrixTSparseDiag : public TMatrixTSparseDiag_const<Element> {
00667 
00668 public:
00669    TMatrixTSparseDiag() {}
00670    TMatrixTSparseDiag(TMatrixTSparse<Element> &matrix);
00671    TMatrixTSparseDiag(const TMatrixTSparseDiag<Element> &md);
00672 
00673    inline Element *GetDataPtr() const { return const_cast<Element *>(this->fDataPtr); }
00674 
00675                 Element  operator()(Int_t i) const;
00676                 Element &operator()(Int_t i);
00677    inline       Element  operator[](Int_t i) const { return (*(const TMatrixTSparseDiag<Element> *)this)(i); }
00678    inline       Element &operator[](Int_t i)       { return (*(TMatrixTSparseDiag<Element> *)this)(i); }
00679 
00680    void operator= (Element val);
00681    void operator+=(Element val);
00682    void operator*=(Element val);
00683 
00684    void operator=(const TMatrixTSparseDiag_const<Element> &d);
00685    TMatrixTSparseDiag<Element>& operator=(const TMatrixTSparseDiag <Element> &d) { operator=((TMatrixTSparseDiag_const<Element> &)d); return *this;}
00686    void operator=(const TVectorT                <Element> &vec);
00687 
00688    void operator+=(const TMatrixTSparseDiag_const<Element> &d);
00689    void operator*=(const TMatrixTSparseDiag_const<Element> &d);
00690 
00691    ClassDef(TMatrixTSparseDiag,0)  // Template of Sparse Matrix Diagonal Access class
00692 };
00693 
00694 Double_t Drand(Double_t &ix);
00695 #endif