]>
Commit | Line | Data |
---|---|---|
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 | //___________________________________________________________ |
14 | ClassImp(AliMatrixSparse) | |
15 | ||
16 | //___________________________________________________________ | |
17 | AliMatrixSparse::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 | //___________________________________________________________ | |
28 | AliMatrixSparse::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 | //___________________________________________________________ | |
37 | AliVectorSparse* 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 | //___________________________________________________________ | |
53 | AliMatrixSparse& 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 | //___________________________________________________________ | |
69 | void 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 | 78 | void 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 | 91 | void 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 | //___________________________________________________________ | |
125 | void 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 | //___________________________________________________________ | |
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; | |
8a9ab0eb | 174 | } |
de34b538 | 175 |