f22d5d54ea227ebbda942cd31956a4358dfc7fd4
[u/mrichter/AliRoot.git] / 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   Double_t*     GetRow(Int_t r);
32   //
33   void          Print(Option_t* option="")                       const;
34   void          AddRows(int nrows=1);
35   //
36   void          Scale(Double_t coeff);
37   void          MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
38   void          MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const;
39   void          AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
40   //
41   // ---------------------------------- Dummy methods of MatrixBase
42   virtual       const Double_t   *GetMatrixArray  () const {return fElems;};
43   virtual             Double_t   *GetMatrixArray  ()       {return (Double_t*)fElems;}
44   virtual       const Int_t      *GetRowIndexArray() const {Error("GetRowIndexArray","Dummy"); return 0;};
45   virtual             Int_t      *GetRowIndexArray()       {Error("GetRowIndexArray","Dummy"); return 0;};
46   virtual       const Int_t      *GetColIndexArray() const {Error("GetColIndexArray","Dummy"); return 0;};
47   virtual             Int_t      *GetColIndexArray()       {Error("GetColIndexArray","Dummy"); return 0;};
48   virtual             TMatrixDBase &SetRowIndexArray(Int_t *) {Error("SetRowIndexArray","Dummy"); return *this;}
49   virtual             TMatrixDBase &SetColIndexArray(Int_t *) {Error("SetColIndexArray","Dummy"); return *this;}
50   virtual             TMatrixDBase &GetSub(Int_t,Int_t,Int_t,Int_t,TMatrixDBase &,Option_t *) const {Error("GetSub","Dummy"); return *((TMatrixDBase*)this);}
51   virtual             TMatrixDBase &SetSub(Int_t,Int_t,const TMatrixDBase &) {Error("GetSub","Dummy"); return *this;}
52   virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
53   virtual             TMatrixDBase &ResizeTo (Int_t,Int_t,Int_t,Int_t,Int_t) {Error("ResizeTo","Dummy"); return *this;}
54   //
55   // ----------------------------- Choleski methods ----------------------------------------
56   AliSymMatrix*       DecomposeChol();                  //Obtain Cholesky decomposition L matrix 
57   void                InvertChol(AliSymMatrix* mchol);  //Invert using provided Choleski decomposition
58   Bool_t              InvertChol();                     //Invert 
59   Bool_t              SolveChol(Double_t *brhs, Bool_t invert=kFALSE);
60   Bool_t              SolveChol(Double_t *brhs, Double_t *bsol,Bool_t invert=kFALSE);
61   Bool_t              SolveChol(TVectorD &brhs, Bool_t invert=kFALSE);
62   Bool_t              SolveChol(TVectorD &brhs, TVectorD& bsol,Bool_t invert=kFALSE);
63   //
64   int                 SolveSpmInv(double *vecB, Bool_t stabilize=kTRUE);
65
66  protected:
67   virtual             Int_t GetIndex(Int_t row,Int_t col)      const;
68   Double_t            GetEl(Int_t row, Int_t col)              const {return operator()(row,col);}
69   void                SetEl(Int_t row, Int_t col,Double_t val)       {operator()(row,col) = val;}
70   //
71  protected:
72   Double_t*  fElems;     //   Elements booked by constructor
73   Double_t** fElemsAdd;  //   Elements (rows) added dynamicaly
74   //
75   static AliSymMatrix* fgBuffer;  // buffer for fast solution
76   static Int_t         fgCopyCnt;  // matrix copy counter
77   ClassDef(AliSymMatrix,0) //Symmetric Matrix Class
78 };
79
80
81 //___________________________________________________________
82 inline Int_t AliSymMatrix::GetIndex(Int_t row,Int_t col) const
83 {
84   // lower triangle is actually filled
85   return ((row*(row+1))>>1) + col;
86 }
87
88 //___________________________________________________________
89 inline Double_t AliSymMatrix::operator()(Int_t row, Int_t col) const
90 {
91   //
92   if (row<col) Swap(row,col);
93   if (row>=fNrowIndex) return 0;
94   return (const Double_t&) (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
95 }
96
97 //___________________________________________________________
98 inline Double_t& AliSymMatrix::operator()(Int_t row, Int_t col)
99 {
100   if (row<col) Swap(row,col);
101   if (row>=fNrowIndex) AddRows(row-fNrowIndex+1);
102   return (row<fNcols ? fElems[GetIndex(row,col)] : (fElemsAdd[row-fNcols])[col]);
103 }
104
105 //___________________________________________________________
106 inline void AliSymMatrix::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const
107 {
108   MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
109 }
110
111 //___________________________________________________________
112 inline void AliSymMatrix::Scale(Double_t coeff)
113 {
114   for (int i=fNrowIndex;i--;) for (int j=i;j--;) { double& el = operator()(i,j); if (el) el *= coeff;}
115 }
116
117 //___________________________________________________________
118 inline void AliSymMatrix::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
119 {
120   for (int i=n;i--;) (*this)(indc[i],r) += valc[i];
121 }
122
123 #endif