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