STEER/AliMatrixSparse.cxx

   1 #include "AliMatrixSparse.h"
   2 #include "TStopwatch.h"
   3
   4 /**********************************************************************************************/
   5 /* Sparse matrix class, used as a global matrix for AliMillePede2                             */
   6 /*                                                                                            */
   7 /* Author: ruben.shahoyan@cern.ch                                                             */
   8 /*                                                                                            */
   9 /**********************************************************************************************/
  10
  11 //___________________________________________________________
  12 ClassImp(AliMatrixSparse)
  13
  14 //___________________________________________________________
  15 AliMatrixSparse::AliMatrixSparse(Int_t sz)
  16 : AliMatrixSq(),fVecs(0)
  17 {
  18   fNcols=fNrows=sz;
  19   //
  20   fVecs = new AliVectorSparse*[sz];
  21   for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse();
  22 }
  23
  24 //___________________________________________________________
  25 AliMatrixSparse::AliMatrixSparse(const AliMatrixSparse& src)
  26   : AliMatrixSq(src),fVecs(0)
  27 {
  28   fVecs = new AliVectorSparse*[fNcols];
  29   for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
  30 }
  31
  32 //___________________________________________________________
  33 AliVectorSparse* AliMatrixSparse::GetRowAdd(Int_t ir)
  34 {
  35   if (ir>=fNcols) {
  36     AliVectorSparse** arrv = new AliVectorSparse*[ir+1];
  37     for (int i=GetSize();i--;) arrv[i] = fVecs[i];
  38     delete fVecs;
  39     fVecs = arrv;
  40     for (int i=GetSize();i<=ir;i++) fVecs[i] = new AliVectorSparse();
  41     fNcols = fNrows = ir+1;
  42   }
  43   return fVecs[ir];
  44 }
  45
  46 //___________________________________________________________
  47 AliMatrixSparse& AliMatrixSparse::operator=(const AliMatrixSparse& src)
  48 {
  49   if (*this == src) return *this;
  50   Clear();
  51   fNcols = fNrows = src.GetSize();
  52   SetSymmetric(src.IsSymmetric());
  53   fVecs = new AliVectorSparse*[fNcols];
  54   for (int i=fNcols;i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
  55   return *this;
  56 }
  57
  58 //___________________________________________________________
  59 void AliMatrixSparse::Clear(Option_t*)
  60 {
  61   for (int i=fNcols;i--;) delete GetRow(i);
  62   delete[] fVecs;
  63   fNcols = fNrows = 0;
  64 }
  65
  66 //___________________________________________________________
  67 void AliMatrixSparse::Print(Option_t*)  const
  68 {
  69   printf("Sparse Matrix of size %d %s\n",fNcols,IsSymmetric() ? " (Symmetric)":"");
  70   for (int i=0;i<fNcols;i++) {
  71     AliVectorSparse* row = GetRow(i);
  72     if (!row->GetNElems()) continue;
  73     printf("%3d: ",i);
  74     row->Print();
  75   }
  76 }
  77
  78 //___________________________________________________________
  79 void AliMatrixSparse::MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const
  80 {
  81   // fill vecOut by matrix*vecIn
  82   // vector should be of the same size as the matrix
  83   //
  84   memset(vecOut,0,GetSize()*sizeof(Double_t));
  85   //
  86   for (int rw=GetSize();rw--;) {  // loop over rows >>>
  87     const AliVectorSparse* rowV = GetRow(rw);
  88     Int_t nel  = rowV->GetNElems();
  89     if (!nel) continue;
  90     //
  91     UShort_t* indV = rowV->GetIndices();
  92     Double_t* elmV = rowV->GetElems();
  93     //
  94     if (IsSymmetric()) {
  95       // treat diagonal term separately. If filled, it should be the last one
  96       if (indV[--nel]==rw) vecOut[rw] += vecIn[rw]*elmV[nel];
  97       else nel = rowV->GetNElems(); // diag elem was not filled
  98       //
  99       for (int iel=nel;iel--;) {          // less element retrieval for symmetric case
 100         if (elmV[iel]) {
 101           vecOut[rw]        += vecIn[indV[iel]]*elmV[iel];
 102           vecOut[indV[iel]] += vecIn[rw]*elmV[iel];
 103         }
 104       }
 105     }
 106     else for (int iel=nel;iel--;) if (elmV[iel]) vecOut[rw] += vecIn[indV[iel]]*elmV[iel];
 107     //
 108   } // loop over rows <<<
 109   //
 110 }
 111
 112 //___________________________________________________________
 113 void AliMatrixSparse::SortIndices(Bool_t valuesToo)
 114 {
 115   TStopwatch sw;
 116   sw.Start();
 117   printf("AliMatrixSparse:sort>>\n");
 118   // sort columns in increasing order. Used to fix the matrix after ILUk decompostion
 119   for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo);
 120   sw.Stop();
 121   sw.Print();
 122   printf("AliMatrixSparse:sort<<\n");
 123 }
 124
 125 //___________________________________________________________
 126 void AliMatrixSparse::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n)
 127 {
 128   // for sym. matrix count how many elems to add have row>=col and assign excplicitly
 129   // those which have row<col
 130   //
 131   // range in increasing order of indices
 132   for (int i=n;i--;) {
 133     for (int j=i;j>=0;j--) {
 134       if (indc[j]>indc[i]) { // swap
 135         int    ti = indc[i]; indc[i] = indc[j]; indc[j] = ti;
 136         double tv = valc[i]; valc[i] = valc[j]; valc[j] = tv;
 137       }
 138     }
 139   }
 140   //
 141   int ni=n;
 142   if (IsSymmetric()) {
 143     while(ni--) {
 144       if (indc[ni]>r) (*this)(indc[ni],r) += valc[ni];
 145       else break;  // use the fact that the indices are ranged in increasing order
 146     }
 147   }
 148   //
 149   if (ni<0) return;
 150   AliVectorSparse* row = GetRowAdd(r);
 151   row->Add(valc,indc,ni+1);
 152 }
 153
 154 //___________________________________________________________
 155 Float_t AliMatrixSparse::GetDensity() const
 156 {
 157   // get fraction of non-zero elements
 158   Int_t nel = 0;
 159   for (int i=GetSize();i--;) nel += GetRow(i)->GetNElems();
 160   int den = IsSymmetric() ? (GetSize()+1)*GetSize()/2 : GetSize()*GetSize();
 161   return float(nel)/den;
 162 }
 163