Changed Zero to SetToZero to avoid hiding of the base class method
[u/mrichter/AliRoot.git] / STEER / AliMatrixSparse.cxx
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 //___________________________________________________________
66 void AliVectorSparse::SetToZero(Int_t ind)
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 }