]>
Commit | Line | Data |
---|---|---|
1 | #include "AliMatrixSparse.h" | |
2 | #include "TStopwatch.h" | |
3 | ||
4 | /**********************************************************************************************/ | |
5 | /* Sparse matrix class, used as a global matrix for AliMillePede2 */ | |
6 | /* */ | |
7 | /* Author: ruben.shahoyan@cern.ch */ | |
8 | /* */ | |
9 | /* */ | |
10 | /* */ | |
11 | /**********************************************************************************************/ | |
12 | ||
13 | //___________________________________________________________ | |
14 | ClassImp(AliMatrixSparse) | |
15 | ||
16 | //___________________________________________________________ | |
17 | AliMatrixSparse::AliMatrixSparse(Int_t sz) | |
18 | : AliMatrixSq(),fVecs(0) | |
19 | { | |
20 | // constructor | |
21 | fNcols=fNrows=sz; | |
22 | // | |
23 | fVecs = new AliVectorSparse*[sz]; | |
24 | for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse(); | |
25 | } | |
26 | ||
27 | //___________________________________________________________ | |
28 | AliMatrixSparse::AliMatrixSparse(const AliMatrixSparse& src) | |
29 | : AliMatrixSq(src),fVecs(0) | |
30 | { | |
31 | // copy c-tor | |
32 | fVecs = new AliVectorSparse*[src.GetSize()]; | |
33 | for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i)); | |
34 | } | |
35 | ||
36 | //___________________________________________________________ | |
37 | AliVectorSparse* AliMatrixSparse::GetRowAdd(Int_t ir) | |
38 | { | |
39 | // get row, add if needed | |
40 | if (ir>=fNrows) { | |
41 | AliVectorSparse** arrv = new AliVectorSparse*[ir+1]; | |
42 | for (int i=GetSize();i--;) arrv[i] = fVecs[i]; | |
43 | delete [] fVecs; | |
44 | fVecs = arrv; | |
45 | for (int i=GetSize();i<=ir;i++) fVecs[i] = new AliVectorSparse(); | |
46 | fNrows = ir+1; | |
47 | if (IsSymmetric() && fNcols<fNrows) fNcols = fNrows; | |
48 | } | |
49 | return fVecs[ir]; | |
50 | } | |
51 | ||
52 | //___________________________________________________________ | |
53 | AliMatrixSparse& AliMatrixSparse::operator=(const AliMatrixSparse& src) | |
54 | { | |
55 | // assignment op-r | |
56 | if (this == &src) return *this; | |
57 | AliMatrixSq::operator=(src); | |
58 | ||
59 | Clear(); | |
60 | fNcols = src.GetNCols(); | |
61 | fNrows = src.GetNRows(); | |
62 | SetSymmetric(src.IsSymmetric()); | |
63 | fVecs = new AliVectorSparse*[fNrows]; | |
64 | for (int i=fNrows;i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i)); | |
65 | return *this; | |
66 | } | |
67 | ||
68 | //___________________________________________________________ | |
69 | void AliMatrixSparse::Clear(Option_t*) | |
70 | { | |
71 | // clear | |
72 | for (int i=fNrows;i--;) delete GetRow(i); | |
73 | delete [] fVecs; | |
74 | fNcols = fNrows = 0; | |
75 | } | |
76 | ||
77 | //___________________________________________________________ | |
78 | void AliMatrixSparse::Print(Option_t* opt) const | |
79 | { | |
80 | // print itself | |
81 | printf("Sparse Matrix of size %d x %d %s\n",fNrows,fNcols,IsSymmetric() ? " (Symmetric)":""); | |
82 | for (int i=0;i<fNrows;i++) { | |
83 | AliVectorSparse* row = GetRow(i); | |
84 | if (!row->GetNElems()) continue; | |
85 | printf("%3d: ",i); | |
86 | row->Print(opt); | |
87 | } | |
88 | } | |
89 | ||
90 | //___________________________________________________________ | |
91 | void AliMatrixSparse::MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const | |
92 | { | |
93 | // fill vecOut by matrix*vecIn | |
94 | // vector should be of the same size as the matrix | |
95 | // | |
96 | memset(vecOut,0,GetSize()*sizeof(Double_t)); | |
97 | // | |
98 | for (int rw=GetSize();rw--;) { // loop over rows >>> | |
99 | const AliVectorSparse* rowV = GetRow(rw); | |
100 | Int_t nel = rowV->GetNElems(); | |
101 | if (!nel) continue; | |
102 | // | |
103 | UShort_t* indV = rowV->GetIndices(); | |
104 | Double_t* elmV = rowV->GetElems(); | |
105 | // | |
106 | if (IsSymmetric()) { | |
107 | // treat diagonal term separately. If filled, it should be the last one | |
108 | if (indV[--nel]==rw) vecOut[rw] += vecIn[rw]*elmV[nel]; | |
109 | else nel = rowV->GetNElems(); // diag elem was not filled | |
110 | // | |
111 | for (int iel=nel;iel--;) { // less element retrieval for symmetric case | |
112 | if (elmV[iel]) { | |
113 | vecOut[rw] += vecIn[indV[iel]]*elmV[iel]; | |
114 | vecOut[indV[iel]] += vecIn[rw]*elmV[iel]; | |
115 | } | |
116 | } | |
117 | } | |
118 | else for (int iel=nel;iel--;) if (elmV[iel]) vecOut[rw] += vecIn[indV[iel]]*elmV[iel]; | |
119 | // | |
120 | } // loop over rows <<< | |
121 | // | |
122 | } | |
123 | ||
124 | //___________________________________________________________ | |
125 | void AliMatrixSparse::SortIndices(Bool_t valuesToo) | |
126 | { | |
127 | // sort columns in increasing order. Used to fix the matrix after ILUk decompostion | |
128 | TStopwatch sw; | |
129 | sw.Start(); | |
130 | printf("AliMatrixSparse:sort>>\n"); | |
131 | for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo); | |
132 | sw.Stop(); | |
133 | sw.Print(); | |
134 | printf("AliMatrixSparse:sort<<\n"); | |
135 | } | |
136 | ||
137 | //___________________________________________________________ | |
138 | void AliMatrixSparse::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n) | |
139 | { | |
140 | // for sym. matrix count how many elems to add have row>=col and assign excplicitly | |
141 | // those which have row<col | |
142 | // | |
143 | // range in increasing order of indices | |
144 | for (int i=n;i--;) { | |
145 | for (int j=i;j>=0;j--) { | |
146 | if (indc[j]>indc[i]) { // swap | |
147 | int ti = indc[i]; indc[i] = indc[j]; indc[j] = ti; | |
148 | double tv = valc[i]; valc[i] = valc[j]; valc[j] = tv; | |
149 | } | |
150 | } | |
151 | } | |
152 | // | |
153 | int ni=n; | |
154 | if (IsSymmetric()) { | |
155 | while(ni--) { | |
156 | if (indc[ni]>r) (*this)(indc[ni],r) += valc[ni]; | |
157 | else break; // use the fact that the indices are ranged in increasing order | |
158 | } | |
159 | } | |
160 | // | |
161 | if (ni<0) return; | |
162 | AliVectorSparse* row = GetRowAdd(r); | |
163 | row->Add(valc,indc,ni+1); | |
164 | } | |
165 | ||
166 | //___________________________________________________________ | |
167 | Float_t AliMatrixSparse::GetDensity() const | |
168 | { | |
169 | // get fraction of non-zero elements | |
170 | Int_t nel = 0; | |
171 | for (int i=GetSize();i--;) nel += GetRow(i)->GetNElems(); | |
172 | int den = IsSymmetric() ? (GetSize()+1)*GetSize()/2 : GetSize()*GetSize(); | |
173 | return float(nel)/den; | |
174 | } | |
175 |