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   virtual Int_t   GetSize()                            const {return fNrows;}
  29   virtual Int_t   GetNRows()                           const {return fNrows;}
  30   virtual Int_t   GetNCols()                           const {return fNcols;}
  31   //
  32   void Clear(Option_t* option="");
  33   void Reset()                            {for (int i=fNcols;i--;) GetRow(i)->Reset();}
  34   void Print(Option_t* option="")                      const;
  35   AliMatrixSparse& operator=(const AliMatrixSparse& src);
  36   Double_t&        operator()(Int_t row,Int_t col);
  37   Double_t         operator()(Int_t row,Int_t col)     const;
  38   void             SetToZero(Int_t row,Int_t col);
  39   Float_t          GetDensity()                        const;
  40   //
  41   Double_t         DiagElem(Int_t r)                   const;
  42   Double_t&        DiagElem(Int_t r);
  43   void             SortIndices(Bool_t valuesToo=kFALSE);
  44   //
  45   void MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const;
  46   void MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const;
  47   //
  48   void AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
  49   //
  50  protected:
  51   //
  52   AliVectorSparse** fVecs;
  53   //
  54   ClassDef(AliMatrixSparse,0)
  55 };
  56
  57 //___________________________________________________
  58 inline void AliMatrixSparse::MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const
  59 {
  60   MultiplyByVec((Double_t*)vecIn.GetMatrixArray(),(Double_t*)vecOut.GetMatrixArray());
  61 }
  62
  63 //___________________________________________________
  64 inline void AliMatrixSparse::SetToZero(Int_t row,Int_t col)
  65 {
  66   //  set existing element to 0
  67   if (IsSymmetric() && col>row) Swap(row,col);
  68   AliVectorSparse* rowv = GetRow(row);
  69   if (rowv) rowv->SetToZero(col);
  70 }
  71
  72 //___________________________________________________
  73 inline Double_t AliMatrixSparse::operator()(Int_t row,Int_t col) const
  74 {
  75   //  printf("M: find\n");
  76   if (IsSymmetric() && col>row) Swap(row,col);
  77   AliVectorSparse* rowv = GetRow(row);
  78   if (!rowv) return 0;
  79   return rowv->FindIndex(col);
  80 }
  81
  82 //___________________________________________________
  83 inline Double_t& AliMatrixSparse::operator()(Int_t row,Int_t col)
  84 {
  85   //  printf("M: findindexAdd\n");
  86   if (IsSymmetric() && col>row) Swap(row,col);
  87   AliVectorSparse* rowv = GetRowAdd(row);
  88   if (col>=fNcols) fNcols = col+1;
  89   return rowv->FindIndexAdd(col);
  90 }
  91
  92 //___________________________________________________
  93 inline Double_t AliMatrixSparse::DiagElem(Int_t row) const
  94 {
  95   AliVectorSparse* rowv = GetRow(row);
  96   if (!rowv) return 0;
  97   if (IsSymmetric()) return (rowv->GetNElems()>0 && rowv->GetLastIndex()==row) ? rowv->GetLastElem() : 0.;
  98   else return rowv->FindIndex(row);
  99   //
 100 }
 101
 102 //___________________________________________________
 103 inline Double_t &AliMatrixSparse::DiagElem(Int_t row)
 104 {
 105   AliVectorSparse* rowv = GetRowAdd(row);
 106   if (row>=fNcols) fNcols = row+1;
 107   if (IsSymmetric()) {
 108     return (rowv->GetNElems()>0 && rowv->GetLastIndex()==row) ?
 109                        rowv->GetLastElem() : rowv->FindIndexAdd(row);
 110   }
 111   else return rowv->FindIndexAdd(row);
 112   //
 113 }
 114
 115
 116 #endif
 117