Fix for raw ctp decoding (Marek)
[u/mrichter/AliRoot.git] / STEER / STEER / AliVectorSparse.cxx
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>
19 #include "AliVectorSparse.h"
20 #include <TString.h>
21
22 /**********************************************************************************************/
23 /* Sparse vector class, used as row of the AliMatrixSparse class                              */
24 /*                                                                                            */ 
25 /* Author: ruben.shahoyan@cern.ch                                                             */
26 /*                                                                                            */ 
27 /**********************************************************************************************/
28
29 ClassImp(AliVectorSparse)
30
31 //___________________________________________________________
32 AliVectorSparse::AliVectorSparse()
33   : fNElems(0),fIndex(0),fElems(0) {}
34
35 //___________________________________________________________
36 AliVectorSparse::AliVectorSparse(const AliVectorSparse& src)
37   : TObject(src),fNElems(src.fNElems),fIndex(0),fElems(0)
38 {
39   // copy c-tor
40   fIndex = new UShort_t[fNElems];
41   fElems = new Double_t[fNElems];
42   memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
43   memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
44 }
45
46 //___________________________________________________________
47 void AliVectorSparse::Clear(Option_t*)
48 {
49   // clear all
50   delete[] fIndex; fIndex = 0;
51   delete[] fElems; fElems = 0;
52   fNElems = 0;
53 }
54
55 //___________________________________________________________
56 AliVectorSparse& AliVectorSparse::operator=(const AliVectorSparse& src)
57 {
58   // assignment op-tor
59   if (&src==this) return *this;
60   Clear();
61   TObject::operator=(src);
62   fNElems = src.fNElems;
63   fIndex = new UShort_t[fNElems];
64   fElems = new Double_t[fNElems];
65   memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
66   memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
67   //
68   return *this;
69 }
70
71 //___________________________________________________________
72 Double_t AliVectorSparse::FindIndex(Int_t ind) const
73 {
74   // return an element with given index
75   //printf("V: findindex\n");
76   int first = 0;
77   int last = fNElems-1;
78   while (first<=last) {
79     int mid = (first+last)>>1;
80     if (ind>fIndex[mid]) first = mid+1;
81     else if (ind<fIndex[mid]) last = mid-1;
82     else return fElems[mid];
83   }
84   return 0.0;
85 }
86
87 //___________________________________________________________
88 void AliVectorSparse::SetToZero(Int_t ind)
89 {
90   // set element to 0 if it was already defined
91   int first = 0;
92   int last = fNElems-1;
93   while (first<=last) {
94     int mid = (first+last)>>1;
95     if (ind>fIndex[mid]) first = mid+1;
96     else if (ind<fIndex[mid]) last = mid-1;
97     else {fElems[mid] = 0.; return;}
98   }
99 }
100
101 //___________________________________________________________
102 Double_t& AliVectorSparse::FindIndexAdd(Int_t ind)
103 {
104   // increment an element with given index
105   //printf("V: findindexAdd\n");
106   int first = 0;
107   int last = fNElems-1;
108   while (first<=last) {
109     int mid = (first+last)>>1;
110     if (ind>fIndex[mid]) first = mid+1;
111     else if (ind<fIndex[mid]) last = mid-1;
112     else return fElems[mid];
113   }
114   // need to insert a new element
115   UShort_t *arrI = new UShort_t[fNElems+1];
116   memcpy(arrI,fIndex,first*sizeof(UShort_t));
117   arrI[first] = ind;
118   memcpy(arrI+first+1,fIndex+first,(fNElems-first)*sizeof(UShort_t));
119   delete[] fIndex;
120   fIndex = arrI;
121   //
122   Double_t   *arrE = new Double_t[fNElems+1];
123   memcpy(arrE,fElems,first*sizeof(Double_t));
124   arrE[first] = 0;
125   memcpy(arrE+first+1,fElems+first,(fNElems-first)*sizeof(Double_t));
126   delete[] fElems;
127   fElems = arrE;
128   //
129   fNElems++;
130   return fElems[first];
131   //
132 }
133
134 //__________________________________________________________
135 void AliVectorSparse::ReSize(Int_t sz,Bool_t copy)
136 {
137   // change the size
138   if (sz<1) {Clear(); return;}
139     // need to insert a new element
140   UShort_t *arrI = new UShort_t[sz];
141   Double_t *arrE = new Double_t[sz];
142   memset(arrI,0,sz*sizeof(UShort_t));
143   memset(arrE,0,sz*sizeof(Double_t));
144   //
145   if (copy && fIndex) {
146     int cpsz = TMath::Min(fNElems,sz);
147     memcpy(arrI,fIndex,cpsz*sizeof(UShort_t));
148     memcpy(arrE,fElems,cpsz*sizeof(Double_t));
149   }
150   delete[] fIndex;
151   delete[] fElems;
152   fIndex = arrI;
153   fElems = arrE;
154   fNElems = sz;
155   //
156 }
157
158 //__________________________________________________________
159 void AliVectorSparse::SortIndices(Bool_t valuesToo)
160 {
161   // sort indices in increasing order. Used to fix the row after ILUk decomposition
162   for (int i=fNElems;i--;) for (int j=i;j--;) if (fIndex[i]<fIndex[j]) { //swap
163         UShort_t tmpI = fIndex[i]; fIndex[i] = fIndex[j]; fIndex[j]=tmpI;
164         if (valuesToo) {Double_t tmpV = fElems[i];fElems[i]=fElems[j];fElems[j]=tmpV;}
165       }
166 }
167
168 //__________________________________________________________
169 void AliVectorSparse::Print(Option_t* opt)  const
170 {
171   // print itself
172   TString sopt = opt; sopt.ToLower();
173   int ndig = sopt.Atoi();
174   if (ndig<=1) ndig = 2;
175   sopt = "%2d:%+.";
176   sopt += ndig;
177   sopt += "e |";
178   printf("|");
179   for (int i=0;i<fNElems;i++) printf(sopt.Data(),fIndex[i],fElems[i]);
180   printf("\n");
181 }
182
183 //___________________________________________________________
184 void AliVectorSparse::Add(Double_t *valc,Int_t *indc,Int_t n)
185 {
186   // add indiced array to row. Indices must be in increasing order
187   int indx;
188   int nadd = 0;
189   //
190   int last = fNElems-1;
191   int mid = 0;
192   for (int i=n;i--;) {
193     // if the element with this index is already defined, just add the value
194     int first = 0;
195     Bool_t toAdd = kTRUE;
196     indx = indc[i];
197     while (first<=last) {
198       mid = (first+last)>>1;
199       if (indx>fIndex[mid]) first = mid+1;
200       else if (indx<fIndex[mid]) last = mid-1;
201       else {
202         fElems[mid] += valc[i];
203         indc[i] = -1;
204         toAdd = kFALSE;
205         last = mid-1;   // profit from the indices being ordered
206         break;
207       }
208     }
209     if (toAdd) nadd++;
210   }
211   //
212   if (nadd<1) return; // nothing to do anymore
213   // 
214   // need to expand the row
215   UShort_t *arrI = new UShort_t[fNElems+nadd];
216   Double_t *arrE = new Double_t[fNElems+nadd];
217   // copy old elems embedding the new ones
218   int inew=0,iold=0;
219   for (int i=0;i<n;i++) {  
220     if ( (indx=indc[i])<0) continue;
221     while (iold<fNElems && fIndex[iold]<indx) {
222       arrI[inew]   = fIndex[iold];
223       arrE[inew++] = fElems[iold++];
224     }
225     arrI[inew] = indx;
226     arrE[inew++] = valc[i];
227   }
228   // copy the rest
229   while (iold<fNElems) {
230     arrI[inew]   = fIndex[iold];
231     arrE[inew++] = fElems[iold++];
232   }
233   //
234   delete[] fIndex;
235   delete[] fElems;
236   fIndex = arrI;
237   fElems = arrE;
238   //
239   fNElems += nadd;
240   //
241 }
242