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                         */
7 /*                                                                               */
8 /*                                                                               */
9 /* Author: ruben.shahoyan@cern.ch                                                */
10 /*                                                                               */
11 /*********************************************************************************/
13 //#include <string.h>
14 #include <TObject.h>
15 #include <TVectorD.h>
16 #include "AliMatrixSq.h"
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;
50   void          SetDecomposed()                                        {SetBit(kDecomposedBit);}
51   Bool_t        IsDecomposed()                                   const {return TestBit(kDecomposedBit);}
52   //
53   void          MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const;
54   void          MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const;
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 };
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 }
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 }
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 }
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 }
105 //___________________________________________________________
106 inline Double_t AliSymBDMatrix::operator()(Int_t row) const
107 {
108   // query diagonal
109   return (const Double_t&) fElems[GetIndex(row)];
110 }
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 }
119 //___________________________________________________________
120 inline void AliSymBDMatrix::MultiplyByVec(TVectorD &vecIn, TVectorD &vecOut) const
121 {
122   MultiplyByVec(vecIn.GetMatrixArray(), vecOut.GetMatrixArray());
123 }
125 #endif