]>
Commit | Line | Data |
---|---|---|
8a9ab0eb | 1 | #include "AliMatrixSparse.h" |
2 | ||
3 | //___________________________________________________________ | |
4 | AliVectorSparse::AliVectorSparse() | |
5 | : fNElems(0),fIndex(0),fElems(0) {} | |
6 | ||
7 | //___________________________________________________________ | |
8 | AliVectorSparse::AliVectorSparse(const AliVectorSparse& src) | |
9 | : fNElems(src.fNElems),fIndex(0),fElems(0) | |
10 | { | |
11 | fIndex = new UShort_t[fNElems]; | |
12 | fElems = new Double_t[fNElems]; | |
13 | memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t)); | |
14 | memcpy(fElems,src.fElems,fNElems*sizeof(Double_t)); | |
15 | } | |
16 | ||
17 | //___________________________________________________________ | |
18 | void AliVectorSparse::Clear() | |
19 | { | |
20 | delete[] fIndex; fIndex = 0; | |
21 | delete[] fElems; fElems = 0; | |
22 | fNElems = 0; | |
23 | } | |
24 | ||
25 | //___________________________________________________________ | |
26 | Float_t AliMatrixSparse::GetDensity() const | |
27 | { | |
28 | // get fraction of non-zero elements | |
29 | Int_t nel = 0; | |
30 | for (int i=GetSize();i--;) nel += GetRow(i)->GetNElems(); | |
31 | int den = IsSymmetric() ? (GetSize()+1)*GetSize()/2 : GetSize()*GetSize(); | |
32 | return float(nel)/den; | |
33 | } | |
34 | ||
35 | //___________________________________________________________ | |
36 | AliVectorSparse& AliVectorSparse::operator=(const AliVectorSparse& src) | |
37 | { | |
38 | if (&src==this) return *this; | |
39 | Clear(); | |
40 | fNElems = src.fNElems; | |
41 | fIndex = new UShort_t[fNElems]; | |
42 | fElems = new Double_t[fNElems]; | |
43 | memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t)); | |
44 | memcpy(fElems,src.fElems,fNElems*sizeof(Double_t)); | |
45 | // | |
46 | return *this; | |
47 | } | |
48 | ||
49 | //___________________________________________________________ | |
50 | Double_t AliVectorSparse::FindIndex(Int_t ind) const | |
51 | { | |
52 | // return an element with given index | |
53 | //printf("V: findindex\n"); | |
54 | int first = 0; | |
55 | int last = fNElems-1; | |
56 | while (first<=last) { | |
57 | int mid = (first+last)>>1; | |
58 | if (ind>fIndex[mid]) first = mid+1; | |
59 | else if (ind<fIndex[mid]) last = mid-1; | |
60 | else return fElems[mid]; | |
61 | } | |
62 | return 0.0; | |
63 | } | |
64 | ||
65 | //___________________________________________________________ | |
f748efd7 | 66 | void AliVectorSparse::SetToZero(Int_t ind) |
8a9ab0eb | 67 | { |
68 | // set element to 0 if it was already defined | |
69 | int first = 0; | |
70 | int last = fNElems-1; | |
71 | while (first<=last) { | |
72 | int mid = (first+last)>>1; | |
73 | if (ind>fIndex[mid]) first = mid+1; | |
74 | else if (ind<fIndex[mid]) last = mid-1; | |
75 | else {fElems[mid] = 0.; return;} | |
76 | } | |
77 | } | |
78 | ||
79 | //___________________________________________________________ | |
80 | Double_t& AliVectorSparse::FindIndexAdd(Int_t ind) | |
81 | { | |
82 | // increment an element with given index | |
83 | //printf("V: findindexAdd\n"); | |
84 | int first = 0; | |
85 | int last = fNElems-1; | |
86 | while (first<=last) { | |
87 | int mid = (first+last)>>1; | |
88 | if (ind>fIndex[mid]) first = mid+1; | |
89 | else if (ind<fIndex[mid]) last = mid-1; | |
90 | else return fElems[mid]; | |
91 | } | |
92 | // need to insert a new element | |
93 | UShort_t *arrI = new UShort_t[fNElems+1]; | |
94 | memcpy(arrI,fIndex,first*sizeof(UShort_t)); | |
95 | arrI[first] = ind; | |
96 | memcpy(arrI+first+1,fIndex+first,(fNElems-first)*sizeof(UShort_t)); | |
97 | delete[] fIndex; | |
98 | fIndex = arrI; | |
99 | // | |
100 | Double_t *arrE = new Double_t[fNElems+1]; | |
101 | memcpy(arrE,fElems,first*sizeof(Double_t)); | |
102 | arrE[first] = 0; | |
103 | memcpy(arrE+first+1,fElems+first,(fNElems-first)*sizeof(Double_t)); | |
104 | delete[] fElems; | |
105 | fElems = arrE; | |
106 | // | |
107 | fNElems++; | |
108 | return fElems[first]; | |
109 | // | |
110 | } | |
111 | ||
112 | //__________________________________________________________ | |
113 | void AliVectorSparse::ReSize(Int_t sz,Bool_t copy) | |
114 | { | |
115 | if (sz<1) {Clear(); return;} | |
116 | // need to insert a new element | |
117 | UShort_t *arrI = new UShort_t[sz]; | |
118 | Double_t *arrE = new Double_t[sz]; | |
119 | memset(arrI,0,sz*sizeof(UShort_t)); | |
120 | memset(arrE,0,sz*sizeof(Double_t)); | |
121 | // | |
122 | if (copy && fIndex) { | |
123 | int cpsz = TMath::Min(fNElems,sz); | |
124 | memcpy(arrI,fIndex,cpsz*sizeof(UShort_t)); | |
125 | memcpy(arrE,fElems,cpsz*sizeof(Double_t)); | |
126 | } | |
127 | delete[] fIndex; | |
128 | delete[] fElems; | |
129 | fIndex = arrI; | |
130 | fElems = arrE; | |
131 | fNElems = sz; | |
132 | // | |
133 | } | |
134 | ||
135 | //__________________________________________________________ | |
136 | void AliVectorSparse::SortIndices(Bool_t valuesToo) | |
137 | { | |
138 | // sort indices in increasing order. Used to fix the row after ILUk decomposition | |
139 | for (int i=fNElems;i--;) for (int j=i;j--;) if (fIndex[i]<fIndex[j]) { //swap | |
140 | UShort_t tmpI = fIndex[i]; fIndex[i] = fIndex[j]; fIndex[j]=tmpI; | |
141 | if (valuesToo) {Double_t tmpV = fElems[i];fElems[i]=fElems[j];fElems[j]=tmpV;} | |
142 | } | |
143 | } | |
144 | ||
145 | //__________________________________________________________ | |
146 | void AliVectorSparse::Print(Option_t* ) const | |
147 | { | |
148 | printf("|"); | |
149 | for (int i=0;i<fNElems;i++) printf("%2d:%+.2e|",fIndex[i],fElems[i]); | |
150 | printf("|\n"); | |
151 | } | |
152 | ||
153 | ||
154 | /////////////////////////////////////////////////////////////////////////////////////////// | |
155 | //___________________________________________________________ | |
156 | ClassImp(AliMatrixSparse) | |
157 | ||
158 | //___________________________________________________________ | |
159 | AliMatrixSparse::AliMatrixSparse(Int_t sz) | |
160 | : AliMatrixSq(),fVecs(0) | |
161 | { | |
162 | fNcols=fNrows=sz; | |
163 | // | |
164 | fVecs = new AliVectorSparse*[sz]; | |
165 | for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse(); | |
166 | } | |
167 | ||
168 | //___________________________________________________________ | |
169 | AliMatrixSparse::AliMatrixSparse(const AliMatrixSparse& src) | |
170 | : AliMatrixSq(src),fVecs(0) | |
171 | { | |
172 | fVecs = new AliVectorSparse*[fNcols]; | |
173 | for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i)); | |
174 | } | |
175 | ||
176 | //___________________________________________________________ | |
177 | AliVectorSparse* AliMatrixSparse::GetRowAdd(Int_t ir) | |
178 | { | |
179 | if (ir>=fNcols) { | |
180 | AliVectorSparse** arrv = new AliVectorSparse*[ir+1]; | |
181 | for (int i=GetSize();i--;) arrv[i] = fVecs[i]; | |
182 | delete fVecs; | |
183 | fVecs = arrv; | |
184 | for (int i=GetSize();i<=ir;i++) fVecs[i] = new AliVectorSparse(); | |
185 | fNcols = fNrows = ir+1; | |
186 | } | |
187 | return fVecs[ir]; | |
188 | } | |
189 | ||
190 | //___________________________________________________________ | |
191 | AliMatrixSparse& AliMatrixSparse::operator=(const AliMatrixSparse& src) | |
192 | { | |
193 | if (*this == src) return *this; | |
194 | Clear(); | |
195 | fNcols = fNrows = src.GetSize(); | |
196 | SetSymmetric(src.IsSymmetric()); | |
197 | fVecs = new AliVectorSparse*[fNcols]; | |
198 | for (int i=fNcols;i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i)); | |
199 | return *this; | |
200 | } | |
201 | ||
202 | //___________________________________________________________ | |
203 | void AliMatrixSparse::Clear(Option_t*) | |
204 | { | |
205 | for (int i=fNcols;i--;) delete GetRow(i); | |
206 | delete[] fVecs; | |
207 | fNcols = fNrows = 0; | |
208 | } | |
209 | ||
210 | //___________________________________________________________ | |
211 | void AliMatrixSparse::Print(Option_t*) const | |
212 | { | |
213 | printf("Sparse Matrix of size %d %s\n",fNcols,IsSymmetric() ? " (Symmetric)":""); | |
214 | for (int i=0;i<fNcols;i++) { | |
215 | AliVectorSparse* row = GetRow(i); | |
216 | if (!row->GetNElems()) continue; | |
217 | printf("%3d: ",i); | |
218 | row->Print(); | |
219 | } | |
220 | } | |
221 | ||
222 | //___________________________________________________________ | |
223 | void AliMatrixSparse::MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const | |
224 | { | |
225 | // fill vecOut by matrix*vecIn | |
226 | // vector should be of the same size as the matrix | |
227 | // | |
228 | memset(vecOut,0,GetSize()*sizeof(Double_t)); | |
229 | // | |
230 | for (int rw=GetSize();rw--;) { // loop over rows >>> | |
231 | const AliVectorSparse* rowV = GetRow(rw); | |
232 | Int_t nel = rowV->GetNElems(); | |
233 | if (!nel) continue; | |
234 | // | |
235 | UShort_t* indV = rowV->GetIndices(); | |
236 | Double_t* elmV = rowV->GetElems(); | |
237 | // | |
238 | if (IsSymmetric()) { | |
239 | // treat diagonal term separately. If filled, it should be the last one | |
240 | if (indV[--nel]==rw) vecOut[rw] += vecIn[rw]*elmV[nel]; | |
241 | else nel = rowV->GetNElems(); // diag elem was not filled | |
242 | // | |
243 | for (int iel=nel;iel--;) { // less element retrieval for symmetric case | |
244 | if (elmV[iel]) { | |
245 | vecOut[rw] += vecIn[indV[iel]]*elmV[iel]; | |
246 | vecOut[indV[iel]] += vecIn[rw]*elmV[iel]; | |
247 | } | |
248 | } | |
249 | } | |
250 | else for (int iel=nel;iel--;) if (elmV[iel]) vecOut[rw] += vecIn[indV[iel]]*elmV[iel]; | |
251 | // | |
252 | } // loop over rows <<< | |
253 | // | |
254 | } | |
255 | ||
256 | //___________________________________________________________ | |
257 | void AliMatrixSparse::SortIndices(Bool_t valuesToo) | |
258 | { | |
259 | // sort columns in increasing order. Used to fix the matrix after ILUk decompostion | |
260 | for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo); | |
261 | } |