Additional changes for report #70680 AliROOT Coverity DELETE_ARRAY checker fix
[u/mrichter/AliRoot.git] / 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