]>
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 | /* */ | |
9 | /**********************************************************************************************/ | |
8a9ab0eb | 10 | |
8a9ab0eb | 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*[fNcols]; | |
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>=fNcols) { | |
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 | fNcols = fNrows = ir+1; | |
42 | } | |
43 | return fVecs[ir]; | |
44 | } | |
45 | ||
46 | //___________________________________________________________ | |
47 | AliMatrixSparse& AliMatrixSparse::operator=(const AliMatrixSparse& src) | |
48 | { | |
49 | if (*this == src) return *this; | |
50 | Clear(); | |
51 | fNcols = fNrows = src.GetSize(); | |
52 | SetSymmetric(src.IsSymmetric()); | |
53 | fVecs = new AliVectorSparse*[fNcols]; | |
54 | for (int i=fNcols;i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i)); | |
55 | return *this; | |
56 | } | |
57 | ||
58 | //___________________________________________________________ | |
59 | void AliMatrixSparse::Clear(Option_t*) | |
60 | { | |
61 | for (int i=fNcols;i--;) delete GetRow(i); | |
62 | delete[] fVecs; | |
63 | fNcols = fNrows = 0; | |
64 | } | |
65 | ||
66 | //___________________________________________________________ | |
67 | void AliMatrixSparse::Print(Option_t*) const | |
68 | { | |
69 | printf("Sparse Matrix of size %d %s\n",fNcols,IsSymmetric() ? " (Symmetric)":""); | |
70 | for (int i=0;i<fNcols;i++) { | |
71 | AliVectorSparse* row = GetRow(i); | |
72 | if (!row->GetNElems()) continue; | |
73 | printf("%3d: ",i); | |
74 | row->Print(); | |
75 | } | |
76 | } | |
77 | ||
78 | //___________________________________________________________ | |
79 | void AliMatrixSparse::MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const | |
80 | { | |
81 | // fill vecOut by matrix*vecIn | |
82 | // vector should be of the same size as the matrix | |
83 | // | |
84 | memset(vecOut,0,GetSize()*sizeof(Double_t)); | |
85 | // | |
86 | for (int rw=GetSize();rw--;) { // loop over rows >>> | |
87 | const AliVectorSparse* rowV = GetRow(rw); | |
88 | Int_t nel = rowV->GetNElems(); | |
89 | if (!nel) continue; | |
90 | // | |
91 | UShort_t* indV = rowV->GetIndices(); | |
92 | Double_t* elmV = rowV->GetElems(); | |
93 | // | |
94 | if (IsSymmetric()) { | |
95 | // treat diagonal term separately. If filled, it should be the last one | |
96 | if (indV[--nel]==rw) vecOut[rw] += vecIn[rw]*elmV[nel]; | |
97 | else nel = rowV->GetNElems(); // diag elem was not filled | |
98 | // | |
99 | for (int iel=nel;iel--;) { // less element retrieval for symmetric case | |
100 | if (elmV[iel]) { | |
101 | vecOut[rw] += vecIn[indV[iel]]*elmV[iel]; | |
102 | vecOut[indV[iel]] += vecIn[rw]*elmV[iel]; | |
103 | } | |
104 | } | |
105 | } | |
106 | else for (int iel=nel;iel--;) if (elmV[iel]) vecOut[rw] += vecIn[indV[iel]]*elmV[iel]; | |
107 | // | |
108 | } // loop over rows <<< | |
109 | // | |
110 | } | |
111 | ||
112 | //___________________________________________________________ | |
113 | void AliMatrixSparse::SortIndices(Bool_t valuesToo) | |
114 | { | |
de34b538 | 115 | TStopwatch sw; |
116 | sw.Start(); | |
117 | printf("AliMatrixSparse:sort>>\n"); | |
8a9ab0eb | 118 | // sort columns in increasing order. Used to fix the matrix after ILUk decompostion |
119 | for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo); | |
de34b538 | 120 | sw.Stop(); |
121 | sw.Print(); | |
122 | printf("AliMatrixSparse:sort<<\n"); | |
123 | } | |
124 | ||
125 | //___________________________________________________________ | |
126 | void AliMatrixSparse::AddToRow(Int_t r, Double_t *valc,Int_t *indc,Int_t n) | |
127 | { | |
128 | // for sym. matrix count how many elems to add have row>=col and assign excplicitly | |
129 | // those which have row<col | |
130 | // | |
131 | // range in increasing order of indices | |
132 | for (int i=n;i--;) { | |
133 | for (int j=i;j>=0;j--) { | |
134 | if (indc[j]>indc[i]) { // swap | |
135 | int ti = indc[i]; indc[i] = indc[j]; indc[j] = ti; | |
136 | double tv = valc[i]; valc[i] = valc[j]; valc[j] = tv; | |
137 | } | |
138 | } | |
139 | } | |
140 | // | |
141 | int ni=n; | |
142 | if (IsSymmetric()) { | |
143 | while(ni--) { | |
144 | if (indc[ni]>r) (*this)(indc[ni],r) += valc[ni]; | |
145 | else break; // use the fact that the indices are ranged in increasing order | |
146 | } | |
147 | } | |
148 | // | |
149 | if (ni<0) return; | |
150 | AliVectorSparse* row = GetRowAdd(r); | |
151 | row->Add(valc,indc,ni+1); | |
152 | } | |
153 | ||
154 | //___________________________________________________________ | |
155 | Float_t AliMatrixSparse::GetDensity() const | |
156 | { | |
157 | // get fraction of non-zero elements | |
158 | Int_t nel = 0; | |
159 | for (int i=GetSize();i--;) nel += GetRow(i)->GetNElems(); | |
160 | int den = IsSymmetric() ? (GetSize()+1)*GetSize()/2 : GetSize()*GetSize(); | |
161 | return float(nel)/den; | |
8a9ab0eb | 162 | } |
de34b538 | 163 |