Event server (Mihai)
[u/mrichter/AliRoot.git] / STEER / 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/* */
339fbe23 9/* */
10/* */
de34b538 11/**********************************************************************************************/
8a9ab0eb 12
8a9ab0eb 13//___________________________________________________________
14ClassImp(AliMatrixSparse)
15
16//___________________________________________________________
17AliMatrixSparse::AliMatrixSparse(Int_t sz)
18: AliMatrixSq(),fVecs(0)
19{
339fbe23 20 // constructor
8a9ab0eb 21 fNcols=fNrows=sz;
22 //
23 fVecs = new AliVectorSparse*[sz];
24 for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse();
25}
26
27//___________________________________________________________
28AliMatrixSparse::AliMatrixSparse(const AliMatrixSparse& src)
29 : AliMatrixSq(src),fVecs(0)
30{
339fbe23 31 // copy c-tor
3bc5dd7d 32 fVecs = new AliVectorSparse*[src.GetSize()];
8a9ab0eb 33 for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
34}
35
36//___________________________________________________________
37AliVectorSparse* AliMatrixSparse::GetRowAdd(Int_t ir)
38{
339fbe23 39 // get row, add if needed
3bc5dd7d 40 if (ir>=fNrows) {
8a9ab0eb 41 AliVectorSparse** arrv = new AliVectorSparse*[ir+1];
42 for (int i=GetSize();i--;) arrv[i] = fVecs[i];
6fefa5ce 43 delete [] fVecs;
8a9ab0eb 44 fVecs = arrv;
45 for (int i=GetSize();i<=ir;i++) fVecs[i] = new AliVectorSparse();
3bc5dd7d 46 fNrows = ir+1;
47 if (IsSymmetric() && fNcols<fNrows) fNcols = fNrows;
8a9ab0eb 48 }
49 return fVecs[ir];
50}
51
52//___________________________________________________________
53AliMatrixSparse& AliMatrixSparse::operator=(const AliMatrixSparse& src)
54{
339fbe23 55 // assignment op-r
e99fb5c9 56 if (this == &src) return *this;
57 AliMatrixSq::operator=(src);
58
8a9ab0eb 59 Clear();
3bc5dd7d 60 fNcols = src.GetNCols();
61 fNrows = src.GetNRows();
8a9ab0eb 62 SetSymmetric(src.IsSymmetric());
3bc5dd7d 63 fVecs = new AliVectorSparse*[fNrows];
64 for (int i=fNrows;i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
8a9ab0eb 65 return *this;
66}
67
68//___________________________________________________________
69void AliMatrixSparse::Clear(Option_t*)
70{
339fbe23 71 // clear
3bc5dd7d 72 for (int i=fNrows;i--;) delete GetRow(i);
6fefa5ce 73 delete [] fVecs;
8a9ab0eb 74 fNcols = fNrows = 0;
75}
76
77//___________________________________________________________
3bc5dd7d 78void AliMatrixSparse::Print(Option_t* opt) const
8a9ab0eb 79{
339fbe23 80 // print itself
3bc5dd7d 81 printf("Sparse Matrix of size %d x %d %s\n",fNrows,fNcols,IsSymmetric() ? " (Symmetric)":"");
82 for (int i=0;i<fNrows;i++) {
8a9ab0eb 83 AliVectorSparse* row = GetRow(i);
84 if (!row->GetNElems()) continue;
85 printf("%3d: ",i);
3bc5dd7d 86 row->Print(opt);
8a9ab0eb 87 }
88}
89
90//___________________________________________________________
551c9e69 91void AliMatrixSparse::MultiplyByVec(const Double_t* vecIn, Double_t* vecOut) const
8a9ab0eb 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//___________________________________________________________
125void AliMatrixSparse::SortIndices(Bool_t valuesToo)
126{
339fbe23 127 // sort columns in increasing order. Used to fix the matrix after ILUk decompostion
de34b538 128 TStopwatch sw;
129 sw.Start();
130 printf("AliMatrixSparse:sort>>\n");
8a9ab0eb 131 for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo);
de34b538 132 sw.Stop();
133 sw.Print();
134 printf("AliMatrixSparse:sort<<\n");
135}
136
137//___________________________________________________________
138void 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//___________________________________________________________
167Float_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;
8a9ab0eb 174}
de34b538 175