Update for changes in Root TExMap, TExMapIter in trunk
[u/mrichter/AliRoot.git] / STEER / AliSymMatrix.h
CommitLineData
8a9ab0eb 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
11class 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 //
de34b538 31 Double_t* GetRow(Int_t r);
32 //
8a9ab0eb 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;
de34b538 39 void AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
8a9ab0eb 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//___________________________________________________________
82inline 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//___________________________________________________________
89inline 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//___________________________________________________________
98inline 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//___________________________________________________________
106inline void AliSymMatrix::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const
107{
108 MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
109}
110
111//___________________________________________________________
112inline 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
de34b538 117//___________________________________________________________
118inline 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}
8a9ab0eb 122
123#endif