1 #ifndef ALISYMBDMATRIX_H
2 #define ALISYMBDMATRIX_H
4 /*********************************************************************************/
5 /* Symmetric Band Diagonal matrix with half band width W (+1 for diagonal) */
6 /* Only lower triangle is stored in the "profile" format */
9 /* Author: ruben.shahoyan@cern.ch */
11 /*********************************************************************************/
16 #include "AliMatrixSq.h"
19 class AliSymBDMatrix : public AliMatrixSq {
22 enum {kDecomposedBit = 0x1};
25 AliSymBDMatrix(Int_t size, Int_t w = 0);
26 AliSymBDMatrix(const AliSymBDMatrix &mat);
27 virtual ~AliSymBDMatrix();
29 Int_t GetBandHWidth() const {return fNrows;}
30 Int_t GetNElemsStored() const {return fNelems;}
31 void Clear(Option_t* option="");
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);
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);}
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());}
49 void Print(Option_t* option="") const;
50 void SetDecomposed(Bool_t v=kTRUE) {SetBit(kDecomposedBit,v);}
51 Bool_t IsDecomposed() const {return TestBit(kDecomposedBit);}
53 void MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const;
54 void MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const;
55 void AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
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;}
64 Double_t* fElems; // Elements booked by constructor
66 ClassDef(AliSymBDMatrix,0) //Symmetric Matrix Class
70 //___________________________________________________________
71 inline Int_t AliSymBDMatrix::GetIndex(Int_t row,Int_t col) const
73 // lower triangle band is actually filled
74 if (row<col) Swap(row,col);
76 if (col < -GetBandHWidth()) return -1;
77 return GetIndex(row) + col;
80 //___________________________________________________________
81 inline Int_t AliSymBDMatrix::GetIndex(Int_t diagID) const
83 // Get index of the diagonal element on row diagID
84 return (diagID+1)*fRowLwb-1;
87 //___________________________________________________________
88 inline Double_t AliSymBDMatrix::operator()(Int_t row, Int_t col) const
91 int idx = GetIndex(row,col);
92 return (const Double_t&) idx<0 ? 0.0 : fElems[idx];
95 //___________________________________________________________
96 inline Double_t& AliSymBDMatrix::operator()(Int_t row, Int_t col)
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];
105 //___________________________________________________________
106 inline Double_t AliSymBDMatrix::operator()(Int_t row) const
109 return (const Double_t&) fElems[GetIndex(row)];
112 //___________________________________________________________
113 inline Double_t& AliSymBDMatrix::operator()(Int_t row)
115 // get diagonal for assingment; assignment outside of the stored range has no effect
116 return fElems[GetIndex(row)];
119 //___________________________________________________________
120 inline void AliSymBDMatrix::MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const
122 MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());