STEER/STEER/AliSymBDMatrix.h

   1 #ifndef ALISYMBDMATRIX_H
   2 #define ALISYMBDMATRIX_H
   3
   4 /*********************************************************************************/
   5 /* Symmetric Band Diagonal matrix with half band width W (+1 for diagonal)       */
   6 /* Only lower triangle is stored in the "profile" format                         */
   7 /*                                                                               */
   8 /*                                                                               */
   9 /* Author: ruben.shahoyan@cern.ch                                                */
  10 /*                                                                               */
  11 /*********************************************************************************/
  12
  13 //#include <string.h>
  14 #include <TObject.h>
  15 #include <TVectorD.h>
  16 #include "AliMatrixSq.h"
  17
  18
  19 class AliSymBDMatrix : public AliMatrixSq {
  20   //
  21  public:
  22   enum {kDecomposedBit = 0x1};
  23   //
  24   AliSymBDMatrix();
  25   AliSymBDMatrix(Int_t size, Int_t w = 0);
  26   AliSymBDMatrix(const AliSymBDMatrix &mat);
  27   virtual ~AliSymBDMatrix();
  28   //
  29   Int_t         GetBandHWidth()                                  const {return fNrows;}
  30   Int_t         GetNElemsStored()                                const {return fNelems;}
  31   void          Clear(Option_t* option="");
  32   void          Reset();
  33   //
  34   Float_t       GetDensity()                                     const;
  35   AliSymBDMatrix& operator=(const AliSymBDMatrix& src);
  36   Double_t      operator()(Int_t rown, Int_t coln)               const;
  37   Double_t&     operator()(Int_t rown, Int_t coln);
  38   Double_t      operator()(Int_t rown)                           const;
  39   Double_t&     operator()(Int_t rown);
  40   //
  41   Double_t      DiagElem(Int_t r)                                const {return (*(const AliSymBDMatrix*)this)(r,r);}
  42   Double_t&     DiagElem(Int_t r)                                      {return (*this)(r,r);}
  43   void          DecomposeLDLT();
  44   void          Solve(Double_t *rhs);
  45   void          Solve(const Double_t *rhs,Double_t *sol);
  46   void          Solve(TVectorD &rhs)                                   {Solve(rhs.GetMatrixArray());}
  47   void          Solve(const TVectorD &rhs,TVectorD &sol)               {Solve(rhs.GetMatrixArray(),sol.GetMatrixArray());}
  48   //
  49   void          Print(Option_t* option="")                       const;
  50   void          SetDecomposed(Bool_t v=kTRUE)                          {SetBit(kDecomposedBit,v);}
  51   Bool_t        IsDecomposed()                                   const {return TestBit(kDecomposedBit);}
  52   //
  53   void          MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const;
  54   void          MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const;
  55   void          AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
  56   //
  57   // protected:
  58   virtual             Int_t GetIndex(Int_t row,Int_t col)      const;
  59   virtual             Int_t GetIndex(Int_t diagID)             const;
  60   Double_t            GetEl(Int_t row, Int_t col)              const {return operator()(row,col);}
  61   void                SetEl(Int_t row, Int_t col,Double_t val)       {operator()(row,col) = val;}
  62   //
  63  protected:
  64   Double_t*  fElems;     //   Elements booked by constructor
  65   //
  66   ClassDef(AliSymBDMatrix,0) //Symmetric Matrix Class
  67 };
  68
  69
  70 //___________________________________________________________
  71 inline Int_t AliSymBDMatrix::GetIndex(Int_t row,Int_t col) const
  72 {
  73   // lower triangle band is actually filled
  74   if (row<col) Swap(row,col);
  75   col -= row;
  76   if (col < -GetBandHWidth()) return -1;
  77   return GetIndex(row) + col;
  78 }
  79
  80 //___________________________________________________________
  81 inline Int_t AliSymBDMatrix::GetIndex(Int_t diagID) const
  82 {
  83   // Get index of the diagonal element on row diagID
  84   return (diagID+1)*fRowLwb-1;
  85 }
  86
  87 //___________________________________________________________
  88 inline Double_t AliSymBDMatrix::operator()(Int_t row, Int_t col) const
  89 {
  90   // query element
  91   int idx = GetIndex(row,col);
  92   return (const Double_t&) idx<0 ? 0.0 : fElems[idx];
  93 }
  94
  95 //___________________________________________________________
  96 inline Double_t& AliSymBDMatrix::operator()(Int_t row, Int_t col)
  97 {
  98   // get element for assingment; assignment outside of the stored range has no effect
  99   int idx = GetIndex(row,col);
 100   if (idx>=0) return fElems[idx];
 101   fTol = 0;
 102   return fTol;
 103 }
 104
 105 //___________________________________________________________
 106 inline Double_t AliSymBDMatrix::operator()(Int_t row) const
 107 {
 108   // query diagonal
 109   return (const Double_t&) fElems[GetIndex(row)];
 110 }
 111
 112 //___________________________________________________________
 113 inline Double_t& AliSymBDMatrix::operator()(Int_t row)
 114 {
 115   // get diagonal for assingment; assignment outside of the stored range has no effect
 116   return fElems[GetIndex(row)];
 117 }
 118
 119 //___________________________________________________________
 120 inline void AliSymBDMatrix::MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const
 121 {
 122   MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
 123 }
 124
 125 #endif