Fixes for coverity.
[u/mrichter/AliRoot.git] / 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