]>
Commit | Line | Data |
---|---|---|
7c3070ec | 1 | #ifndef ALISYMBDMATRIX_H |
2 | #define ALISYMBDMATRIX_H | |
3 | ||
4 | /*********************************************************************************/ | |
5 | /* Symmetric Band Diagonal matrix with half band width W (+1 for diagonal) */ | |
6 | /* Only lower triangle is stored in the "profile" format */ | |
7 | /* */ | |
8 | /* */ | |
9 | /* Author: ruben.shahoyan@cern.ch */ | |
10 | /* */ | |
11 | /*********************************************************************************/ | |
12 | ||
13 | //#include <string.h> | |
14 | #include <TObject.h> | |
15 | #include <TVectorD.h> | |
16 | #include "AliMatrixSq.h" | |
17 | ||
18 | ||
19 | class AliSymBDMatrix : public AliMatrixSq { | |
20 | // | |
21 | public: | |
22 | enum {kDecomposedBit = 0x1}; | |
23 | // | |
24 | AliSymBDMatrix(); | |
25 | AliSymBDMatrix(Int_t size, Int_t w = 0); | |
26 | AliSymBDMatrix(const AliSymBDMatrix &mat); | |
27 | virtual ~AliSymBDMatrix(); | |
28 | // | |
29 | Int_t GetBandHWidth() const {return fNrows;} | |
30 | Int_t GetNElemsStored() const {return fNelems;} | |
31 | void Clear(Option_t* option=""); | |
32 | void Reset(); | |
33 | // | |
34 | Float_t GetDensity() const; | |
35 | AliSymBDMatrix& operator=(const AliSymBDMatrix& src); | |
36 | Double_t operator()(Int_t rown, Int_t coln) const; | |
37 | Double_t& operator()(Int_t rown, Int_t coln); | |
38 | Double_t operator()(Int_t rown) const; | |
39 | Double_t& operator()(Int_t rown); | |
40 | // | |
41 | Double_t DiagElem(Int_t r) const {return (*(const AliSymBDMatrix*)this)(r,r);} | |
42 | Double_t& DiagElem(Int_t r) {return (*this)(r,r);} | |
43 | void DecomposeLDLT(); | |
44 | void Solve(Double_t *rhs); | |
45 | void Solve(const Double_t *rhs,Double_t *sol); | |
46 | void Solve(TVectorD &rhs) {Solve(rhs.GetMatrixArray());} | |
47 | void Solve(const TVectorD &rhs,TVectorD &sol) {Solve(rhs.GetMatrixArray(),sol.GetMatrixArray());} | |
48 | // | |
49 | void Print(Option_t* option="") const; | |
c130d912 | 50 | void SetDecomposed(Bool_t v=kTRUE) {SetBit(kDecomposedBit,v);} |
7c3070ec | 51 | Bool_t IsDecomposed() const {return TestBit(kDecomposedBit);} |
52 | // | |
551c9e69 | 53 | void MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const; |
54 | void MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const; | |
7c3070ec | 55 | void AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n); |
56 | // | |
57 | // protected: | |
58 | virtual Int_t GetIndex(Int_t row,Int_t col) const; | |
59 | virtual Int_t GetIndex(Int_t diagID) const; | |
60 | Double_t GetEl(Int_t row, Int_t col) const {return operator()(row,col);} | |
61 | void SetEl(Int_t row, Int_t col,Double_t val) {operator()(row,col) = val;} | |
62 | // | |
63 | protected: | |
64 | Double_t* fElems; // Elements booked by constructor | |
65 | // | |
66 | ClassDef(AliSymBDMatrix,0) //Symmetric Matrix Class | |
67 | }; | |
68 | ||
69 | ||
70 | //___________________________________________________________ | |
71 | inline Int_t AliSymBDMatrix::GetIndex(Int_t row,Int_t col) const | |
72 | { | |
73 | // lower triangle band is actually filled | |
74 | if (row<col) Swap(row,col); | |
75 | col -= row; | |
76 | if (col < -GetBandHWidth()) return -1; | |
77 | return GetIndex(row) + col; | |
78 | } | |
79 | ||
80 | //___________________________________________________________ | |
81 | inline Int_t AliSymBDMatrix::GetIndex(Int_t diagID) const | |
82 | { | |
83 | // Get index of the diagonal element on row diagID | |
84 | return (diagID+1)*fRowLwb-1; | |
85 | } | |
86 | ||
87 | //___________________________________________________________ | |
88 | inline Double_t AliSymBDMatrix::operator()(Int_t row, Int_t col) const | |
89 | { | |
90 | // query element | |
91 | int idx = GetIndex(row,col); | |
92 | return (const Double_t&) idx<0 ? 0.0 : fElems[idx]; | |
93 | } | |
94 | ||
95 | //___________________________________________________________ | |
96 | inline Double_t& AliSymBDMatrix::operator()(Int_t row, Int_t col) | |
97 | { | |
98 | // get element for assingment; assignment outside of the stored range has no effect | |
99 | int idx = GetIndex(row,col); | |
100 | if (idx>=0) return fElems[idx]; | |
101 | fTol = 0; | |
102 | return fTol; | |
103 | } | |
104 | ||
105 | //___________________________________________________________ | |
106 | inline Double_t AliSymBDMatrix::operator()(Int_t row) const | |
107 | { | |
108 | // query diagonal | |
109 | return (const Double_t&) fElems[GetIndex(row)]; | |
110 | } | |
111 | ||
112 | //___________________________________________________________ | |
113 | inline Double_t& AliSymBDMatrix::operator()(Int_t row) | |
114 | { | |
115 | // get diagonal for assingment; assignment outside of the stored range has no effect | |
116 | return fElems[GetIndex(row)]; | |
117 | } | |
118 | ||
119 | //___________________________________________________________ | |
551c9e69 | 120 | inline void AliSymBDMatrix::MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const |
7c3070ec | 121 | { |
122 | MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray()); | |
123 | } | |
124 | ||
125 | #endif |