Fixes for bug #52499: Field polarities inconsistiency
[u/mrichter/AliRoot.git] / STEER / AliVectorSparse.cxx
CommitLineData
81b4a4ef 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id: AliRun.cxx 30859 2009-02-02 16:24:37Z fca $ */
17
18#include <string.h>
de34b538 19#include "AliVectorSparse.h"
20
21/**********************************************************************************************/
22/* Sparse vector class, used as row of the AliMatrixSparse class */
23/* */
24/* Author: ruben.shahoyan@cern.ch */
25/* */
26/**********************************************************************************************/
27
28ClassImp(AliVectorSparse)
29
30//___________________________________________________________
31AliVectorSparse::AliVectorSparse()
32 : fNElems(0),fIndex(0),fElems(0) {}
33
34//___________________________________________________________
35AliVectorSparse::AliVectorSparse(const AliVectorSparse& src)
36 : TObject(src),fNElems(src.fNElems),fIndex(0),fElems(0)
37{
38 fIndex = new UShort_t[fNElems];
39 fElems = new Double_t[fNElems];
40 memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
41 memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
42}
43
44//___________________________________________________________
45void AliVectorSparse::Clear(Option_t*)
46{
47 delete[] fIndex; fIndex = 0;
48 delete[] fElems; fElems = 0;
49 fNElems = 0;
50}
51
52//___________________________________________________________
53AliVectorSparse& AliVectorSparse::operator=(const AliVectorSparse& src)
54{
55 if (&src==this) return *this;
56 Clear();
57 TObject::operator=(src);
58 fNElems = src.fNElems;
59 fIndex = new UShort_t[fNElems];
60 fElems = new Double_t[fNElems];
61 memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
62 memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
63 //
64 return *this;
65}
66
67//___________________________________________________________
68Double_t AliVectorSparse::FindIndex(Int_t ind) const
69{
70 // return an element with given index
71 //printf("V: findindex\n");
72 int first = 0;
73 int last = fNElems-1;
74 while (first<=last) {
75 int mid = (first+last)>>1;
76 if (ind>fIndex[mid]) first = mid+1;
77 else if (ind<fIndex[mid]) last = mid-1;
78 else return fElems[mid];
79 }
80 return 0.0;
81}
82
83//___________________________________________________________
84void AliVectorSparse::SetToZero(Int_t ind)
85{
86 // set element to 0 if it was already defined
87 int first = 0;
88 int last = fNElems-1;
89 while (first<=last) {
90 int mid = (first+last)>>1;
91 if (ind>fIndex[mid]) first = mid+1;
92 else if (ind<fIndex[mid]) last = mid-1;
93 else {fElems[mid] = 0.; return;}
94 }
95}
96
97//___________________________________________________________
98Double_t& AliVectorSparse::FindIndexAdd(Int_t ind)
99{
100 // increment an element with given index
101 //printf("V: findindexAdd\n");
102 int first = 0;
103 int last = fNElems-1;
104 while (first<=last) {
105 int mid = (first+last)>>1;
106 if (ind>fIndex[mid]) first = mid+1;
107 else if (ind<fIndex[mid]) last = mid-1;
108 else return fElems[mid];
109 }
110 // need to insert a new element
111 UShort_t *arrI = new UShort_t[fNElems+1];
112 memcpy(arrI,fIndex,first*sizeof(UShort_t));
113 arrI[first] = ind;
114 memcpy(arrI+first+1,fIndex+first,(fNElems-first)*sizeof(UShort_t));
115 delete[] fIndex;
116 fIndex = arrI;
117 //
118 Double_t *arrE = new Double_t[fNElems+1];
119 memcpy(arrE,fElems,first*sizeof(Double_t));
120 arrE[first] = 0;
121 memcpy(arrE+first+1,fElems+first,(fNElems-first)*sizeof(Double_t));
122 delete[] fElems;
123 fElems = arrE;
124 //
125 fNElems++;
126 return fElems[first];
127 //
128}
129
130//__________________________________________________________
131void AliVectorSparse::ReSize(Int_t sz,Bool_t copy)
132{
133 if (sz<1) {Clear(); return;}
134 // need to insert a new element
135 UShort_t *arrI = new UShort_t[sz];
136 Double_t *arrE = new Double_t[sz];
137 memset(arrI,0,sz*sizeof(UShort_t));
138 memset(arrE,0,sz*sizeof(Double_t));
139 //
140 if (copy && fIndex) {
141 int cpsz = TMath::Min(fNElems,sz);
142 memcpy(arrI,fIndex,cpsz*sizeof(UShort_t));
143 memcpy(arrE,fElems,cpsz*sizeof(Double_t));
144 }
145 delete[] fIndex;
146 delete[] fElems;
147 fIndex = arrI;
148 fElems = arrE;
149 fNElems = sz;
150 //
151}
152
153//__________________________________________________________
154void AliVectorSparse::SortIndices(Bool_t valuesToo)
155{
156 // sort indices in increasing order. Used to fix the row after ILUk decomposition
157 for (int i=fNElems;i--;) for (int j=i;j--;) if (fIndex[i]<fIndex[j]) { //swap
158 UShort_t tmpI = fIndex[i]; fIndex[i] = fIndex[j]; fIndex[j]=tmpI;
159 if (valuesToo) {Double_t tmpV = fElems[i];fElems[i]=fElems[j];fElems[j]=tmpV;}
160 }
161}
162
163//__________________________________________________________
164void AliVectorSparse::Print(Option_t* ) const
165{
166 printf("|");
167 for (int i=0;i<fNElems;i++) printf("%2d:%+.2e|",fIndex[i],fElems[i]);
168 printf("|\n");
169}
170
171//___________________________________________________________
172void AliVectorSparse::Add(Double_t *valc,Int_t *indc,Int_t n)
173{
174 // add indiced array to row. Indices must be in increasing order
175 int indx;
176 int nadd = 0;
177 //
178 int last = fNElems-1;
179 int mid = 0;
180 for (int i=n;i--;) {
181 // if the element with this index is already defined, just add the value
182 int first = 0;
183 Bool_t toAdd = kTRUE;
184 indx = indc[i];
185 while (first<=last) {
186 mid = (first+last)>>1;
187 if (indx>fIndex[mid]) first = mid+1;
188 else if (indx<fIndex[mid]) last = mid-1;
189 else {
190 fElems[mid] += valc[i];
191 indc[i] = -1;
192 toAdd = kFALSE;
193 last = mid-1; // profit from the indices being ordered
194 break;
195 }
196 }
197 if (toAdd) nadd++;
198 }
199 //
200 if (nadd<1) return; // nothing to do anymore
201 //
202 // need to expand the row
203 UShort_t *arrI = new UShort_t[fNElems+nadd];
204 Double_t *arrE = new Double_t[fNElems+nadd];
205 // copy old elems embedding the new ones
206 int inew=0,iold=0;
207 for (int i=0;i<n;i++) {
208 if ( (indx=indc[i])<0) continue;
209 while (iold<fNElems && fIndex[iold]<indx) {
210 arrI[inew] = fIndex[iold];
211 arrE[inew++] = fElems[iold++];
212 }
213 arrI[inew] = indx;
214 arrE[inew++] = valc[i];
215 }
216 // copy the rest
217 while (iold<fNElems) {
218 arrI[inew] = fIndex[iold];
219 arrE[inew++] = fElems[iold++];
220 }
221 //
222 delete[] fIndex;
223 delete[] fElems;
224 fIndex = arrI;
225 fElems = arrE;
226 //
227 fNElems += nadd;
228 //
229}
230