STEER/AliSymMatrix.h

   1 #ifndef ALISYMMATRIX_H
   2 #define ALISYMMATRIX_H
   3
   4 #include <string.h>
   5 #include <TObject.h>
   6 #include <TVectorD.h>
   7 #include "AliMatrixSq.h"
   8
   9
  10
  11 class AliSymMatrix : public AliMatrixSq {
  12   //
  13  public:
  14   AliSymMatrix();
  15   AliSymMatrix(Int_t size);
  16   AliSymMatrix(const AliSymMatrix &mat);
  17   virtual ~AliSymMatrix();
  18   //
  19   void          Clear(Option_t* option="");
  20   void          Reset();
  21   //
  22   Int_t         GetSize()                                        const {return fNrowIndex;}
  23   Float_t       GetDensity()                                     const;
  24   AliSymMatrix& operator=(const AliSymMatrix& src);
  25   Double_t      operator()(Int_t rown, Int_t coln)               const;
  26   Double_t&     operator()(Int_t rown, Int_t coln);
  27   //
  28   Double_t      DiagElem(Int_t r)                                const {return (*(const AliSymMatrix*)this)(r,r);}
  29   Double_t&     DiagElem(Int_t r)                                      {return (*this)(r,r);}
  30   //
  31   void          Print(Option_t* option="")                       const;
  32   void          AddRows(int nrows=1);
  33   //
  34   void          Scale(Double_t coeff);
  35   void          MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
  36   void          MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const;
  37   //
  38   // ---------------------------------- Dummy methods of MatrixBase
  39   virtual       const Double_t   *GetMatrixArray  () const {return fElems;};
  40   virtual             Double_t   *GetMatrixArray  ()       {return (Double_t*)fElems;}
  41   virtual       const Int_t      *GetRowIndexArray() const {Error("GetRowIndexArray","Dummy"); return 0;};
  42   virtual             Int_t      *GetRowIndexArray()       {Error("GetRowIndexArray","Dummy"); return 0;};
  43   virtual       const Int_t      *GetColIndexArray() const {Error("GetColIndexArray","Dummy"); return 0;};
  44   virtual             Int_t      *GetColIndexArray()       {Error("GetColIndexArray","Dummy"); return 0;};
  45   virtual             TMatrixDBase &SetRowIndexArray(Int_t *) {Error("SetRowIndexArray","Dummy"); return *this;}
  46   virtual             TMatrixDBase &SetColIndexArray(Int_t *) {Error("SetColIndexArray","Dummy"); return *this;}
  47   virtual             TMatrixDBase &GetSub(Int_t,Int_t,Int_t,Int_t,TMatrixDBase &,Option_t *) const {Error("GetSub","Dummy"); return *((TMatrixDBase*)this);}
  48   virtual             TMatrixDBase &SetSub(Int_t,Int_t,const TMatrixDBase &) {Error("GetSub","Dummy"); return *this;}
  49   virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
  50   virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
  51   //
  52   // ----------------------------- Choleski methods ----------------------------------------
  53   AliSymMatrix*       DecomposeChol();                  //Obtain Cholesky decomposition L matrix
  54   void                InvertChol(AliSymMatrix* mchol);  //Invert using provided Choleski decomposition
  55   Bool_t              InvertChol();                     //Invert
  56   Bool_t              SolveChol(Double_t *brhs, Bool_t invert=kFALSE);
  57   Bool_t              SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert=kFALSE);
  58   Bool_t              SolveChol(TVectorD &brhs, Bool_t invert=kFALSE);
  59   Bool_t              SolveChol(TVectorD &brhs, TVectorD& bsol,Bool_t invert=kFALSE);
  60   //
  61   int                 SolveSpmInv(double *vecB, Bool_t stabilize=kTRUE);
  62
  63  protected:
  64   virtual             Int_t GetIndex(Int_t row,Int_t col)      const;
  65   Double_t            GetEl(Int_t row, Int_t col)              const {return operator()(row,col);}
  66   void                SetEl(Int_t row, Int_t col,Double_t val)       {operator()(row,col) = val;}
  67   //
  68  protected:
  69   Double_t*  fElems;     //   Elements booked by constructor
  70   Double_t** fElemsAdd;  //   Elements (rows) added dynamicaly
  71   //
  72   static AliSymMatrix* fgBuffer;  // buffer for fast solution
  73   static Int_t         fgCopyCnt;  // matrix copy counter
  74   ClassDef(AliSymMatrix,0) //Symmetric Matrix Class
  75 };
  76
  77
  78 //___________________________________________________________
  79 inline Int_t AliSymMatrix::GetIndex(Int_t row,Int_t col) const
  80 {
  81   // lower triangle is actually filled
  82   return ((row*(row+1))>>1) + col;
  83 }
  84
  85 //___________________________________________________________
  86 inline Double_t AliSymMatrix::operator()(Int_t row, Int_t col) const
  87 {
  88   //
  89   if (row<col) Swap(row,col);
  90   if (row>=fNrowIndex) return 0;
  91   return (const Double_t&) (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
  92 }
  93
  94 //___________________________________________________________
  95 inline Double_t& AliSymMatrix::operator()(Int_t row, Int_t col)
  96 {
  97   if (row<col) Swap(row,col);
  98   if (row>=fNrowIndex) AddRows(row-fNrowIndex+1);
  99   return (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
 100 }
 101
 102 //___________________________________________________________
 103 inline void AliSymMatrix::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const
 104 {
 105   MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
 106 }
 107
 108 //___________________________________________________________
 109 inline void AliSymMatrix::Scale(Double_t coeff)
 110 {
 111   for (int i=fNrowIndex;i--;) for (int j=i;j--;) { double& el = operator()(i,j); if (el) el *= coeff;}
 112 }
 113
 114
 115 #endif