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