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