STEER/AliMatrixSparse.h

   1 #ifndef ALIMATRIXSPARSE_H
   2 #define ALIMATRIXSPARSE_H
   3
   4 #include "AliMatrixSq.h"
   5 #include "AliVectorSparse.h"
   6
   7
   8 /**********************************************************************************************/
   9 /* Sparse matrix class, used as a global matrix for AliMillePede2                             */
  10 /*                                                                                            */
  11 /* Author: ruben.shahoyan@cern.ch                                                             */
  12 /*                                                                                            */
  13 /**********************************************************************************************/
  14
  15
  16 //
  17 class AliMatrixSparse : public AliMatrixSq
  18 {
  19  public:
  20   AliMatrixSparse() : fVecs(0) {}
  21   AliMatrixSparse(Int_t size);
  22   AliMatrixSparse(const AliMatrixSparse &mat);
  23   virtual ~AliMatrixSparse() {Clear();}
  24   //
  25   AliVectorSparse* GetRow(Int_t ir)        const {return (ir<fNcols) ? fVecs[ir] : 0;}
  26   AliVectorSparse* GetRowAdd(Int_t ir);
  27
  28   void Clear(Option_t* option="");
  29   void Reset()                            {for (int i=fNcols;i--;) GetRow(i)->Reset();}
  30   void Print(Option_t* option="")                      const;
  31   AliMatrixSparse& operator=(const AliMatrixSparse& src);
  32   Double_t&        operator()(Int_t row,Int_t col);
  33   Double_t         operator()(Int_t row,Int_t col)     const;
  34   void             SetToZero(Int_t row,Int_t col);
  35   Float_t          GetDensity()                        const;
  36   //
  37   Double_t         DiagElem(Int_t r)                   const;
  38   Double_t&        DiagElem(Int_t r);
  39   void             SortIndices(Bool_t valuesToo=kFALSE);
  40   //
  41   void MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const;
  42   void MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
  43   //
  44   void AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
  45   //
  46  protected:
  47   //
  48   AliVectorSparse** fVecs;
  49   //
  50   ClassDef(AliMatrixSparse,0)
  51 };
  52
  53 //___________________________________________________
  54 inline void AliMatrixSparse::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const
  55 {
  56   MultiplyByVec((Double_t*)vecIn.GetMatrixArray(),(Double_t*)vecOut.GetMatrixArray());
  57 }
  58
  59 //___________________________________________________
  60 inline void AliMatrixSparse::SetToZero(Int_t row,Int_t col)
  61 {
  62   //  set existing element to 0
  63   if (IsSymmetric() && col>row) Swap(row,col);
  64   AliVectorSparse* rowv = GetRow(row);
  65   if (rowv) rowv->SetToZero(col);
  66 }
  67
  68 //___________________________________________________
  69 inline Double_t AliMatrixSparse::operator()(Int_t row,Int_t col) const
  70 {
  71   //  printf("M: find\n");
  72   if (IsSymmetric() && col>row) Swap(row,col);
  73   AliVectorSparse* rowv = GetRow(row);
  74   if (!rowv) return 0;
  75   return rowv->FindIndex(col);
  76 }
  77
  78 //___________________________________________________
  79 inline Double_t& AliMatrixSparse::operator()(Int_t row,Int_t col)
  80 {
  81   //  printf("M: findindexAdd\n");
  82   if (IsSymmetric() && col>row) Swap(row,col);
  83   AliVectorSparse* rowv = GetRowAdd(row);
  84   return rowv->FindIndexAdd(col);
  85 }
  86
  87 //___________________________________________________
  88 inline Double_t AliMatrixSparse::DiagElem(Int_t row) const
  89 {
  90   AliVectorSparse* rowv = GetRow(row);
  91   if (!rowv) return 0;
  92   if (IsSymmetric()) return (rowv->GetNElems()>0 && rowv->GetLastIndex()==row) ? rowv->GetLastElem() : 0.;
  93   else return rowv->FindIndex(row);
  94   //
  95 }
  96
  97 //___________________________________________________
  98 inline Double_t &AliMatrixSparse::DiagElem(Int_t row)
  99 {
 100   AliVectorSparse* rowv = GetRowAdd(row);
 101   if (IsSymmetric()) return (rowv->GetNElems()>0 && rowv->GetLastIndex()==row) ?
 102                        rowv->GetLastElem() : rowv->FindIndexAdd(row);
 103   else return rowv->FindIndexAdd(row);
 104   //
 105 }
 106
 107
 108 #endif
 109