1 #include "AliMatrixSparse.h"
2 #include "TStopwatch.h"
4 /**********************************************************************************************/
5 /* Sparse matrix class, used as a global matrix for AliMillePede2 */
7 /* Author: ruben.shahoyan@cern.ch */
9 /**********************************************************************************************/
11 //___________________________________________________________
12 ClassImp(AliMatrixSparse)
14 //___________________________________________________________
15 AliMatrixSparse::AliMatrixSparse(Int_t sz)
16 : AliMatrixSq(),fVecs(0)
20 fVecs = new AliVectorSparse*[sz];
21 for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse();
24 //___________________________________________________________
25 AliMatrixSparse::AliMatrixSparse(const AliMatrixSparse& src)
26 : AliMatrixSq(src),fVecs(0)
28 fVecs = new AliVectorSparse*[fNcols];
29 for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
32 //___________________________________________________________
33 AliVectorSparse* AliMatrixSparse::GetRowAdd(Int_t ir)
36 AliVectorSparse** arrv = new AliVectorSparse*[ir+1];
37 for (int i=GetSize();i--;) arrv[i] = fVecs[i];
40 for (int i=GetSize();i<=ir;i++) fVecs[i] = new AliVectorSparse();
41 fNcols = fNrows = ir+1;
46 //___________________________________________________________
47 AliMatrixSparse& AliMatrixSparse::operator=(const AliMatrixSparse& src)
49 if (*this == src) return *this;
51 fNcols = fNrows = src.GetSize();
52 SetSymmetric(src.IsSymmetric());
53 fVecs = new AliVectorSparse*[fNcols];
54 for (int i=fNcols;i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
58 //___________________________________________________________
59 void AliMatrixSparse::Clear(Option_t*)
61 for (int i=fNcols;i--;) delete GetRow(i);
66 //___________________________________________________________
67 void AliMatrixSparse::Print(Option_t*) const
69 printf("Sparse Matrix of size %d %s\n",fNcols,IsSymmetric() ? " (Symmetric)":"");
70 for (int i=0;i<fNcols;i++) {
71 AliVectorSparse* row = GetRow(i);
72 if (!row->GetNElems()) continue;
78 //___________________________________________________________
79 void AliMatrixSparse::MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const
81 // fill vecOut by matrix*vecIn
82 // vector should be of the same size as the matrix
84 memset(vecOut,0,GetSize()*sizeof(Double_t));
86 for (int rw=GetSize();rw--;) { // loop over rows >>>
87 const AliVectorSparse* rowV = GetRow(rw);
88 Int_t nel = rowV->GetNElems();
91 UShort_t* indV = rowV->GetIndices();
92 Double_t* elmV = rowV->GetElems();
95 // treat diagonal term separately. If filled, it should be the last one
96 if (indV[--nel]==rw) vecOut[rw] += vecIn[rw]*elmV[nel];
97 else nel = rowV->GetNElems(); // diag elem was not filled
99 for (int iel=nel;iel--;) { // less element retrieval for symmetric case
101 vecOut[rw] += vecIn[indV[iel]]*elmV[iel];
102 vecOut[indV[iel]] += vecIn[rw]*elmV[iel];
106 else for (int iel=nel;iel--;) if (elmV[iel]) vecOut[rw] += vecIn[indV[iel]]*elmV[iel];
108 } // loop over rows <<<
112 //___________________________________________________________
113 void AliMatrixSparse::SortIndices(Bool_t valuesToo)
117 printf("AliMatrixSparse:sort>>\n");
118 // sort columns in increasing order. Used to fix the matrix after ILUk decompostion
119 for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo);
122 printf("AliMatrixSparse:sort<<\n");
125 //___________________________________________________________
126 void AliMatrixSparse::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
128 // for sym. matrix count how many elems to add have row>=col and assign excplicitly
129 // those which have row<col
131 // range in increasing order of indices
133 for (int j=i;j>=0;j--) {
134 if (indc[j]>indc[i]) { // swap
135 int ti = indc[i]; indc[i] = indc[j]; indc[j] = ti;
136 double tv = valc[i]; valc[i] = valc[j]; valc[j] = tv;
144 if (indc[ni]>r) (*this)(indc[ni],r) += valc[ni];
145 else break; // use the fact that the indices are ranged in increasing order
150 AliVectorSparse* row = GetRowAdd(r);
151 row->Add(valc,indc,ni+1);
154 //___________________________________________________________
155 Float_t AliMatrixSparse::GetDensity() const
157 // get fraction of non-zero elements
159 for (int i=GetSize();i--;) nel += GetRow(i)->GetNElems();
160 int den = IsSymmetric() ? (GetSize()+1)*GetSize()/2 : GetSize()*GetSize();
161 return float(nel)/den;