Coverity 16571
[u/mrichter/AliRoot.git] / STEER / AliMatrixSparse.cxx
CommitLineData
8a9ab0eb 1#include "AliMatrixSparse.h"
de34b538 2#include "TStopwatch.h"
8a9ab0eb 3
de34b538 4/**********************************************************************************************/
5/* Sparse matrix class, used as a global matrix for AliMillePede2 */
6/* */
7/* Author: ruben.shahoyan@cern.ch */
8/* */
9/**********************************************************************************************/
8a9ab0eb 10
8a9ab0eb 11//___________________________________________________________
12ClassImp(AliMatrixSparse)
13
14//___________________________________________________________
15AliMatrixSparse::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//___________________________________________________________
25AliMatrixSparse::AliMatrixSparse(const AliMatrixSparse& src)
26 : AliMatrixSq(src),fVecs(0)
27{
3bc5dd7d 28 fVecs = new AliVectorSparse*[src.GetSize()];
8a9ab0eb 29 for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
30}
31
32//___________________________________________________________
33AliVectorSparse* AliMatrixSparse::GetRowAdd(Int_t ir)
34{
3bc5dd7d 35 if (ir>=fNrows) {
8a9ab0eb 36 AliVectorSparse** arrv = new AliVectorSparse*[ir+1];
37 for (int i=GetSize();i--;) arrv[i] = fVecs[i];
6fefa5ce 38 delete [] fVecs;
8a9ab0eb 39 fVecs = arrv;
40 for (int i=GetSize();i<=ir;i++) fVecs[i] = new AliVectorSparse();
3bc5dd7d 41 fNrows = ir+1;
42 if (IsSymmetric() && fNcols<fNrows) fNcols = fNrows;
8a9ab0eb 43 }
44 return fVecs[ir];
45}
46
47//___________________________________________________________
48AliMatrixSparse& AliMatrixSparse::operator=(const AliMatrixSparse& src)
49{
50 if (*this == src) return *this;
51 Clear();
3bc5dd7d 52 fNcols = src.GetNCols();
53 fNrows = src.GetNRows();
8a9ab0eb 54 SetSymmetric(src.IsSymmetric());
3bc5dd7d 55 fVecs = new AliVectorSparse*[fNrows];
56 for (int i=fNrows;i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
8a9ab0eb 57 return *this;
58}
59
60//___________________________________________________________
61void AliMatrixSparse::Clear(Option_t*)
62{
3bc5dd7d 63 for (int i=fNrows;i--;) delete GetRow(i);
6fefa5ce 64 delete [] fVecs;
8a9ab0eb 65 fNcols = fNrows = 0;
66}
67
68//___________________________________________________________
3bc5dd7d 69void AliMatrixSparse::Print(Option_t* opt) const
8a9ab0eb 70{
3bc5dd7d 71 printf("Sparse Matrix of size %d x %d %s\n",fNrows,fNcols,IsSymmetric() ? " (Symmetric)":"");
72 for (int i=0;i<fNrows;i++) {
8a9ab0eb 73 AliVectorSparse* row = GetRow(i);
74 if (!row->GetNElems()) continue;
75 printf("%3d: ",i);
3bc5dd7d 76 row->Print(opt);
8a9ab0eb 77 }
78}
79
80//___________________________________________________________
551c9e69 81void AliMatrixSparse::MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const
8a9ab0eb 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//___________________________________________________________
115void AliMatrixSparse::SortIndices(Bool_t valuesToo)
116{
de34b538 117 TStopwatch sw;
118 sw.Start();
119 printf("AliMatrixSparse:sort>>\n");
8a9ab0eb 120 // sort columns in increasing order. Used to fix the matrix after ILUk decompostion
121 for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo);
de34b538 122 sw.Stop();
123 sw.Print();
124 printf("AliMatrixSparse:sort<<\n");
125}
126
127//___________________________________________________________
128void 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//___________________________________________________________
157Float_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;
8a9ab0eb 164}
de34b538 165