STEER/AliMatrixSparse.cxx

   1 #include "AliMatrixSparse.h"
   2
   3 //___________________________________________________________
   4 AliVectorSparse::AliVectorSparse()
   5   : fNElems(0),fIndex(0),fElems(0) {}
   6
   7 //___________________________________________________________
   8 AliVectorSparse::AliVectorSparse(const AliVectorSparse& src)
   9   : fNElems(src.fNElems),fIndex(0),fElems(0)
  10 {
  11   fIndex = new UShort_t[fNElems];
  12   fElems = new Double_t[fNElems];
  13   memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
  14   memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
  15 }
  16
  17 //___________________________________________________________
  18 void AliVectorSparse::Clear()
  19 {
  20   delete[] fIndex; fIndex = 0;
  21   delete[] fElems; fElems = 0;
  22   fNElems = 0;
  23 }
  24
  25 //___________________________________________________________
  26 Float_t AliMatrixSparse::GetDensity() const
  27 {
  28   // get fraction of non-zero elements
  29   Int_t nel = 0;
  30   for (int i=GetSize();i--;) nel += GetRow(i)->GetNElems();
  31   int den = IsSymmetric() ? (GetSize()+1)*GetSize()/2 : GetSize()*GetSize();
  32   return float(nel)/den;
  33 }
  34
  35 //___________________________________________________________
  36 AliVectorSparse& AliVectorSparse::operator=(const AliVectorSparse& src)
  37 {
  38   if (&src==this) return *this;
  39   Clear();
  40   fNElems = src.fNElems;
  41   fIndex = new UShort_t[fNElems];
  42   fElems = new Double_t[fNElems];
  43   memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
  44   memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
  45   //
  46   return *this;
  47 }
  48
  49 //___________________________________________________________
  50 Double_t AliVectorSparse::FindIndex(Int_t ind) const
  51 {
  52   // return an element with given index
  53   //printf("V: findindex\n");
  54   int first = 0;
  55   int last = fNElems-1;
  56   while (first<=last) {
  57     int mid = (first+last)>>1;
  58     if (ind>fIndex[mid]) first = mid+1;
  59     else if (ind<fIndex[mid]) last = mid-1;
  60     else return fElems[mid];
  61   }
  62   return 0.0;
  63 }
  64
  65 //___________________________________________________________
  66 void AliVectorSparse::SetToZero(Int_t ind)
  67 {
  68   // set element to 0 if it was already defined
  69   int first = 0;
  70   int last = fNElems-1;
  71   while (first<=last) {
  72     int mid = (first+last)>>1;
  73     if (ind>fIndex[mid]) first = mid+1;
  74     else if (ind<fIndex[mid]) last = mid-1;
  75     else {fElems[mid] = 0.; return;}
  76   }
  77 }
  78
  79 //___________________________________________________________
  80 Double_t& AliVectorSparse::FindIndexAdd(Int_t ind)
  81 {
  82   // increment an element with given index
  83   //printf("V: findindexAdd\n");
  84   int first = 0;
  85   int last = fNElems-1;
  86   while (first<=last) {
  87     int mid = (first+last)>>1;
  88     if (ind>fIndex[mid]) first = mid+1;
  89     else if (ind<fIndex[mid]) last = mid-1;
  90     else return fElems[mid];
  91   }
  92   // need to insert a new element
  93   UShort_t *arrI = new UShort_t[fNElems+1];
  94   memcpy(arrI,fIndex,first*sizeof(UShort_t));
  95   arrI[first] = ind;
  96   memcpy(arrI+first+1,fIndex+first,(fNElems-first)*sizeof(UShort_t));
  97   delete[] fIndex;
  98   fIndex = arrI;
  99   //
 100   Double_t   *arrE = new Double_t[fNElems+1];
 101   memcpy(arrE,fElems,first*sizeof(Double_t));
 102   arrE[first] = 0;
 103   memcpy(arrE+first+1,fElems+first,(fNElems-first)*sizeof(Double_t));
 104   delete[] fElems;
 105   fElems = arrE;
 106   //
 107   fNElems++;
 108   return fElems[first];
 109   //
 110 }
 111
 112 //__________________________________________________________
 113 void AliVectorSparse::ReSize(Int_t sz,Bool_t copy)
 114 {
 115   if (sz<1) {Clear(); return;}
 116     // need to insert a new element
 117   UShort_t *arrI = new UShort_t[sz];
 118   Double_t *arrE = new Double_t[sz];
 119   memset(arrI,0,sz*sizeof(UShort_t));
 120   memset(arrE,0,sz*sizeof(Double_t));
 121   //
 122   if (copy && fIndex) {
 123     int cpsz = TMath::Min(fNElems,sz);
 124     memcpy(arrI,fIndex,cpsz*sizeof(UShort_t));
 125     memcpy(arrE,fElems,cpsz*sizeof(Double_t));
 126   }
 127   delete[] fIndex;
 128   delete[] fElems;
 129   fIndex = arrI;
 130   fElems = arrE;
 131   fNElems = sz;
 132   //
 133 }
 134
 135 //__________________________________________________________
 136 void AliVectorSparse::SortIndices(Bool_t valuesToo)
 137 {
 138   // sort indices in increasing order. Used to fix the row after ILUk decomposition
 139   for (int i=fNElems;i--;) for (int j=i;j--;) if (fIndex[i]<fIndex[j]) { //swap
 140         UShort_t tmpI = fIndex[i]; fIndex[i] = fIndex[j]; fIndex[j]=tmpI;
 141         if (valuesToo) {Double_t tmpV = fElems[i];fElems[i]=fElems[j];fElems[j]=tmpV;}
 142       }
 143 }
 144
 145 //__________________________________________________________
 146 void AliVectorSparse::Print(Option_t* )  const
 147 {
 148   printf("|");
 149   for (int i=0;i<fNElems;i++) printf("%2d:%+.2e|",fIndex[i],fElems[i]);
 150   printf("|\n");
 151 }
 152
 153
 154 ///////////////////////////////////////////////////////////////////////////////////////////
 155 //___________________________________________________________
 156 ClassImp(AliMatrixSparse)
 157
 158 //___________________________________________________________
 159 AliMatrixSparse::AliMatrixSparse(Int_t sz)
 160 : AliMatrixSq(),fVecs(0)
 161 {
 162   fNcols=fNrows=sz;
 163   //
 164   fVecs = new AliVectorSparse*[sz];
 165   for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse();
 166 }
 167
 168 //___________________________________________________________
 169 AliMatrixSparse::AliMatrixSparse(const AliMatrixSparse& src)
 170   : AliMatrixSq(src),fVecs(0)
 171 {
 172   fVecs = new AliVectorSparse*[fNcols];
 173   for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
 174 }
 175
 176 //___________________________________________________________
 177 AliVectorSparse* AliMatrixSparse::GetRowAdd(Int_t ir)
 178 {
 179   if (ir>=fNcols) {
 180     AliVectorSparse** arrv = new AliVectorSparse*[ir+1];
 181     for (int i=GetSize();i--;) arrv[i] = fVecs[i];
 182     delete fVecs;
 183     fVecs = arrv;
 184     for (int i=GetSize();i<=ir;i++) fVecs[i] = new AliVectorSparse();
 185     fNcols = fNrows = ir+1;
 186   }
 187   return fVecs[ir];
 188 }
 189
 190 //___________________________________________________________
 191 AliMatrixSparse& AliMatrixSparse::operator=(const AliMatrixSparse& src)
 192 {
 193   if (*this == src) return *this;
 194   Clear();
 195   fNcols = fNrows = src.GetSize();
 196   SetSymmetric(src.IsSymmetric());
 197   fVecs = new AliVectorSparse*[fNcols];
 198   for (int i=fNcols;i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
 199   return *this;
 200 }
 201
 202 //___________________________________________________________
 203 void AliMatrixSparse::Clear(Option_t*)
 204 {
 205   for (int i=fNcols;i--;) delete GetRow(i);
 206   delete[] fVecs;
 207   fNcols = fNrows = 0;
 208 }
 209
 210 //___________________________________________________________
 211 void AliMatrixSparse::Print(Option_t*)  const
 212 {
 213   printf("Sparse Matrix of size %d %s\n",fNcols,IsSymmetric() ? " (Symmetric)":"");
 214   for (int i=0;i<fNcols;i++) {
 215     AliVectorSparse* row = GetRow(i);
 216     if (!row->GetNElems()) continue;
 217     printf("%3d: ",i);
 218     row->Print();
 219   }
 220 }
 221
 222 //___________________________________________________________
 223 void AliMatrixSparse::MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const
 224 {
 225   // fill vecOut by matrix*vecIn
 226   // vector should be of the same size as the matrix
 227   //
 228   memset(vecOut,0,GetSize()*sizeof(Double_t));
 229   //
 230   for (int rw=GetSize();rw--;) {  // loop over rows >>>
 231     const AliVectorSparse* rowV = GetRow(rw);
 232     Int_t nel  = rowV->GetNElems();
 233     if (!nel) continue;
 234     //
 235     UShort_t* indV = rowV->GetIndices();
 236     Double_t* elmV = rowV->GetElems();
 237     //
 238     if (IsSymmetric()) {
 239       // treat diagonal term separately. If filled, it should be the last one
 240       if (indV[--nel]==rw) vecOut[rw] += vecIn[rw]*elmV[nel];
 241       else nel = rowV->GetNElems(); // diag elem was not filled
 242       //
 243       for (int iel=nel;iel--;) {          // less element retrieval for symmetric case
 244         if (elmV[iel]) {
 245           vecOut[rw]        += vecIn[indV[iel]]*elmV[iel];
 246           vecOut[indV[iel]] += vecIn[rw]*elmV[iel];
 247         }
 248       }
 249     }
 250     else for (int iel=nel;iel--;) if (elmV[iel]) vecOut[rw] += vecIn[indV[iel]]*elmV[iel];
 251     //
 252   } // loop over rows <<<
 253   //
 254 }
 255
 256 //___________________________________________________________
 257 void AliMatrixSparse::SortIndices(Bool_t valuesToo)
 258 {
 259   // sort columns in increasing order. Used to fix the matrix after ILUk decompostion
 260   for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo);
 261 }