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