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
  21 /**********************************************************************************************/
  22 /* Sparse vector class, used as row of the AliMatrixSparse class                              */
  23 /*                                                                                            */
  24 /* Author: ruben.shahoyan@cern.ch                                                             */
  25 /*                                                                                            */
  26 /**********************************************************************************************/
  27
  28 ClassImp(AliVectorSparse)
  29
  30 //___________________________________________________________
  31 AliVectorSparse::AliVectorSparse()
  32   : fNElems(0),fIndex(0),fElems(0) {}
  33
  34 //___________________________________________________________
  35 AliVectorSparse::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 //___________________________________________________________
  45 void AliVectorSparse::Clear(Option_t*)
  46 {
  47   delete[] fIndex; fIndex = 0;
  48   delete[] fElems; fElems = 0;
  49   fNElems = 0;
  50 }
  51
  52 //___________________________________________________________
  53 AliVectorSparse& 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 //___________________________________________________________
  68 Double_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 //___________________________________________________________
  84 void 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 //___________________________________________________________
  98 Double_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 //__________________________________________________________
 131 void 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 //__________________________________________________________
 154 void 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 //__________________________________________________________
 164 void 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 //___________________________________________________________
 172 void 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