STEER/AliSymMatrix.h

   1 #ifndef ALISYMMATRIX_H
   2 #define ALISYMMATRIX_H
   3 /**********************************************************************************************/
   4 /* Fast symmetric matrix with dynamically expandable size.                                    */
   5 /* Only part can be used for matrix operations. It is defined as:                             */
   6 /* fNCols: rows built by constructor (GetSizeBooked)                                          */
   7 /* fNRows: number of rows added dynamically (automatically added on assignment to row)        */
   8 /*         GetNRowAdded                                                                       */
   9 /* fNRowIndex: total size (fNCols+fNRows), GetSize                                            */
  10 /* fRowLwb   : actual size to used for given operation, by default = total size, GetSizeUsed  */
  11 /*                                                                                            */
  12 /* Author: ruben.shahoyan@cern.ch                                                             */
  13 /*                                                                                            */
  14 /**********************************************************************************************/
  15
  16 #include <string.h>
  17 #include <TObject.h>
  18 #include <TVectorD.h>
  19 #include "AliMatrixSq.h"
  20
  21
  22
  23 class AliSymMatrix : public AliMatrixSq {
  24   //
  25  public:
  26   AliSymMatrix();
  27   AliSymMatrix(Int_t size);
  28   AliSymMatrix(const AliSymMatrix &mat);
  29   virtual ~AliSymMatrix();
  30   //
  31   void          Clear(Option_t* option="");
  32   void          Reset();
  33   //
  34   Int_t         GetSize()                                        const {return fNrowIndex;}
  35   Int_t         GetSizeUsed()                                    const {return fRowLwb;}
  36   Int_t         GetSizeBooked()                                  const {return fNcols;}
  37   Int_t         GetSizeAdded()                                   const {return fNrows;}
  38   Float_t       GetDensity()                                     const;
  39   AliSymMatrix& operator=(const AliSymMatrix& src);
  40   AliSymMatrix& operator+=(const AliSymMatrix& src);
  41   Double_t      operator()(Int_t rown, Int_t coln)               const;
  42   Double_t&     operator()(Int_t rown, Int_t coln);
  43   //
  44   Double_t      DiagElem(Int_t r)                                const {return (*(const AliSymMatrix*)this)(r,r);}
  45   Double_t&     DiagElem(Int_t r)                                      {return (*this)(r,r);}
  46   //
  47   Double_t*     GetRow(Int_t r);
  48   //
  49   void          Print(Option_t* option="")                       const;
  50   void          AddRows(int nrows=1);
  51   void          SetSizeUsed(Int_t sz)                                  {fRowLwb = sz;}
  52   //
  53   void          Scale(Double_t coeff);
  54   void          MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
  55   void          MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const;
  56   void          AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
  57   //
  58   // ---------------------------------- Dummy methods of MatrixBase
  59   virtual       const Double_t   *GetMatrixArray  () const {return fElems;};
  60   virtual             Double_t   *GetMatrixArray  ()       {return (Double_t*)fElems;}
  61   virtual       const Int_t      *GetRowIndexArray() const {Error("GetRowIndexArray","Dummy"); return 0;};
  62   virtual             Int_t      *GetRowIndexArray()       {Error("GetRowIndexArray","Dummy"); return 0;};
  63   virtual       const Int_t      *GetColIndexArray() const {Error("GetColIndexArray","Dummy"); return 0;};
  64   virtual             Int_t      *GetColIndexArray()       {Error("GetColIndexArray","Dummy"); return 0;};
  65   virtual             TMatrixDBase &SetRowIndexArray(Int_t *) {Error("SetRowIndexArray","Dummy"); return *this;}
  66   virtual             TMatrixDBase &SetColIndexArray(Int_t *) {Error("SetColIndexArray","Dummy"); return *this;}
  67   virtual             TMatrixDBase &GetSub(Int_t,Int_t,Int_t,Int_t,TMatrixDBase &,Option_t *) const {Error("GetSub","Dummy"); return *((TMatrixDBase*)this);}
  68   virtual             TMatrixDBase &SetSub(Int_t,Int_t,const TMatrixDBase &) {Error("GetSub","Dummy"); return *this;}
  69   virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
  70   virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
  71   //
  72   // ----------------------------- Choleski methods ----------------------------------------
  73   AliSymMatrix*       DecomposeChol();                  //Obtain Cholesky decomposition L matrix
  74   void                InvertChol(AliSymMatrix* mchol);  //Invert using provided Choleski decomposition
  75   Bool_t              InvertChol();                     //Invert
  76   Bool_t              SolveChol(Double_t *brhs, Bool_t invert=kFALSE);
  77   Bool_t              SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert=kFALSE);
  78   Bool_t              SolveChol(TVectorD &brhs, Bool_t invert=kFALSE);
  79   Bool_t              SolveChol(TVectorD &brhs, TVectorD& bsol,Bool_t invert=kFALSE);
  80   //
  81   int                 SolveSpmInv(double *vecB, Bool_t stabilize=kTRUE);
  82
  83  protected:
  84   virtual             Int_t GetIndex(Int_t row,Int_t col)      const;
  85   Double_t            GetEl(Int_t row, Int_t col)              const {return operator()(row,col);}
  86   void                SetEl(Int_t row, Int_t col,Double_t val)       {operator()(row,col) = val;}
  87   //
  88  protected:
  89   Double_t*  fElems;     //   Elements booked by constructor
  90   Double_t** fElemsAdd;  //   Elements (rows) added dynamicaly
  91   //
  92   static AliSymMatrix* fgBuffer;  // buffer for fast solution
  93   static Int_t         fgCopyCnt;  // matrix copy counter
  94   ClassDef(AliSymMatrix,0) //Symmetric Matrix Class
  95 };
  96
  97
  98 //___________________________________________________________
  99 inline Int_t AliSymMatrix::GetIndex(Int_t row,Int_t col) const
 100 {
 101   // lower triangle is actually filled
 102   return ((row*(row+1))>>1) + col;
 103 }
 104
 105 //___________________________________________________________
 106 inline Double_t AliSymMatrix::operator()(Int_t row, Int_t col) const
 107 {
 108   //
 109   if (row<col) Swap(row,col);
 110   if (row>=fNrowIndex) return 0;
 111   return (const Double_t&) (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
 112 }
 113
 114 //___________________________________________________________
 115 inline Double_t& AliSymMatrix::operator()(Int_t row, Int_t col)
 116 {
 117   if (row<col) Swap(row,col);
 118   if (row>=fNrowIndex) AddRows(row-fNrowIndex+1);
 119   return (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
 120 }
 121
 122 //___________________________________________________________
 123 inline void AliSymMatrix::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const
 124 {
 125   MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
 126 }
 127
 128 //___________________________________________________________
 129 inline void AliSymMatrix::Scale(Double_t coeff)
 130 {
 131   for (int i=fNrowIndex;i--;) for (int j=i;j--;) { double& el = operator()(i,j); if (el) el *= coeff;}
 132 }
 133
 134 //___________________________________________________________
 135 inline void AliSymMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
 136 {
 137   for (int i=n;i--;) (*this)(indc[i],r) += valc[i];
 138 }
 139
 140 #endif