]>
Commit | Line | Data |
---|---|---|
8a9ab0eb | 1 | #include "AliMinResSolve.h" |
2 | #include "AliLog.h" | |
3 | #include <TStopwatch.h> | |
4 | ||
5 | using namespace std; | |
6 | ||
7 | ClassImp(AliMinResSolve) | |
8 | ||
9 | //______________________________________________________ | |
10 | AliMinResSolve::AliMinResSolve() : | |
11 | fSize(0),fPrecon(0),fMatrix(0),fRHS(0), | |
12 | fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0), | |
13 | fPvv(0),fPvz(0),fPhh(0), | |
14 | fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0) | |
15 | {} | |
16 | ||
17 | //______________________________________________________ | |
18 | AliMinResSolve::AliMinResSolve(const AliMinResSolve& src) : | |
19 | TObject(src), | |
20 | fSize(src.fSize),fPrecon(src.fPrecon),fMatrix(src.fMatrix),fRHS(src.fRHS), | |
21 | fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0), | |
22 | fPvv(0),fPvz(0),fPhh(0), | |
23 | fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0) | |
24 | {} | |
25 | ||
26 | //______________________________________________________ | |
27 | AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const TVectorD* rhs) : | |
28 | fSize(mat->GetSize()),fPrecon(0),fMatrix(mat),fRHS(rhs->GetMatrixArray()), | |
29 | fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0), | |
30 | fPvv(0),fPvz(0),fPhh(0), | |
31 | fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0) | |
32 | {} | |
33 | ||
34 | //______________________________________________________ | |
35 | AliMinResSolve::AliMinResSolve(const AliMatrixSq *mat, const double* rhs) : | |
36 | fSize(mat->GetSize()),fPrecon(0),fMatrix(mat),fRHS(rhs), | |
37 | fPVecY(0),fPVecR1(0),fPVecR2(0),fPVecV(0),fPVecW(0),fPVecW1(0),fPVecW2(0), | |
38 | fPvv(0),fPvz(0),fPhh(0), | |
39 | fPrecDiag(0),fPrecAux(0),fMatL(0),fMatU(0) | |
40 | {} | |
41 | ||
42 | //______________________________________________________ | |
43 | AliMinResSolve::~AliMinResSolve() | |
44 | { | |
45 | ClearAux(); | |
46 | } | |
47 | ||
48 | //______________________________________________________ | |
49 | AliMinResSolve& AliMinResSolve::operator=(const AliMinResSolve& src) | |
50 | { | |
51 | if (this != &src) { | |
52 | fSize = src.fSize; | |
53 | fPrecon = src.fPrecon; | |
54 | fMatrix = src.fMatrix; | |
55 | fRHS = src.fRHS; | |
56 | } | |
57 | return *this; | |
58 | } | |
59 | ||
60 | //_______________________________________________________________ | |
61 | Int_t AliMinResSolve::BuildPrecon(Int_t prec) | |
62 | { | |
63 | const Double_t kTiny = 1E-12; | |
64 | fPrecon = prec; | |
65 | // | |
66 | if (fPrecon == kPreconDiag) return 0; // nothing to do | |
67 | // | |
68 | else if (fPrecon == kPreconDILU) { | |
69 | fPrecDiag = new Double_t[ fSize ]; | |
70 | fPrecAux = new Double_t[ fSize ]; | |
71 | // | |
72 | // copy inverse diagonal | |
73 | for (int i=0;i<fSize;i++) fPrecDiag[i] = fMatrix->QuerryDiag(i); | |
74 | // | |
75 | for (int i=0;i<fSize;i++) { | |
76 | fPrecDiag[i] = TMath::Abs(fPrecDiag[i])>kTiny ? 1./fPrecDiag[i] : 1./TMath::Sign(kTiny,fPrecDiag[i]); | |
77 | for (int j=i+1;j<fSize;j++) { | |
78 | double vl = fMatrix->Querry(j,i); | |
79 | if (vl!=0) fPrecDiag[j] -= fPrecDiag[i]*vl*vl; | |
80 | } | |
81 | } | |
82 | return 0; | |
83 | } // precon DILU | |
84 | // | |
85 | if (fPrecon>=kPreconILU0 && fPrecon<=kPreconILU10) { | |
86 | if (fMatrix->InheritsFrom("AliMatrixSparse")) return BuildPreconILUK(fPrecon-kPreconILU0); | |
87 | else return BuildPreconILUKDense(fPrecon-kPreconILU0); | |
88 | } | |
89 | // | |
90 | return -1; | |
91 | } | |
92 | ||
93 | ||
94 | //________________________________ FGMRES METHODS ________________________________ | |
95 | Bool_t AliMinResSolve::SolveFGMRES(TVectorD& VecSol,Int_t precon,int itnlim,double rtol,int nkrylov) | |
96 | { | |
97 | return SolveFGMRES(VecSol.GetMatrixArray(),precon,itnlim,rtol,nkrylov); | |
98 | } | |
99 | ||
100 | //________________________________________________________________________________ | |
101 | Bool_t AliMinResSolve::SolveFGMRES(double* VecSol,Int_t precon,int itnlim,double rtol,int nkrylov) | |
102 | { | |
103 | // Adapted from Y.Saad fgmrs.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/ | |
104 | /*---------------------------------------------------------------------- | |
105 | | *** Preconditioned FGMRES *** | |
106 | +----------------------------------------------------------------------- | |
107 | | This is a simple version of the ARMS preconditioned FGMRES algorithm. | |
108 | +----------------------------------------------------------------------- | |
109 | | Y. S. Dec. 2000. -- Apr. 2008 | |
110 | +----------------------------------------------------------------------- | |
111 | | VecSol = real vector of length n containing an initial guess to the | |
112 | | precon = precondtioner id (0 = no precon) | |
113 | | itnlim = max n of iterations | |
114 | | rtol = tolerance for stopping iteration | |
115 | | nkrylov = N of Krylov vectors to store | |
116 | +---------------------------------------------------------------------*/ | |
117 | int l; | |
118 | double status = kTRUE; | |
119 | double t,beta,eps1=0; | |
120 | const double epsmac = 2.22E-16; | |
121 | // | |
122 | AliInfo(Form("Solution by FGMRes: Preconditioner#%d Max.iter.: %d Tol.: %e NKrylov: %d", | |
123 | precon,itnlim,rtol,nkrylov)); | |
124 | // | |
125 | int its = 0; | |
126 | if (nkrylov<10) { | |
127 | AliInfo(Form("Changing N Krylov vectors from %d to %d",nkrylov,10)); | |
128 | nkrylov = 10; | |
129 | } | |
130 | // | |
131 | if (precon>0) { | |
132 | if (precon>=kPreconsTot) { | |
133 | AliInfo(Form("Unknown preconditioner identifier %d, ignore",precon)); | |
134 | } | |
135 | else { | |
136 | if (BuildPrecon(precon)<0) { | |
137 | ClearAux(); | |
138 | AliInfo("FGMRES failed to build the preconditioner"); | |
139 | return kFALSE; | |
140 | } | |
141 | } | |
142 | } | |
143 | // | |
144 | if (!InitAuxFGMRES(nkrylov)) return kFALSE; | |
145 | // | |
146 | for (l=fSize;l--;) VecSol[l] = 0; | |
147 | // | |
148 | //-------------------- outer loop starts here | |
149 | TStopwatch timer; timer.Start(); | |
150 | while(1) { | |
151 | // | |
152 | //-------------------- compute initial residual vector | |
153 | fMatrix->MultiplyByVec(VecSol,fPvv[0]); | |
154 | for (l=fSize;l--;) fPvv[0][l] = fRHS[l] - fPvv[0][l]; // fPvv[0]= initial residual | |
155 | beta = 0; | |
156 | for (l=fSize;l--;) beta += fPvv[0][l]*fPvv[0][l]; | |
157 | beta = TMath::Sqrt(beta); | |
158 | // | |
159 | if (beta == 0.0) break; // success? | |
160 | t = 1.0 / beta; | |
161 | //-------------------- normalize: fPvv[0] = fPvv[0] / beta | |
162 | for (l=fSize;l--;) fPvv[0][l] *= t; | |
163 | if (its == 0) eps1 = rtol*beta; | |
164 | // | |
165 | // ** initialize 1-st term of rhs of hessenberg system | |
166 | fPVecV[0] = beta; | |
167 | int i = -1; | |
168 | do { | |
169 | i++; | |
170 | its++; | |
171 | int i1 = i+1; | |
172 | // | |
173 | // (Right) Preconditioning Operation z_{j} = M^{-1} v_{j} | |
174 | // | |
175 | if (precon>0) ApplyPrecon( fPvv[i], fPvz[i]); | |
176 | else for (l=fSize;l--;) fPvz[i][l] = fPvv[i][l]; | |
177 | // | |
178 | //-------------------- matvec operation w = A z_{j} = A M^{-1} v_{j} | |
179 | fMatrix->MultiplyByVec(fPvz[i],fPvv[i1]); | |
180 | // | |
181 | // modified gram - schmidt... | |
182 | // h_{i,j} = (w,v_{i}) | |
183 | // w = w - h_{i,j} v_{i} | |
184 | // | |
185 | for (int j=0; j<=i; j++) { | |
186 | for (t=0, l=fSize;l--;) t+= fPvv[j][l]*fPvv[i1][l]; | |
187 | fPhh[i][j] = t; | |
188 | for (l=fSize;l--;) fPvv[i1][l] -= t*fPvv[j][l]; | |
189 | } | |
190 | // -------------------- h_{j+1,j} = ||w||_{2} | |
191 | for (t=0,l=fSize;l--;) t += fPvv[i1][l]*fPvv[i1][l]; t = TMath::Sqrt(t); | |
192 | fPhh[i][i1] = t; | |
193 | if (t != 0.0) for (t=1./t, l=0; l<fSize; l++) fPvv[i1][l] *= t; // v_{j+1} = w / h_{j+1,j} | |
194 | // | |
195 | // done with modified gram schimdt and arnoldi step | |
196 | // now update factorization of fPhh | |
197 | // | |
198 | // perform previous transformations on i-th column of h | |
199 | // | |
200 | for (l=1; l<=i; l++) { | |
201 | int l1 = l-1; | |
202 | t = fPhh[i][l1]; | |
203 | fPhh[i][l1] = fPVecR1[l1]*t + fPVecR2[l1]*fPhh[i][l]; | |
204 | fPhh[i][l] = -fPVecR2[l1]*t + fPVecR1[l1]*fPhh[i][l]; | |
205 | } | |
206 | double gam = TMath::Sqrt( fPhh[i][i]*fPhh[i][i] + fPhh[i][i1]*fPhh[i][i1]); | |
207 | // | |
208 | // if gamma is zero then any small value will do... | |
209 | // will affect only residual estimate | |
210 | if (gam == 0.0) gam = epsmac; | |
211 | // get next plane rotation | |
212 | fPVecR1[i] = fPhh[i][i]/gam; | |
213 | fPVecR2[i] = fPhh[i][i1]/gam; | |
214 | fPVecV[i1] =-fPVecR2[i]*fPVecV[i]; | |
215 | fPVecV[i] *= fPVecR1[i]; | |
216 | // | |
217 | // determine residual norm and test for convergence | |
218 | fPhh[i][i] = fPVecR1[i]*fPhh[i][i] + fPVecR2[i]*fPhh[i][i1]; | |
219 | beta = TMath::Abs(fPVecV[i1]); | |
220 | // | |
221 | } while ( (i < nkrylov-1) && (beta > eps1) && (its < itnlim) ); | |
222 | // | |
223 | // now compute solution. 1st, solve upper triangular system | |
224 | fPVecV[i] = fPVecV[i]/fPhh[i][i]; | |
225 | for (int j=1; j<=i; j++) { | |
226 | int k=i-j; | |
227 | for (t=fPVecV[k],l=k+1; l<=i; l++) t -= fPhh[l][k]*fPVecV[l]; | |
228 | fPVecV[k] = t/fPhh[k][k]; | |
229 | } | |
230 | // -------------------- linear combination of v[i]'s to get sol. | |
231 | for (int j=0; j<=i; j++) for (t=fPVecV[j],l=0; l<fSize; l++) VecSol[l] += t*fPvz[j][l]; | |
232 | // | |
233 | // -------------------- restart outer loop if needed | |
234 | // | |
235 | if (beta <= eps1) { | |
236 | timer.Stop(); | |
237 | AliInfo(Form("FGMRES converged in %d iterations, CPU time: %.1f sec",its,timer.CpuTime())); | |
238 | break; // success | |
239 | } | |
240 | // | |
241 | if (its >= itnlim) { | |
242 | timer.Stop(); | |
243 | AliInfo(Form("%d iterations limit exceeded, CPU time: %.1f sec",itnlim,timer.CpuTime())); | |
244 | status = kFALSE; | |
245 | break; | |
246 | } | |
247 | } | |
248 | // | |
249 | ClearAux(); | |
250 | return status; | |
251 | // | |
252 | } | |
253 | ||
254 | ||
255 | //________________________________ MINRES METHODS ________________________________ | |
256 | Bool_t AliMinResSolve::SolveMinRes(TVectorD& VecSol,Int_t precon,int itnlim,double rtol) | |
257 | { | |
258 | return SolveMinRes(VecSol.GetMatrixArray(), precon, itnlim, rtol); | |
259 | } | |
260 | ||
261 | //________________________________________________________________________________ | |
262 | Bool_t AliMinResSolve::SolveMinRes(double *VecSol,Int_t precon,int itnlim,double rtol) | |
263 | { | |
264 | /* | |
265 | Adapted from author's Fortran code: | |
266 | Michael A. Saunders na.msaunders@na-net.ornl.gov | |
267 | ||
268 | MINRES is an implementation of the algorithm described in the following reference: | |
269 | C. C. Paige and M. A. Saunders (1975), | |
270 | Solution of sparse indefinite systems of linear equations, | |
271 | SIAM J. Numer. Anal. 12(4), pp. 617-629. | |
272 | ||
273 | */ | |
274 | if (!fMatrix->IsSymmetric()) { | |
275 | AliInfo("MinRes cannot solve asymmetric matrices, use FGMRes instead"); | |
276 | return kFALSE; | |
277 | } | |
278 | // | |
279 | ClearAux(); | |
280 | const double eps = 2.22E-16; | |
281 | double beta1; | |
282 | // | |
283 | if (precon>0) { | |
284 | if (precon>=kPreconsTot) { | |
285 | AliInfo(Form("Unknown preconditioner identifier %d, ignore",precon)); | |
286 | } | |
287 | else { | |
288 | if (BuildPrecon(precon)<0) { | |
289 | ClearAux(); | |
290 | AliInfo("MinRes failed to build the preconditioner"); | |
291 | return kFALSE; | |
292 | } | |
293 | } | |
294 | } | |
295 | AliInfo(Form("Solution by MinRes: Preconditioner#%d Max.iter.: %d Tol.: %e",precon,itnlim,rtol)); | |
296 | // | |
297 | // ------------------------ initialization ---------------------->>>> | |
298 | memset(VecSol,0,fSize*sizeof(double)); | |
299 | int status=0, itn=0; | |
300 | double Anorm = 0; | |
301 | double Acond = 0; | |
302 | double ynorm = 0; | |
303 | double rnorm = 0; | |
304 | double gam,gmax=1,gmin=1,gbar,oldeps,epsa,epsx,epsr,diag, delta,phi,denom,z; | |
305 | // | |
306 | if (!InitAuxMinRes()) return kFALSE; | |
307 | // | |
308 | memset(VecSol,0,fSize*sizeof(double)); | |
309 | // | |
310 | // ------------ init aux -------------------------<<<< | |
311 | // Set up y and v for the first Lanczos vector v1. | |
312 | // y = beta1 P' v1, where P = C**(-1). v is really P' v1. | |
313 | // | |
314 | for (int i=fSize;i--;) fPVecY[i] = fPVecR1[i] = fRHS[i]; | |
315 | // | |
316 | if ( precon>0 ) ApplyPrecon( fRHS, fPVecY); | |
317 | beta1 = 0; for (int i=fSize;i--;) beta1 += fRHS[i]*fPVecY[i]; // | |
318 | // | |
319 | if (beta1 < 0) { | |
320 | AliInfo(Form("Preconditioner is indefinite (init) (%e).",beta1)); | |
321 | ClearAux(); | |
322 | status = 7; | |
323 | return kFALSE; | |
324 | } | |
325 | // | |
326 | if (beta1 == 0) { | |
327 | AliInfo(Form("RHS is zero or is the nullspace of the Preconditioner: Solution is {0}")); | |
328 | ClearAux(); | |
329 | return kTRUE; | |
330 | } | |
331 | // | |
332 | beta1 = TMath::Sqrt( beta1 ); // Normalize y to get v1 later. | |
333 | // | |
334 | // See if Msolve is symmetric. //RS: Skept | |
335 | // See if Aprod is symmetric. //RS: Skept | |
336 | // | |
337 | double oldb = 0; | |
338 | double beta = beta1; | |
339 | double dbar = 0; | |
340 | double epsln = 0; | |
341 | double qrnorm = beta1; | |
342 | double phibar = beta1; | |
343 | double rhs1 = beta1; | |
344 | double rhs2 = 0; | |
345 | double tnorm2 = 0; | |
346 | double ynorm2 = 0; | |
347 | double cs = -1; | |
348 | double sn = 0; | |
349 | for (int i=fSize;i--;) fPVecR2[i] = fPVecR1[i]; | |
350 | // | |
351 | TStopwatch timer; timer.Start(); | |
352 | while(status==0) { //----------------- Main iteration loop ---------------------->>>> | |
353 | // | |
354 | itn++; | |
355 | /*----------------------------------------------------------------- | |
356 | Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,... | |
357 | The general iteration is similar to the case k = 1 with v0 = 0: | |
358 | p1 = Operator * v1 - beta1 * v0, | |
359 | alpha1 = v1'p1, | |
360 | q2 = p2 - alpha1 * v1, | |
361 | beta2^2 = q2'q2, | |
362 | v2 = (1/beta2) q2. | |
363 | Again, y = betak P vk, where P = C**(-1). | |
364 | .... more description needed. | |
365 | -----------------------------------------------------------------*/ | |
366 | // | |
367 | double s = 1./beta; // Normalize previous vector (in y). | |
368 | for (int i=fSize;i--;) fPVecV[i] = s*fPVecY[i]; // v = vk if P = I | |
369 | // | |
370 | fMatrix->MultiplyByVec(fPVecV,fPVecY); // APROD (VecV, VecY); | |
371 | // | |
372 | if (itn>=2) { | |
373 | double btrat = beta/oldb; | |
374 | for (int i=fSize;i--;) fPVecY[i] -= btrat*fPVecR1[i]; | |
375 | } | |
376 | double alfa = 0; for (int i=fSize;i--;) alfa += fPVecV[i]*fPVecY[i]; // alphak | |
377 | // | |
378 | double alf2bt = alfa/beta; | |
379 | for (int i=fSize;i--;) { | |
380 | fPVecY[i] -= alf2bt*fPVecR2[i]; | |
381 | fPVecR1[i] = fPVecR2[i]; | |
382 | fPVecR2[i] = fPVecY[i]; | |
383 | } | |
384 | // | |
385 | if ( precon>0 ) ApplyPrecon(fPVecR2, fPVecY); | |
386 | // | |
387 | oldb = beta; // oldb = betak | |
388 | beta = 0; for (int i=fSize;i--;) beta += fPVecR2[i]*fPVecY[i]; // beta = betak+1^2 | |
389 | // | |
390 | if (beta<0) { | |
391 | AliInfo(Form("Preconditioner is indefinite (%e)",beta)); | |
392 | status = 7; | |
393 | break; | |
394 | } | |
395 | // | |
396 | beta = TMath::Sqrt(beta); // beta = betak+1 | |
397 | tnorm2 += alfa*alfa + oldb*oldb + beta*beta; | |
398 | // | |
399 | if (itn == 1) { // Initialize a few things. | |
400 | if (beta/beta1 <= 10.0*eps) { | |
401 | status = 0; //-1 //????? beta2 = 0 or ~ 0, terminate later. | |
402 | AliInfo("RHS is eigenvector"); | |
403 | } | |
404 | // !tnorm2 = alfa**2 | |
405 | gmax = TMath::Abs(alfa); // alpha1 | |
406 | gmin = gmax; // alpha1 | |
407 | } | |
408 | // | |
409 | /* | |
410 | Apply previous rotation Qk-1 to get | |
411 | [deltak epslnk+1] = [cs sn][dbark 0 ] | |
412 | [gbar k dbar k+1] [sn -cs][alfak betak+1]. | |
413 | */ | |
414 | // | |
415 | oldeps = epsln; | |
416 | delta = cs * dbar + sn * alfa; // delta1 = 0 deltak | |
417 | gbar = sn * dbar - cs * alfa; // gbar 1 = alfa1 gbar k | |
418 | epsln = sn * beta; // epsln2 = 0 epslnk+1 | |
419 | dbar = - cs * beta; // dbar 2 = beta2 dbar k+1 | |
420 | // | |
421 | // Compute the next plane rotation Qk | |
422 | // | |
423 | gam = TMath::Sqrt( gbar*gbar + beta*beta ); // gammak | |
424 | cs = gbar / gam; // ck | |
425 | sn = beta / gam; // sk | |
426 | phi = cs * phibar; // phik | |
427 | phibar = sn * phibar; // phibark+1 | |
428 | // | |
429 | // Update x. | |
430 | denom = 1./gam; | |
431 | // | |
432 | for (int i=fSize;i--;) { | |
433 | fPVecW1[i] = fPVecW2[i]; | |
434 | fPVecW2[i] = fPVecW[i]; | |
435 | fPVecW[i] = denom*(fPVecV[i]-oldeps*fPVecW1[i]-delta*fPVecW2[i]); | |
436 | VecSol[i] += phi*fPVecW[i]; | |
437 | } | |
438 | // | |
439 | // Go round again. | |
440 | // | |
441 | gmax = TMath::Max( gmax, gam ); | |
442 | gmin = TMath::Min( gmin, gam ); | |
443 | z = rhs1 / gam; | |
444 | ynorm2 += z*z; | |
445 | rhs1 = rhs2 - delta * z; | |
446 | rhs2 = - epsln * z; | |
447 | // | |
448 | // Estimate various norms and test for convergence. | |
449 | Anorm = TMath::Sqrt( tnorm2 ); | |
450 | ynorm = TMath::Sqrt( ynorm2 ); | |
451 | epsa = Anorm * eps; | |
452 | epsx = Anorm * ynorm * eps; | |
453 | epsr = Anorm * ynorm * rtol; | |
454 | diag = gbar; | |
455 | if (diag == 0) diag = epsa; | |
456 | // | |
457 | qrnorm = phibar; | |
458 | rnorm = qrnorm; | |
459 | /* | |
460 | Estimate cond(A). | |
461 | In this version we look at the diagonals of R in the | |
462 | factorization of the lower Hessenberg matrix, Q * H = R, | |
463 | where H is the tridiagonal matrix from Lanczos with one | |
464 | extra row, beta(k+1) e_k^T. | |
465 | */ | |
466 | Acond = gmax / gmin; | |
467 | // | |
468 | // See if any of the stopping criteria are satisfied. | |
469 | // In rare cases, istop is already -1 from above (Abar = const*I). | |
470 | // | |
471 | if (status == 0) { | |
472 | if (itn >= itnlim ) {status = 5; AliInfo(Form("%d iterations limit exceeded",itnlim));} | |
473 | if (Acond >= 0.1/eps ) {status = 4; AliInfo(Form("Matrix condition nymber %e exceeds limit %e",Acond,0.1/eps));} | |
474 | if (epsx >= beta1 ) {status = 3; AliInfo(Form("Approximate convergence"));} | |
475 | if (qrnorm <= epsx ) {status = 2; AliInfo(Form("Converged within machine precision"));} | |
476 | if (qrnorm <= epsr ) {status = 1; AliInfo(Form("Converged"));} | |
477 | } | |
478 | // | |
479 | } //----------------- Main iteration loop ----------------------<<< | |
480 | // | |
481 | ClearAux(); | |
482 | // | |
483 | timer.Stop(); | |
484 | AliInfo(Form("Exit from MinRes: CPU time: %.2f sec\n" | |
485 | "Status : %2d\n" | |
486 | "Iterations: %4d\n" | |
487 | "Norm : %+e\n" | |
488 | "Condition : %+e\n" | |
489 | "Res.Norm : %+e\n" | |
490 | "Sol.Norm : %+e", | |
491 | timer.CpuTime(),status,itn,Anorm,Acond,rnorm,ynorm)); | |
492 | // | |
493 | return status>=0 && status<=3; | |
494 | // | |
495 | } | |
496 | ||
497 | //______________________________________________________________ | |
498 | void AliMinResSolve::ApplyPrecon(const TVectorD& vecRHS, TVectorD& vecOut) const | |
499 | { | |
500 | ApplyPrecon(vecRHS.GetMatrixArray(), vecOut.GetMatrixArray()); | |
501 | } | |
502 | ||
503 | //______________________________________________________________ | |
504 | void AliMinResSolve::ApplyPrecon(const double* vecRHS, double* vecOut) const | |
505 | { | |
506 | // Application of the preconditioner matrix: | |
507 | // implicitly defines the matrix solving the M*VecOut = VecRHS | |
508 | const Double_t kTiny = 1E-12; | |
509 | if (fPrecon==kPreconDiag) { | |
510 | for (int i=fSize;i--;) { | |
511 | double d = fMatrix->QuerryDiag(i); | |
512 | vecOut[i] = vecRHS[i]/(TMath::Abs(d)>kTiny ? d : kTiny); | |
513 | } | |
514 | // return; | |
515 | } | |
516 | // | |
517 | else if (fPrecon==kPreconDILU) { | |
518 | for (int i=0;i<fSize;i++) { | |
519 | double el = 0; | |
520 | for (int j=i;j--;) {double vl = fMatrix->Querry(i,j);if (vl!=0) el += fPrecAux[j]*vl;} | |
521 | fPrecAux[i] = fPrecDiag[i]*(vecRHS[i]-el); | |
522 | } | |
523 | for (int i=fSize;i--;) { | |
524 | double el = 0; | |
525 | for (int j=i+1;j<fSize;j++) {double vl = fMatrix->Querry(i,j);if (vl!=0) el += vl*vecOut[j];} | |
526 | vecOut[i] = fPrecAux[i] - fPrecDiag[i]*el; | |
527 | } | |
528 | // return; | |
529 | } | |
530 | // | |
531 | else if (fPrecon>=kPreconILU0 && fPrecon<=kPreconILU10) { | |
532 | // | |
533 | for(int i=0; i<fSize; i++) { // Block L solve | |
534 | vecOut[i] = vecRHS[i]; | |
535 | AliVectorSparse &rowLi = *fMatL->GetRow(i); | |
536 | int n = rowLi.GetNElems(); | |
537 | for(int j=0;j<n;j++) vecOut[i] -= vecOut[ rowLi.GetIndex(j) ] * rowLi.GetElem(j); | |
538 | ||
539 | } | |
540 | // | |
541 | for(int i=fSize; i--; ) { // Block -- U solve | |
542 | AliVectorSparse &rowUi = *fMatU->GetRow(i); | |
543 | int n = rowUi.GetNElems(); | |
544 | for(int j=0;j<n;j++ ) vecOut[i] -= vecOut[ rowUi.GetIndex(j) ] * rowUi.GetElem(j); | |
545 | vecOut[i] *= fPrecDiag[i]; | |
546 | } | |
547 | // | |
548 | } | |
549 | // | |
550 | } | |
551 | ||
552 | ||
553 | //___________________________________________________________ | |
554 | Bool_t AliMinResSolve::InitAuxMinRes() | |
555 | { | |
556 | try { | |
557 | fPVecY = new double[fSize]; | |
558 | fPVecR1 = new double[fSize]; | |
559 | fPVecR2 = new double[fSize]; | |
560 | fPVecV = new double[fSize]; | |
561 | fPVecW = new double[fSize]; | |
562 | fPVecW1 = new double[fSize]; | |
563 | fPVecW2 = new double[fSize]; | |
564 | } | |
565 | catch(bad_alloc&) { | |
566 | AliInfo(Form("Failed to allocate the memory for auxialiary arrays for %d equations",fSize)); | |
567 | ClearAux(); | |
568 | return kFALSE; | |
569 | } | |
570 | // | |
571 | for (int i=fSize;i--;) fPVecY[i]=fPVecR1[i]=fPVecR2[i]=fPVecV[i]=fPVecW[i]=fPVecW1[i]=fPVecW2[i]=0.0; | |
572 | // | |
573 | return kTRUE; | |
574 | } | |
575 | ||
576 | ||
577 | //___________________________________________________________ | |
578 | Bool_t AliMinResSolve::InitAuxFGMRES(int nkrylov) | |
579 | { | |
580 | try { | |
581 | fPvv = new double*[nkrylov+1]; | |
582 | fPvz = new double*[nkrylov]; | |
583 | for (int i=0; i<=nkrylov; i++) fPvv[i] = new double[fSize]; | |
584 | fPhh = new double*[nkrylov]; | |
585 | for (int i=0; i<nkrylov; i++) { | |
586 | fPhh[i] = new double[i+2]; | |
587 | fPvz[i] = new double[fSize]; | |
588 | } | |
589 | // | |
590 | fPVecR1 = new double[nkrylov]; | |
591 | fPVecR2 = new double[nkrylov]; | |
592 | fPVecV = new double[nkrylov+1]; | |
593 | } | |
594 | catch(bad_alloc&) { | |
595 | AliInfo(Form("Failed to allocate the memory for auxialiary arrays for %d equations",fSize)); | |
596 | ClearAux(); | |
597 | return kFALSE; | |
598 | } | |
599 | // | |
600 | return kTRUE; | |
601 | } | |
602 | ||
603 | ||
604 | //___________________________________________________________ | |
605 | void AliMinResSolve::ClearAux() | |
606 | { | |
607 | if (fPVecY) delete[] fPVecY; fPVecY = 0; | |
608 | if (fPVecR1) delete[] fPVecR1; fPVecR1 = 0; | |
609 | if (fPVecR2) delete[] fPVecR2; fPVecR2 = 0; | |
610 | if (fPVecV) delete[] fPVecV; fPVecV = 0; | |
611 | if (fPVecW) delete[] fPVecW; fPVecW = 0; | |
612 | if (fPVecW1) delete[] fPVecW1; fPVecW1 = 0; | |
613 | if (fPVecW2) delete[] fPVecW2; fPVecW2 = 0; | |
614 | if (fPrecDiag) delete[] fPrecDiag; fPrecDiag = 0; | |
615 | if (fPrecAux) delete[] fPrecAux; fPrecAux = 0; | |
616 | if (fMatL) delete fMatL; fMatL = 0; | |
617 | if (fMatU) delete fMatU; fMatU = 0; | |
618 | } | |
619 | ||
620 | //___________________________________________________________ | |
621 | Int_t AliMinResSolve::BuildPreconILUK(Int_t lofM) | |
622 | { | |
623 | /*---------------------------------------------------------------------------- | |
624 | * ILUK preconditioner | |
625 | * incomplete LU factorization with level of fill dropping | |
626 | * Adapted from iluk.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/ | |
627 | *----------------------------------------------------------------------------*/ | |
628 | // | |
629 | AliInfo(Form("Building ILU%d preconditioner",lofM)); | |
630 | AliMatrixSparse* Matrix = (AliMatrixSparse*)fMatrix; | |
631 | // | |
632 | TStopwatch sw; sw.Start(); | |
633 | fMatL = new AliMatrixSparse(fSize); | |
634 | fMatU = new AliMatrixSparse(fSize); | |
635 | fMatL->SetSymmetric(kFALSE); | |
636 | fMatU->SetSymmetric(kFALSE); | |
637 | fPrecDiag = new Double_t[fSize]; | |
638 | // | |
639 | // symbolic factorization to calculate level of fill index arrays | |
640 | if ( PreconILUKsymb(lofM)<0 ) { | |
641 | ClearAux(); | |
642 | return -1; | |
643 | } | |
644 | // | |
645 | Int_t *jw = new Int_t[fSize]; | |
646 | for(int j=fSize;j--;) jw[j] = -1; // set indicator array jw to -1 | |
647 | // | |
648 | for(int i=0; i<fSize; i++ ) { // beginning of main loop | |
649 | /* setup array jw[], and initial i-th row */ | |
650 | AliVectorSparse& rowLi = *fMatL->GetRow(i); | |
651 | AliVectorSparse& rowUi = *fMatU->GetRow(i); | |
652 | AliVectorSparse& rowM = *Matrix->GetRow(i); | |
653 | // | |
654 | for(int j=rowLi.GetNElems(); j--;) { // initialize L part | |
655 | int col = rowLi.GetIndex(j); | |
656 | jw[col] = j; | |
657 | rowLi.GetElem(j) = 0.; // do we need this ? | |
658 | } | |
659 | jw[i] = i; | |
660 | fPrecDiag[i] = 0; // initialize diagonal | |
661 | // | |
662 | for(int j=rowUi.GetNElems();j--;) { // initialize U part | |
663 | int col = rowUi.GetIndex(j); | |
664 | jw[col] = j; | |
665 | rowUi.GetElem(j) = 0; | |
666 | } | |
667 | // copy row from csmat into L,U D | |
668 | for(int j=rowM.GetNElems(); j--;) { // L and D part | |
669 | if (rowM.GetElem(j)==0) continue; | |
670 | int col = rowM.GetIndex(j); // (the original matrix stores only lower triangle) | |
671 | if( col < i ) rowLi.GetElem(jw[col]) = rowM.GetElem(j); | |
672 | else if(col==i) fPrecDiag[i] = rowM.GetElem(j); | |
673 | else rowUi.GetElem(jw[col]) = rowM.GetElem(j); | |
674 | } | |
675 | if (fMatrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // part of the row I on the right of diagonal is stored as | |
676 | double vl = fMatrix->Querry(col,i); // the lower part of the column I | |
677 | if (vl==0) continue; | |
678 | rowUi.GetElem(jw[col]) = vl; | |
679 | } | |
680 | // | |
681 | // eliminate previous rows | |
682 | for(int j=0; j<rowLi.GetNElems(); j++) { | |
683 | int jrow = rowLi.GetIndex(j); | |
684 | // get the multiplier for row to be eliminated (jrow) | |
685 | rowLi.GetElem(j) *= fPrecDiag[jrow]; | |
686 | // | |
687 | // combine current row and row jrow | |
688 | AliVectorSparse& rowUj = *fMatU->GetRow(jrow); | |
689 | for(int k=0; k<rowUj.GetNElems(); k++ ) { | |
690 | int col = rowUj.GetIndex(k); | |
691 | int jpos = jw[col]; | |
692 | if( jpos == -1 ) continue; | |
693 | if( col < i ) rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k); | |
694 | else if(col==i) fPrecDiag[i] -= rowLi.GetElem(j) * rowUj.GetElem(k); | |
695 | else rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k); | |
696 | } | |
697 | } | |
698 | // reset double-pointer to -1 ( U-part) | |
699 | for(int j=rowLi.GetNElems(); j--;) jw[ rowLi.GetIndex(j) ] = -1; | |
700 | jw[i] = -1; | |
701 | for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1; | |
702 | // | |
703 | if( fPrecDiag[i] == 0 ) { | |
704 | AliInfo(Form("Fatal error in ILIk: Zero diagonal found...")); | |
705 | delete[] jw; | |
706 | return -1; | |
707 | } | |
708 | fPrecDiag[i] = 1.0 / fPrecDiag[i]; | |
709 | } | |
710 | // | |
711 | delete[] jw; | |
712 | // | |
713 | sw.Stop(); | |
714 | AliInfo(Form("ILU%d preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime())); | |
715 | AliInfo(Form("Densities: M %f L %f U %f",fMatrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity())); | |
716 | // | |
717 | return 0; | |
718 | } | |
719 | ||
720 | //___________________________________________________________ | |
721 | Int_t AliMinResSolve::BuildPreconILUKDense(Int_t lofM) | |
722 | { | |
723 | /*---------------------------------------------------------------------------- | |
724 | * ILUK preconditioner | |
725 | * incomplete LU factorization with level of fill dropping | |
726 | * Adapted from iluk.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/ | |
727 | *----------------------------------------------------------------------------*/ | |
728 | // | |
729 | TStopwatch sw; sw.Start(); | |
730 | AliInfo(Form("Building ILU%d preconditioner for dense matrix",lofM)); | |
731 | // | |
732 | fMatL = new AliMatrixSparse(fSize); | |
733 | fMatU = new AliMatrixSparse(fSize); | |
734 | fMatL->SetSymmetric(kFALSE); | |
735 | fMatU->SetSymmetric(kFALSE); | |
736 | fPrecDiag = new Double_t[fSize]; | |
737 | // | |
738 | // symbolic factorization to calculate level of fill index arrays | |
739 | if ( PreconILUKsymbDense(lofM)<0 ) { | |
740 | ClearAux(); | |
741 | return -1; | |
742 | } | |
743 | // | |
744 | Int_t *jw = new Int_t[fSize]; | |
745 | for(int j=fSize;j--;) jw[j] = -1; // set indicator array jw to -1 | |
746 | // | |
747 | for(int i=0; i<fSize; i++ ) { // beginning of main loop | |
748 | /* setup array jw[], and initial i-th row */ | |
749 | AliVectorSparse& rowLi = *fMatL->GetRow(i); | |
750 | AliVectorSparse& rowUi = *fMatU->GetRow(i); | |
751 | // | |
752 | for(int j=rowLi.GetNElems(); j--;) { // initialize L part | |
753 | int col = rowLi.GetIndex(j); | |
754 | jw[col] = j; | |
755 | rowLi.GetElem(j) = 0.; // do we need this ? | |
756 | } | |
757 | jw[i] = i; | |
758 | fPrecDiag[i] = 0; // initialize diagonal | |
759 | // | |
760 | for(int j=rowUi.GetNElems();j--;) { // initialize U part | |
761 | int col = rowUi.GetIndex(j); | |
762 | jw[col] = j; | |
763 | rowUi.GetElem(j) = 0; | |
764 | } | |
765 | // copy row from csmat into L,U D | |
766 | for(int j=fSize; j--;) { // L and D part | |
767 | double vl = fMatrix->Querry(i,j); | |
768 | if (vl==0) continue; | |
769 | if( j < i ) rowLi.GetElem(jw[j]) = vl; | |
770 | else if(j==i) fPrecDiag[i] = vl; | |
771 | else rowUi.GetElem(jw[j]) = vl; | |
772 | } | |
773 | // eliminate previous rows | |
774 | for(int j=0; j<rowLi.GetNElems(); j++) { | |
775 | int jrow = rowLi.GetIndex(j); | |
776 | // get the multiplier for row to be eliminated (jrow) | |
777 | rowLi.GetElem(j) *= fPrecDiag[jrow]; | |
778 | // | |
779 | // combine current row and row jrow | |
780 | AliVectorSparse& rowUj = *fMatU->GetRow(jrow); | |
781 | for(int k=0; k<rowUj.GetNElems(); k++ ) { | |
782 | int col = rowUj.GetIndex(k); | |
783 | int jpos = jw[col]; | |
784 | if( jpos == -1 ) continue; | |
785 | if( col < i ) rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k); | |
786 | else if(col==i) fPrecDiag[i] -= rowLi.GetElem(j) * rowUj.GetElem(k); | |
787 | else rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k); | |
788 | } | |
789 | } | |
790 | // reset double-pointer to -1 ( U-part) | |
791 | for(int j=rowLi.GetNElems(); j--;) jw[ rowLi.GetIndex(j) ] = -1; | |
792 | jw[i] = -1; | |
793 | for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1; | |
794 | // | |
795 | if( fPrecDiag[i] == 0 ) { | |
796 | AliInfo(Form("Fatal error in ILIk: Zero diagonal found...")); | |
797 | delete[] jw; | |
798 | return -1; | |
799 | } | |
800 | fPrecDiag[i] = 1.0 / fPrecDiag[i]; | |
801 | } | |
802 | // | |
803 | delete[] jw; | |
804 | // | |
805 | sw.Stop(); | |
806 | AliInfo(Form("ILU%d dense preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime())); | |
807 | // AliInfo(Form("Densities: M %f L %f U %f",fMatrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity())); | |
808 | // | |
809 | return 0; | |
810 | } | |
811 | ||
812 | //___________________________________________________________ | |
813 | Int_t AliMinResSolve::PreconILUKsymb(Int_t lofM) | |
814 | { | |
815 | /*---------------------------------------------------------------------------- | |
816 | * ILUK preconditioner | |
817 | * incomplete LU factorization with level of fill dropping | |
818 | * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/ | |
819 | *----------------------------------------------------------------------------*/ | |
820 | // | |
821 | AliMatrixSparse* Matrix = (AliMatrixSparse*)fMatrix; | |
822 | // | |
823 | UChar_t **ulvl=0,*levls=0; | |
824 | UShort_t *jbuf=0; | |
825 | Int_t *iw=0; | |
826 | try { | |
827 | ulvl = new UChar_t*[fSize]; // stores lev-fils for U part of ILU factorization | |
828 | levls = new UChar_t[fSize]; | |
829 | jbuf = new UShort_t[fSize]; | |
830 | iw = new Int_t[fSize]; | |
831 | } | |
832 | // | |
833 | catch(bad_alloc&) { | |
834 | AliInfo("Failed to allocate the memory in AliMinResSolve::PreconILUKsymb"); | |
835 | if (ulvl) delete[] ulvl; | |
836 | if (levls) delete[] levls; | |
837 | if (jbuf) delete[] jbuf; | |
838 | if (iw) delete[] iw; | |
839 | ClearAux(); | |
840 | return -1; | |
841 | } | |
842 | // | |
843 | for(int j=fSize; j--;) iw[j] = -1; // initialize iw | |
844 | for(int i=0; i<fSize; i++) { | |
845 | int incl = 0; | |
846 | int incu = i; | |
847 | AliVectorSparse& row = *Matrix->GetRow(i); | |
848 | // | |
849 | // assign lof = 0 for matrix elements | |
850 | for(int j=0;j<row.GetNElems(); j++) { | |
851 | int col = row.GetIndex(j); | |
852 | if (row.GetElem(j)==0) continue; // !!!! matrix is sparse but sometimes 0 appears | |
853 | if (col<i) { // L-part | |
854 | jbuf[incl] = col; | |
855 | levls[incl] = 0; | |
856 | iw[col] = incl++; | |
857 | } | |
858 | else if (col>i) { // This works only for general matrix | |
859 | jbuf[incu] = col; | |
860 | levls[incu] = 0; | |
861 | iw[col] = incu++; | |
862 | } | |
863 | } | |
864 | if (Matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix | |
865 | if (fMatrix->Querry(col,i)==0) continue; // Due to the symmetry == matrix(i,col) | |
866 | jbuf[incu] = col; | |
867 | levls[incu] = 0; | |
868 | iw[col] = incu++; | |
869 | } | |
870 | // | |
871 | // symbolic k,i,j Gaussian elimination | |
872 | int jpiv = -1; | |
873 | while (++jpiv < incl) { | |
874 | int k = jbuf[jpiv] ; // select leftmost pivot | |
875 | int kmin = k; | |
876 | int jmin = jpiv; | |
877 | for(int j=jpiv+1; j<incl; j++) if( jbuf[j]<kmin ) { kmin = jbuf[j]; jmin = j;} | |
878 | // | |
879 | // ------------------------------------ swap | |
880 | if(jmin!=jpiv) { | |
881 | jbuf[jpiv] = kmin; | |
882 | jbuf[jmin] = k; | |
883 | iw[kmin] = jpiv; | |
884 | iw[k] = jmin; | |
885 | int tj = levls[jpiv] ; | |
886 | levls[jpiv] = levls[jmin]; | |
887 | levls[jmin] = tj; | |
888 | k = kmin; | |
889 | } | |
890 | // ------------------------------------ symbolic linear combinaiton of rows | |
891 | AliVectorSparse& rowU = *fMatU->GetRow(k); | |
892 | for(int j=0; j<rowU.GetNElems(); j++ ) { | |
893 | int col = rowU.GetIndex(j); | |
894 | int it = ulvl[k][j]+levls[jpiv]+1; | |
895 | if( it > lofM ) continue; | |
896 | int ip = iw[col]; | |
897 | if( ip == -1 ) { | |
898 | if( col < i) { | |
899 | jbuf[incl] = col; | |
900 | levls[incl] = it; | |
901 | iw[col] = incl++; | |
902 | } | |
903 | else if( col > i ) { | |
904 | jbuf[incu] = col; | |
905 | levls[incu] = it; | |
906 | iw[col] = incu++; | |
907 | } | |
908 | } | |
909 | else levls[ip] = TMath::Min(levls[ip], it); | |
910 | } | |
911 | // | |
912 | } // end - while loop | |
913 | // | |
914 | // reset iw | |
915 | for (int j=0;j<incl;j++) iw[jbuf[j]] = -1; | |
916 | for (int j=i;j<incu;j++) iw[jbuf[j]] = -1; | |
917 | // | |
918 | // copy L-part | |
919 | AliVectorSparse& rowLi = *fMatL->GetRow(i); | |
920 | rowLi.ReSize(incl); | |
921 | if(incl>0) memcpy(rowLi.GetIndices(), jbuf, sizeof(UShort_t)*incl); | |
922 | // copy U-part | |
923 | int k = incu-i; | |
924 | AliVectorSparse& rowUi = *fMatU->GetRow(i); | |
925 | rowUi.ReSize(k); | |
926 | if( k > 0 ) { | |
927 | memcpy(rowUi.GetIndices(), jbuf+i, sizeof(UShort_t)*k); | |
928 | ulvl[i] = new UChar_t[k]; // update matrix of levels | |
929 | memcpy( ulvl[i], levls+i, k*sizeof(UChar_t) ); | |
930 | } | |
931 | } | |
932 | // | |
933 | // free temp space and leave | |
934 | delete[] levls; | |
935 | delete[] jbuf; | |
936 | for(int i=fSize; i--;) if (fMatU->GetRow(i)->GetNElems()) delete[] ulvl[i]; | |
937 | delete[] ulvl; | |
938 | // | |
939 | fMatL->SortIndices(); | |
940 | fMatU->SortIndices(); | |
941 | return 0; | |
942 | } | |
943 | ||
944 | ||
945 | //___________________________________________________________ | |
946 | Int_t AliMinResSolve::PreconILUKsymbDense(Int_t lofM) | |
947 | { | |
948 | /*---------------------------------------------------------------------------- | |
949 | * ILUK preconditioner | |
950 | * incomplete LU factorization with level of fill dropping | |
951 | * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/ | |
952 | *----------------------------------------------------------------------------*/ | |
953 | // | |
954 | UChar_t **ulvl=0,*levls=0; | |
955 | UShort_t *jbuf=0; | |
956 | Int_t *iw=0; | |
957 | try { | |
958 | ulvl = new UChar_t*[fSize]; // stores lev-fils for U part of ILU factorization | |
959 | levls = new UChar_t[fSize]; | |
960 | jbuf = new UShort_t[fSize]; | |
961 | iw = new Int_t[fSize]; | |
962 | } | |
963 | // | |
964 | catch(bad_alloc&) { | |
965 | AliInfo("Failed to allocate the memory in AliMinResSolve::PreconILUKsymb"); | |
966 | if (ulvl) delete[] ulvl; | |
967 | if (levls) delete[] levls; | |
968 | if (jbuf) delete[] jbuf; | |
969 | if (iw) delete[] iw; | |
970 | ClearAux(); | |
971 | return -1; | |
972 | } | |
973 | // | |
974 | for(int j=fSize; j--;) iw[j] = -1; // initialize iw | |
975 | for(int i=0; i<fSize; i++) { | |
976 | int incl = 0; | |
977 | int incu = i; | |
978 | // | |
979 | // assign lof = 0 for matrix elements | |
980 | for(int j=0;j<fSize; j++) { | |
981 | if (fMatrix->Querry(i,j)==0) continue; | |
982 | if (j<i) { // L-part | |
983 | jbuf[incl] = j; | |
984 | levls[incl] = 0; | |
985 | iw[j] = incl++; | |
986 | } | |
987 | else if (j>i) { // This works only for general matrix | |
988 | jbuf[incu] = j; | |
989 | levls[incu] = 0; | |
990 | iw[j] = incu++; | |
991 | } | |
992 | } | |
993 | // | |
994 | // symbolic k,i,j Gaussian elimination | |
995 | int jpiv = -1; | |
996 | while (++jpiv < incl) { | |
997 | int k = jbuf[jpiv] ; // select leftmost pivot | |
998 | int kmin = k; | |
999 | int jmin = jpiv; | |
1000 | for(int j=jpiv+1; j<incl; j++) if( jbuf[j]<kmin ) { kmin = jbuf[j]; jmin = j;} | |
1001 | // | |
1002 | // ------------------------------------ swap | |
1003 | if(jmin!=jpiv) { | |
1004 | jbuf[jpiv] = kmin; | |
1005 | jbuf[jmin] = k; | |
1006 | iw[kmin] = jpiv; | |
1007 | iw[k] = jmin; | |
1008 | int tj = levls[jpiv] ; | |
1009 | levls[jpiv] = levls[jmin]; | |
1010 | levls[jmin] = tj; | |
1011 | k = kmin; | |
1012 | } | |
1013 | // ------------------------------------ symbolic linear combinaiton of rows | |
1014 | AliVectorSparse& rowU = *fMatU->GetRow(k); | |
1015 | for(int j=0; j<rowU.GetNElems(); j++ ) { | |
1016 | int col = rowU.GetIndex(j); | |
1017 | int it = ulvl[k][j]+levls[jpiv]+1; | |
1018 | if( it > lofM ) continue; | |
1019 | int ip = iw[col]; | |
1020 | if( ip == -1 ) { | |
1021 | if( col < i) { | |
1022 | jbuf[incl] = col; | |
1023 | levls[incl] = it; | |
1024 | iw[col] = incl++; | |
1025 | } | |
1026 | else if( col > i ) { | |
1027 | jbuf[incu] = col; | |
1028 | levls[incu] = it; | |
1029 | iw[col] = incu++; | |
1030 | } | |
1031 | } | |
1032 | else levls[ip] = TMath::Min(levls[ip], it); | |
1033 | } | |
1034 | // | |
1035 | } // end - while loop | |
1036 | // | |
1037 | // reset iw | |
1038 | for (int j=0;j<incl;j++) iw[jbuf[j]] = -1; | |
1039 | for (int j=i;j<incu;j++) iw[jbuf[j]] = -1; | |
1040 | // | |
1041 | // copy L-part | |
1042 | AliVectorSparse& rowLi = *fMatL->GetRow(i); | |
1043 | rowLi.ReSize(incl); | |
1044 | if(incl>0) memcpy(rowLi.GetIndices(), jbuf, sizeof(UShort_t)*incl); | |
1045 | // copy U-part | |
1046 | int k = incu-i; | |
1047 | AliVectorSparse& rowUi = *fMatU->GetRow(i); | |
1048 | rowUi.ReSize(k); | |
1049 | if( k > 0 ) { | |
1050 | memcpy(rowUi.GetIndices(), jbuf+i, sizeof(UShort_t)*k); | |
1051 | ulvl[i] = new UChar_t[k]; // update matrix of levels | |
1052 | memcpy( ulvl[i], levls+i, k*sizeof(UChar_t) ); | |
1053 | } | |
1054 | } | |
1055 | // | |
1056 | // free temp space and leave | |
1057 | delete[] levls; | |
1058 | delete[] jbuf; | |
1059 | for(int i=fSize; i--;) if (fMatU->GetRow(i)->GetNElems()) delete[] ulvl[i]; | |
1060 | delete[] ulvl; | |
1061 | // | |
1062 | fMatL->SortIndices(); | |
1063 | fMatU->SortIndices(); | |
1064 | return 0; | |
1065 | } |