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