1 #ifndef ALIMATRIXSPARSE_H
2 #define ALIMATRIXSPARSE_H
4 #include "AliMatrixSq.h"
5 #include "AliVectorSparse.h"
8 /**********************************************************************************************/
9 /* Sparse matrix class, used as a global matrix for AliMillePede2 */
11 /* Author: ruben.shahoyan@cern.ch */
13 /**********************************************************************************************/
17 class AliMatrixSparse : public AliMatrixSq
20 AliMatrixSparse() : fVecs(0) {}
21 AliMatrixSparse(Int_t size);
22 AliMatrixSparse(const AliMatrixSparse &mat);
23 virtual ~AliMatrixSparse() {Clear();}
25 AliVectorSparse* GetRow(Int_t ir) const {return (ir<fNcols) ? fVecs[ir] : 0;}
26 AliVectorSparse* GetRowAdd(Int_t ir);
28 virtual Int_t GetSize() const {return fNrows;}
29 virtual Int_t GetNRows() const {return fNrows;}
30 virtual Int_t GetNCols() const {return fNcols;}
32 void Clear(Option_t* option="");
33 void Reset() {for (int i=fNcols;i--;) GetRow(i)->Reset();}
34 void Print(Option_t* option="") const;
35 AliMatrixSparse& operator=(const AliMatrixSparse& src);
36 Double_t& operator()(Int_t row,Int_t col);
37 Double_t operator()(Int_t row,Int_t col) const;
38 void SetToZero(Int_t row,Int_t col);
39 Float_t GetDensity() const;
41 Double_t DiagElem(Int_t r) const;
42 Double_t& DiagElem(Int_t r);
43 void SortIndices(Bool_t valuesToo=kFALSE);
45 void MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const;
46 void MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const;
48 void AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n);
52 AliVectorSparse** fVecs;
54 ClassDef(AliMatrixSparse,0)
57 //___________________________________________________
58 inline void AliMatrixSparse::MultiplyByVec(const TVectorD &vecIn, TVectorD &vecOut) const
61 MultiplyByVec((Double_t*)vecIn.GetMatrixArray(),(Double_t*)vecOut.GetMatrixArray());
64 //___________________________________________________
65 inline void AliMatrixSparse::SetToZero(Int_t row,Int_t col)
67 // set existing element to 0
68 if (IsSymmetric() && col>row) Swap(row,col);
69 AliVectorSparse* rowv = GetRow(row);
70 if (rowv) rowv->SetToZero(col);
73 //___________________________________________________
74 inline Double_t AliMatrixSparse::operator()(Int_t row,Int_t col) const
76 // printf("M: find\n");
77 if (IsSymmetric() && col>row) Swap(row,col);
78 AliVectorSparse* rowv = GetRow(row);
80 return rowv->FindIndex(col);
83 //___________________________________________________
84 inline Double_t& AliMatrixSparse::operator()(Int_t row,Int_t col)
86 // printf("M: findindexAdd\n");
87 if (IsSymmetric() && col>row) Swap(row,col);
88 AliVectorSparse* rowv = GetRowAdd(row);
89 if (col>=fNcols) fNcols = col+1;
90 return rowv->FindIndexAdd(col);
93 //___________________________________________________
94 inline Double_t AliMatrixSparse::DiagElem(Int_t row) const
97 AliVectorSparse* rowv = GetRow(row);
99 if (IsSymmetric()) return (rowv->GetNElems()>0 && rowv->GetLastIndex()==row) ? rowv->GetLastElem() : 0.;
100 else return rowv->FindIndex(row);
104 //___________________________________________________
105 inline Double_t &AliMatrixSparse::DiagElem(Int_t row)
108 AliVectorSparse* rowv = GetRowAdd(row);
109 if (row>=fNcols) fNcols = row+1;
111 return (rowv->GetNElems()>0 && rowv->GetLastIndex()==row) ?
112 rowv->GetLastElem() : rowv->FindIndexAdd(row);
114 else return rowv->FindIndexAdd(row);