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