1 #include "AliMatrixSparse.h"
3 //___________________________________________________________
4 AliVectorSparse::AliVectorSparse()
5 : fNElems(0),fIndex(0),fElems(0) {}
7 //___________________________________________________________
8 AliVectorSparse::AliVectorSparse(const AliVectorSparse& src)
9 : fNElems(src.fNElems),fIndex(0),fElems(0)
11 fIndex = new UShort_t[fNElems];
12 fElems = new Double_t[fNElems];
13 memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
14 memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
17 //___________________________________________________________
18 void AliVectorSparse::Clear()
20 delete[] fIndex; fIndex = 0;
21 delete[] fElems; fElems = 0;
25 //___________________________________________________________
26 Float_t AliMatrixSparse::GetDensity() const
28 // get fraction of non-zero elements
30 for (int i=GetSize();i--;) nel += GetRow(i)->GetNElems();
31 int den = IsSymmetric() ? (GetSize()+1)*GetSize()/2 : GetSize()*GetSize();
32 return float(nel)/den;
35 //___________________________________________________________
36 AliVectorSparse& AliVectorSparse::operator=(const AliVectorSparse& src)
38 if (&src==this) return *this;
40 fNElems = src.fNElems;
41 fIndex = new UShort_t[fNElems];
42 fElems = new Double_t[fNElems];
43 memcpy(fIndex,src.fIndex,fNElems*sizeof(UShort_t));
44 memcpy(fElems,src.fElems,fNElems*sizeof(Double_t));
49 //___________________________________________________________
50 Double_t AliVectorSparse::FindIndex(Int_t ind) const
52 // return an element with given index
53 //printf("V: findindex\n");
57 int mid = (first+last)>>1;
58 if (ind>fIndex[mid]) first = mid+1;
59 else if (ind<fIndex[mid]) last = mid-1;
60 else return fElems[mid];
65 //___________________________________________________________
66 void AliVectorSparse::Zero(Int_t ind)
68 // set element to 0 if it was already defined
72 int mid = (first+last)>>1;
73 if (ind>fIndex[mid]) first = mid+1;
74 else if (ind<fIndex[mid]) last = mid-1;
75 else {fElems[mid] = 0.; return;}
79 //___________________________________________________________
80 Double_t& AliVectorSparse::FindIndexAdd(Int_t ind)
82 // increment an element with given index
83 //printf("V: findindexAdd\n");
87 int mid = (first+last)>>1;
88 if (ind>fIndex[mid]) first = mid+1;
89 else if (ind<fIndex[mid]) last = mid-1;
90 else return fElems[mid];
92 // need to insert a new element
93 UShort_t *arrI = new UShort_t[fNElems+1];
94 memcpy(arrI,fIndex,first*sizeof(UShort_t));
96 memcpy(arrI+first+1,fIndex+first,(fNElems-first)*sizeof(UShort_t));
100 Double_t *arrE = new Double_t[fNElems+1];
101 memcpy(arrE,fElems,first*sizeof(Double_t));
103 memcpy(arrE+first+1,fElems+first,(fNElems-first)*sizeof(Double_t));
108 return fElems[first];
112 //__________________________________________________________
113 void AliVectorSparse::ReSize(Int_t sz,Bool_t copy)
115 if (sz<1) {Clear(); return;}
116 // need to insert a new element
117 UShort_t *arrI = new UShort_t[sz];
118 Double_t *arrE = new Double_t[sz];
119 memset(arrI,0,sz*sizeof(UShort_t));
120 memset(arrE,0,sz*sizeof(Double_t));
122 if (copy && fIndex) {
123 int cpsz = TMath::Min(fNElems,sz);
124 memcpy(arrI,fIndex,cpsz*sizeof(UShort_t));
125 memcpy(arrE,fElems,cpsz*sizeof(Double_t));
135 //__________________________________________________________
136 void AliVectorSparse::SortIndices(Bool_t valuesToo)
138 // sort indices in increasing order. Used to fix the row after ILUk decomposition
139 for (int i=fNElems;i--;) for (int j=i;j--;) if (fIndex[i]<fIndex[j]) { //swap
140 UShort_t tmpI = fIndex[i]; fIndex[i] = fIndex[j]; fIndex[j]=tmpI;
141 if (valuesToo) {Double_t tmpV = fElems[i];fElems[i]=fElems[j];fElems[j]=tmpV;}
145 //__________________________________________________________
146 void AliVectorSparse::Print(Option_t* ) const
149 for (int i=0;i<fNElems;i++) printf("%2d:%+.2e|",fIndex[i],fElems[i]);
154 ///////////////////////////////////////////////////////////////////////////////////////////
155 //___________________________________________________________
156 ClassImp(AliMatrixSparse)
158 //___________________________________________________________
159 AliMatrixSparse::AliMatrixSparse(Int_t sz)
160 : AliMatrixSq(),fVecs(0)
164 fVecs = new AliVectorSparse*[sz];
165 for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse();
168 //___________________________________________________________
169 AliMatrixSparse::AliMatrixSparse(const AliMatrixSparse& src)
170 : AliMatrixSq(src),fVecs(0)
172 fVecs = new AliVectorSparse*[fNcols];
173 for (int i=GetSize();i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
176 //___________________________________________________________
177 AliVectorSparse* AliMatrixSparse::GetRowAdd(Int_t ir)
180 AliVectorSparse** arrv = new AliVectorSparse*[ir+1];
181 for (int i=GetSize();i--;) arrv[i] = fVecs[i];
184 for (int i=GetSize();i<=ir;i++) fVecs[i] = new AliVectorSparse();
185 fNcols = fNrows = ir+1;
190 //___________________________________________________________
191 AliMatrixSparse& AliMatrixSparse::operator=(const AliMatrixSparse& src)
193 if (*this == src) return *this;
195 fNcols = fNrows = src.GetSize();
196 SetSymmetric(src.IsSymmetric());
197 fVecs = new AliVectorSparse*[fNcols];
198 for (int i=fNcols;i--;) fVecs[i] = new AliVectorSparse( *src.GetRow(i));
202 //___________________________________________________________
203 void AliMatrixSparse::Clear(Option_t*)
205 for (int i=fNcols;i--;) delete GetRow(i);
210 //___________________________________________________________
211 void AliMatrixSparse::Print(Option_t*) const
213 printf("Sparse Matrix of size %d %s\n",fNcols,IsSymmetric() ? " (Symmetric)":"");
214 for (int i=0;i<fNcols;i++) {
215 AliVectorSparse* row = GetRow(i);
216 if (!row->GetNElems()) continue;
222 //___________________________________________________________
223 void AliMatrixSparse::MultiplyByVec(Double_t* vecIn, Double_t* vecOut) const
225 // fill vecOut by matrix*vecIn
226 // vector should be of the same size as the matrix
228 memset(vecOut,0,GetSize()*sizeof(Double_t));
230 for (int rw=GetSize();rw--;) { // loop over rows >>>
231 const AliVectorSparse* rowV = GetRow(rw);
232 Int_t nel = rowV->GetNElems();
235 UShort_t* indV = rowV->GetIndices();
236 Double_t* elmV = rowV->GetElems();
239 // treat diagonal term separately. If filled, it should be the last one
240 if (indV[--nel]==rw) vecOut[rw] += vecIn[rw]*elmV[nel];
241 else nel = rowV->GetNElems(); // diag elem was not filled
243 for (int iel=nel;iel--;) { // less element retrieval for symmetric case
245 vecOut[rw] += vecIn[indV[iel]]*elmV[iel];
246 vecOut[indV[iel]] += vecIn[rw]*elmV[iel];
250 else for (int iel=nel;iel--;) if (elmV[iel]) vecOut[rw] += vecIn[indV[iel]]*elmV[iel];
252 } // loop over rows <<<
256 //___________________________________________________________
257 void AliMatrixSparse::SortIndices(Bool_t valuesToo)
259 // sort columns in increasing order. Used to fix the matrix after ILUk decompostion
260 for (int i=GetSize();i--;) GetRow(i)->SortIndices(valuesToo);