]>
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 | { | |
3bc5dd7d | 28 | fVecs = new AliVectorSparse*[src.GetSize()]; |
8a9ab0eb | 29 | for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i)); |
30 | } | |
31 | ||
32 | //___________________________________________________________ | |
33 | AliVectorSparse* 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 | //___________________________________________________________ | |
48 | AliMatrixSparse& 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 | //___________________________________________________________ | |
61 | void 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 | 69 | void 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 | 81 | void 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 | //___________________________________________________________ | |
115 | void 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 | //___________________________________________________________ | |
128 | void 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 | //___________________________________________________________ | |
157 | Float_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 |