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