]>
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 |
3c5b4cc8 | 554 | fPVecY = new double[fSize]; |
555 | fPVecR1 = new double[fSize]; | |
556 | fPVecR2 = new double[fSize]; | |
557 | fPVecV = new double[fSize]; | |
558 | fPVecW = new double[fSize]; | |
559 | fPVecW1 = new double[fSize]; | |
560 | fPVecW2 = new double[fSize]; | |
8a9ab0eb | 561 | // |
562 | for (int i=fSize;i--;) fPVecY[i]=fPVecR1[i]=fPVecR2[i]=fPVecV[i]=fPVecW[i]=fPVecW1[i]=fPVecW2[i]=0.0; | |
563 | // | |
564 | return kTRUE; | |
565 | } | |
566 | ||
567 | ||
568 | //___________________________________________________________ | |
569 | Bool_t AliMinResSolve::InitAuxFGMRES(int nkrylov) | |
570 | { | |
7c3070ec | 571 | // init auxiliary space for fgmres |
3c5b4cc8 | 572 | fPvv = new double*[nkrylov+1]; |
573 | fPvz = new double*[nkrylov]; | |
574 | for (int i=0; i<=nkrylov; i++) fPvv[i] = new double[fSize]; | |
575 | fPhh = new double*[nkrylov]; | |
576 | for (int i=0; i<nkrylov; i++) { | |
577 | fPhh[i] = new double[i+2]; | |
578 | fPvz[i] = new double[fSize]; | |
8a9ab0eb | 579 | } |
580 | // | |
3c5b4cc8 | 581 | fPVecR1 = new double[nkrylov]; |
582 | fPVecR2 = new double[nkrylov]; | |
583 | fPVecV = new double[nkrylov+1]; | |
584 | // | |
8a9ab0eb | 585 | return kTRUE; |
586 | } | |
587 | ||
588 | ||
589 | //___________________________________________________________ | |
590 | void AliMinResSolve::ClearAux() | |
591 | { | |
7c3070ec | 592 | // clear aux. space |
8a9ab0eb | 593 | if (fPVecY) delete[] fPVecY; fPVecY = 0; |
594 | if (fPVecR1) delete[] fPVecR1; fPVecR1 = 0; | |
595 | if (fPVecR2) delete[] fPVecR2; fPVecR2 = 0; | |
596 | if (fPVecV) delete[] fPVecV; fPVecV = 0; | |
597 | if (fPVecW) delete[] fPVecW; fPVecW = 0; | |
598 | if (fPVecW1) delete[] fPVecW1; fPVecW1 = 0; | |
599 | if (fPVecW2) delete[] fPVecW2; fPVecW2 = 0; | |
7c3070ec | 600 | if (fDiagLU) delete[] fDiagLU; fDiagLU = 0; |
8a9ab0eb | 601 | if (fMatL) delete fMatL; fMatL = 0; |
602 | if (fMatU) delete fMatU; fMatU = 0; | |
7c3070ec | 603 | if (fMatBD) delete fMatBD; fMatBD = 0; |
604 | } | |
605 | ||
606 | //___________________________________________________________ | |
607 | Int_t AliMinResSolve::BuildPreconBD(Int_t hwidth) | |
608 | { | |
609 | // build preconditioner | |
610 | AliInfo(Form("Building Band-Diagonal preconditioner of half-width = %d",hwidth)); | |
611 | fMatBD = new AliSymBDMatrix( fMatrix->GetSize(), hwidth ); | |
612 | // | |
613 | // fill the band-diagonal part of the matrix | |
614 | if (fMatrix->InheritsFrom("AliMatrixSparse")) { | |
615 | for (int ir=fMatrix->GetSize();ir--;) { | |
616 | int jmin = TMath::Max(0,ir-hwidth); | |
617 | AliVectorSparse& irow = *((AliMatrixSparse*)fMatrix)->GetRow(ir); | |
618 | for(int j=irow.GetNElems();j--;) { | |
619 | int jind = irow.GetIndex(j); | |
620 | if (jind<jmin) break; | |
621 | (*fMatBD)(ir,jind) = irow.GetElem(j); | |
622 | } | |
623 | } | |
624 | } | |
625 | else { | |
626 | for (int ir=fMatrix->GetSize();ir--;) { | |
627 | int jmin = TMath::Max(0,ir-hwidth); | |
628 | for(int jr=jmin;jr<=ir;jr++) (*fMatBD)(ir,jr) = fMatrix->Query(ir,jr); | |
629 | } | |
630 | } | |
631 | // | |
632 | fMatBD->DecomposeLDLT(); | |
633 | // | |
634 | return 0; | |
8a9ab0eb | 635 | } |
636 | ||
637 | //___________________________________________________________ | |
638 | Int_t AliMinResSolve::BuildPreconILUK(Int_t lofM) | |
639 | { | |
640 | /*---------------------------------------------------------------------------- | |
641 | * ILUK preconditioner | |
642 | * incomplete LU factorization with level of fill dropping | |
643 | * Adapted from iluk.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/ | |
644 | *----------------------------------------------------------------------------*/ | |
645 | // | |
646 | AliInfo(Form("Building ILU%d preconditioner",lofM)); | |
8a9ab0eb | 647 | // |
648 | TStopwatch sw; sw.Start(); | |
649 | fMatL = new AliMatrixSparse(fSize); | |
650 | fMatU = new AliMatrixSparse(fSize); | |
651 | fMatL->SetSymmetric(kFALSE); | |
652 | fMatU->SetSymmetric(kFALSE); | |
7c3070ec | 653 | fDiagLU = new Double_t[fSize]; |
654 | AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix; | |
8a9ab0eb | 655 | // |
656 | // symbolic factorization to calculate level of fill index arrays | |
657 | if ( PreconILUKsymb(lofM)<0 ) { | |
658 | ClearAux(); | |
659 | return -1; | |
660 | } | |
661 | // | |
662 | Int_t *jw = new Int_t[fSize]; | |
663 | for(int j=fSize;j--;) jw[j] = -1; // set indicator array jw to -1 | |
664 | // | |
665 | for(int i=0; i<fSize; i++ ) { // beginning of main loop | |
de34b538 | 666 | if ( (i%int(0.1*fSize)) == 0) { |
7c3070ec | 667 | AliInfo(Form("BuildPrecon: row %d of %d",i,fSize)); |
de34b538 | 668 | sw.Stop(); |
669 | sw.Print(); | |
670 | sw.Start(kFALSE); | |
671 | } | |
8a9ab0eb | 672 | /* setup array jw[], and initial i-th row */ |
673 | AliVectorSparse& rowLi = *fMatL->GetRow(i); | |
674 | AliVectorSparse& rowUi = *fMatU->GetRow(i); | |
7c3070ec | 675 | AliVectorSparse& rowM = *matrix->GetRow(i); |
8a9ab0eb | 676 | // |
677 | for(int j=rowLi.GetNElems(); j--;) { // initialize L part | |
678 | int col = rowLi.GetIndex(j); | |
679 | jw[col] = j; | |
680 | rowLi.GetElem(j) = 0.; // do we need this ? | |
681 | } | |
682 | jw[i] = i; | |
7c3070ec | 683 | fDiagLU[i] = 0; // initialize diagonal |
8a9ab0eb | 684 | // |
685 | for(int j=rowUi.GetNElems();j--;) { // initialize U part | |
686 | int col = rowUi.GetIndex(j); | |
687 | jw[col] = j; | |
688 | rowUi.GetElem(j) = 0; | |
689 | } | |
690 | // copy row from csmat into L,U D | |
691 | for(int j=rowM.GetNElems(); j--;) { // L and D part | |
7c185459 | 692 | if (AliMatrixSq::IsZero(rowM.GetElem(j))) continue; |
8a9ab0eb | 693 | int col = rowM.GetIndex(j); // (the original matrix stores only lower triangle) |
694 | if( col < i ) rowLi.GetElem(jw[col]) = rowM.GetElem(j); | |
7c3070ec | 695 | else if(col==i) fDiagLU[i] = rowM.GetElem(j); |
8a9ab0eb | 696 | else rowUi.GetElem(jw[col]) = rowM.GetElem(j); |
697 | } | |
7c3070ec | 698 | if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // part of the row I on the right of diagonal is stored as |
699 | double vl = matrix->Query(col,i); // the lower part of the column I | |
7c185459 | 700 | if (AliMatrixSq::IsZero(vl)) continue; |
8a9ab0eb | 701 | rowUi.GetElem(jw[col]) = vl; |
702 | } | |
703 | // | |
704 | // eliminate previous rows | |
705 | for(int j=0; j<rowLi.GetNElems(); j++) { | |
706 | int jrow = rowLi.GetIndex(j); | |
707 | // get the multiplier for row to be eliminated (jrow) | |
7c3070ec | 708 | rowLi.GetElem(j) *= fDiagLU[jrow]; |
8a9ab0eb | 709 | // |
710 | // combine current row and row jrow | |
711 | AliVectorSparse& rowUj = *fMatU->GetRow(jrow); | |
712 | for(int k=0; k<rowUj.GetNElems(); k++ ) { | |
713 | int col = rowUj.GetIndex(k); | |
714 | int jpos = jw[col]; | |
715 | if( jpos == -1 ) continue; | |
716 | if( col < i ) rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k); | |
7c3070ec | 717 | else if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k); |
8a9ab0eb | 718 | else rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k); |
719 | } | |
720 | } | |
721 | // reset double-pointer to -1 ( U-part) | |
722 | for(int j=rowLi.GetNElems(); j--;) jw[ rowLi.GetIndex(j) ] = -1; | |
723 | jw[i] = -1; | |
724 | for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1; | |
725 | // | |
7c185459 | 726 | if( AliMatrixSq::IsZero(fDiagLU[i]) ) { |
8a9ab0eb | 727 | AliInfo(Form("Fatal error in ILIk: Zero diagonal found...")); |
728 | delete[] jw; | |
729 | return -1; | |
730 | } | |
7c3070ec | 731 | fDiagLU[i] = 1.0 / fDiagLU[i]; |
8a9ab0eb | 732 | } |
733 | // | |
734 | delete[] jw; | |
735 | // | |
736 | sw.Stop(); | |
737 | AliInfo(Form("ILU%d preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime())); | |
7c3070ec | 738 | AliInfo(Form("Densities: M %f L %f U %f",matrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity())); |
8a9ab0eb | 739 | // |
740 | return 0; | |
741 | } | |
742 | ||
743 | //___________________________________________________________ | |
744 | Int_t AliMinResSolve::BuildPreconILUKDense(Int_t lofM) | |
745 | { | |
746 | /*---------------------------------------------------------------------------- | |
747 | * ILUK preconditioner | |
748 | * incomplete LU factorization with level of fill dropping | |
749 | * Adapted from iluk.c of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/ | |
750 | *----------------------------------------------------------------------------*/ | |
751 | // | |
752 | TStopwatch sw; sw.Start(); | |
753 | AliInfo(Form("Building ILU%d preconditioner for dense matrix",lofM)); | |
754 | // | |
755 | fMatL = new AliMatrixSparse(fSize); | |
756 | fMatU = new AliMatrixSparse(fSize); | |
757 | fMatL->SetSymmetric(kFALSE); | |
758 | fMatU->SetSymmetric(kFALSE); | |
7c3070ec | 759 | fDiagLU = new Double_t[fSize]; |
8a9ab0eb | 760 | // |
761 | // symbolic factorization to calculate level of fill index arrays | |
762 | if ( PreconILUKsymbDense(lofM)<0 ) { | |
763 | ClearAux(); | |
764 | return -1; | |
765 | } | |
766 | // | |
767 | Int_t *jw = new Int_t[fSize]; | |
768 | for(int j=fSize;j--;) jw[j] = -1; // set indicator array jw to -1 | |
769 | // | |
770 | for(int i=0; i<fSize; i++ ) { // beginning of main loop | |
771 | /* setup array jw[], and initial i-th row */ | |
772 | AliVectorSparse& rowLi = *fMatL->GetRow(i); | |
773 | AliVectorSparse& rowUi = *fMatU->GetRow(i); | |
774 | // | |
775 | for(int j=rowLi.GetNElems(); j--;) { // initialize L part | |
776 | int col = rowLi.GetIndex(j); | |
777 | jw[col] = j; | |
778 | rowLi.GetElem(j) = 0.; // do we need this ? | |
779 | } | |
780 | jw[i] = i; | |
7c3070ec | 781 | fDiagLU[i] = 0; // initialize diagonal |
8a9ab0eb | 782 | // |
783 | for(int j=rowUi.GetNElems();j--;) { // initialize U part | |
784 | int col = rowUi.GetIndex(j); | |
785 | jw[col] = j; | |
786 | rowUi.GetElem(j) = 0; | |
787 | } | |
788 | // copy row from csmat into L,U D | |
789 | for(int j=fSize; j--;) { // L and D part | |
de34b538 | 790 | double vl = fMatrix->Query(i,j); |
7c185459 | 791 | if (AliMatrixSq::IsZero(vl)) continue; |
8a9ab0eb | 792 | if( j < i ) rowLi.GetElem(jw[j]) = vl; |
7c3070ec | 793 | else if(j==i) fDiagLU[i] = vl; |
8a9ab0eb | 794 | else rowUi.GetElem(jw[j]) = vl; |
795 | } | |
796 | // eliminate previous rows | |
797 | for(int j=0; j<rowLi.GetNElems(); j++) { | |
798 | int jrow = rowLi.GetIndex(j); | |
799 | // get the multiplier for row to be eliminated (jrow) | |
7c3070ec | 800 | rowLi.GetElem(j) *= fDiagLU[jrow]; |
8a9ab0eb | 801 | // |
802 | // combine current row and row jrow | |
803 | AliVectorSparse& rowUj = *fMatU->GetRow(jrow); | |
804 | for(int k=0; k<rowUj.GetNElems(); k++ ) { | |
805 | int col = rowUj.GetIndex(k); | |
806 | int jpos = jw[col]; | |
807 | if( jpos == -1 ) continue; | |
808 | if( col < i ) rowLi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k); | |
7c3070ec | 809 | else if(col==i) fDiagLU[i] -= rowLi.GetElem(j) * rowUj.GetElem(k); |
8a9ab0eb | 810 | else rowUi.GetElem(jpos) -= rowLi.GetElem(j) * rowUj.GetElem(k); |
811 | } | |
812 | } | |
813 | // reset double-pointer to -1 ( U-part) | |
814 | for(int j=rowLi.GetNElems(); j--;) jw[ rowLi.GetIndex(j) ] = -1; | |
815 | jw[i] = -1; | |
816 | for(int j=rowUi.GetNElems(); j--;) jw[ rowUi.GetIndex(j) ] = -1; | |
817 | // | |
7c185459 | 818 | if( AliMatrixSq::IsZero(fDiagLU[i])) { |
8a9ab0eb | 819 | AliInfo(Form("Fatal error in ILIk: Zero diagonal found...")); |
820 | delete[] jw; | |
821 | return -1; | |
822 | } | |
7c3070ec | 823 | fDiagLU[i] = 1.0 / fDiagLU[i]; |
8a9ab0eb | 824 | } |
825 | // | |
826 | delete[] jw; | |
827 | // | |
828 | sw.Stop(); | |
829 | AliInfo(Form("ILU%d dense preconditioner OK, CPU time: %.1f sec",lofM,sw.CpuTime())); | |
830 | // AliInfo(Form("Densities: M %f L %f U %f",fMatrix->GetDensity(),fMatL->GetDensity(),fMatU->GetDensity())); | |
831 | // | |
832 | return 0; | |
833 | } | |
834 | ||
835 | //___________________________________________________________ | |
836 | Int_t AliMinResSolve::PreconILUKsymb(Int_t lofM) | |
837 | { | |
838 | /*---------------------------------------------------------------------------- | |
839 | * ILUK preconditioner | |
840 | * incomplete LU factorization with level of fill dropping | |
841 | * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/ | |
842 | *----------------------------------------------------------------------------*/ | |
843 | // | |
de34b538 | 844 | TStopwatch sw; |
7c3070ec | 845 | AliInfo("PreconILUKsymb>>"); |
846 | AliMatrixSparse* matrix = (AliMatrixSparse*)fMatrix; | |
de34b538 | 847 | sw.Start(); |
8a9ab0eb | 848 | // |
849 | UChar_t **ulvl=0,*levls=0; | |
850 | UShort_t *jbuf=0; | |
851 | Int_t *iw=0; | |
3c5b4cc8 | 852 | ulvl = new UChar_t*[fSize]; // stores lev-fils for U part of ILU factorization |
853 | levls = new UChar_t[fSize]; | |
854 | jbuf = new UShort_t[fSize]; | |
855 | iw = new Int_t[fSize]; | |
8a9ab0eb | 856 | // |
857 | for(int j=fSize; j--;) iw[j] = -1; // initialize iw | |
858 | for(int i=0; i<fSize; i++) { | |
859 | int incl = 0; | |
860 | int incu = i; | |
7c3070ec | 861 | AliVectorSparse& row = *matrix->GetRow(i); |
8a9ab0eb | 862 | // |
863 | // assign lof = 0 for matrix elements | |
864 | for(int j=0;j<row.GetNElems(); j++) { | |
865 | int col = row.GetIndex(j); | |
7c185459 | 866 | if (AliMatrixSq::IsZero(row.GetElem(j))) continue; // !!!! matrix is sparse but sometimes 0 appears |
8a9ab0eb | 867 | if (col<i) { // L-part |
868 | jbuf[incl] = col; | |
869 | levls[incl] = 0; | |
870 | iw[col] = incl++; | |
871 | } | |
872 | else if (col>i) { // This works only for general matrix | |
873 | jbuf[incu] = col; | |
874 | levls[incu] = 0; | |
875 | iw[col] = incu++; | |
876 | } | |
877 | } | |
7c3070ec | 878 | if (matrix->IsSymmetric()) for (int col=i+1;col<fSize;col++) { // U-part of symmetric matrix |
7c185459 | 879 | if (AliMatrixSq::IsZero(matrix->Query(col,i))) continue; // Due to the symmetry == matrix(i,col) |
8a9ab0eb | 880 | jbuf[incu] = col; |
881 | levls[incu] = 0; | |
882 | iw[col] = incu++; | |
883 | } | |
884 | // | |
885 | // symbolic k,i,j Gaussian elimination | |
886 | int jpiv = -1; | |
887 | while (++jpiv < incl) { | |
888 | int k = jbuf[jpiv] ; // select leftmost pivot | |
889 | int kmin = k; | |
890 | int jmin = jpiv; | |
891 | for(int j=jpiv+1; j<incl; j++) if( jbuf[j]<kmin ) { kmin = jbuf[j]; jmin = j;} | |
892 | // | |
893 | // ------------------------------------ swap | |
894 | if(jmin!=jpiv) { | |
895 | jbuf[jpiv] = kmin; | |
896 | jbuf[jmin] = k; | |
897 | iw[kmin] = jpiv; | |
898 | iw[k] = jmin; | |
899 | int tj = levls[jpiv] ; | |
900 | levls[jpiv] = levls[jmin]; | |
901 | levls[jmin] = tj; | |
902 | k = kmin; | |
903 | } | |
904 | // ------------------------------------ symbolic linear combinaiton of rows | |
905 | AliVectorSparse& rowU = *fMatU->GetRow(k); | |
906 | for(int j=0; j<rowU.GetNElems(); j++ ) { | |
907 | int col = rowU.GetIndex(j); | |
908 | int it = ulvl[k][j]+levls[jpiv]+1; | |
909 | if( it > lofM ) continue; | |
910 | int ip = iw[col]; | |
911 | if( ip == -1 ) { | |
912 | if( col < i) { | |
913 | jbuf[incl] = col; | |
914 | levls[incl] = it; | |
915 | iw[col] = incl++; | |
916 | } | |
917 | else if( col > i ) { | |
918 | jbuf[incu] = col; | |
919 | levls[incu] = it; | |
920 | iw[col] = incu++; | |
921 | } | |
922 | } | |
923 | else levls[ip] = TMath::Min(levls[ip], it); | |
924 | } | |
925 | // | |
926 | } // end - while loop | |
927 | // | |
928 | // reset iw | |
929 | for (int j=0;j<incl;j++) iw[jbuf[j]] = -1; | |
930 | for (int j=i;j<incu;j++) iw[jbuf[j]] = -1; | |
931 | // | |
932 | // copy L-part | |
933 | AliVectorSparse& rowLi = *fMatL->GetRow(i); | |
934 | rowLi.ReSize(incl); | |
935 | if(incl>0) memcpy(rowLi.GetIndices(), jbuf, sizeof(UShort_t)*incl); | |
936 | // copy U-part | |
937 | int k = incu-i; | |
938 | AliVectorSparse& rowUi = *fMatU->GetRow(i); | |
939 | rowUi.ReSize(k); | |
940 | if( k > 0 ) { | |
941 | memcpy(rowUi.GetIndices(), jbuf+i, sizeof(UShort_t)*k); | |
942 | ulvl[i] = new UChar_t[k]; // update matrix of levels | |
943 | memcpy( ulvl[i], levls+i, k*sizeof(UChar_t) ); | |
944 | } | |
945 | } | |
946 | // | |
947 | // free temp space and leave | |
948 | delete[] levls; | |
949 | delete[] jbuf; | |
950 | for(int i=fSize; i--;) if (fMatU->GetRow(i)->GetNElems()) delete[] ulvl[i]; | |
951 | delete[] ulvl; | |
f8b1c575 | 952 | delete[] iw; |
8a9ab0eb | 953 | // |
954 | fMatL->SortIndices(); | |
955 | fMatU->SortIndices(); | |
de34b538 | 956 | sw.Stop(); |
957 | sw.Print(); | |
7c3070ec | 958 | AliInfo("PreconILUKsymb<<"); |
8a9ab0eb | 959 | return 0; |
960 | } | |
961 | ||
962 | ||
963 | //___________________________________________________________ | |
964 | Int_t AliMinResSolve::PreconILUKsymbDense(Int_t lofM) | |
965 | { | |
966 | /*---------------------------------------------------------------------------- | |
967 | * ILUK preconditioner | |
968 | * incomplete LU factorization with level of fill dropping | |
969 | * Adapted from iluk.c: lofC of ITSOL_1 package by Y.Saad: http://www-users.cs.umn.edu/~saad/software/ | |
970 | *----------------------------------------------------------------------------*/ | |
971 | // | |
972 | UChar_t **ulvl=0,*levls=0; | |
973 | UShort_t *jbuf=0; | |
974 | Int_t *iw=0; | |
3c5b4cc8 | 975 | ulvl = new UChar_t*[fSize]; // stores lev-fils for U part of ILU factorization |
976 | levls = new UChar_t[fSize]; | |
977 | jbuf = new UShort_t[fSize]; | |
978 | iw = new Int_t[fSize]; | |
8a9ab0eb | 979 | // |
980 | for(int j=fSize; j--;) iw[j] = -1; // initialize iw | |
981 | for(int i=0; i<fSize; i++) { | |
982 | int incl = 0; | |
983 | int incu = i; | |
984 | // | |
985 | // assign lof = 0 for matrix elements | |
986 | for(int j=0;j<fSize; j++) { | |
7c185459 | 987 | if (AliMatrixSq::IsZero(fMatrix->Query(i,j))) continue; |
8a9ab0eb | 988 | if (j<i) { // L-part |
989 | jbuf[incl] = j; | |
990 | levls[incl] = 0; | |
991 | iw[j] = incl++; | |
992 | } | |
993 | else if (j>i) { // This works only for general matrix | |
994 | jbuf[incu] = j; | |
995 | levls[incu] = 0; | |
996 | iw[j] = incu++; | |
997 | } | |
998 | } | |
999 | // | |
1000 | // symbolic k,i,j Gaussian elimination | |
1001 | int jpiv = -1; | |
1002 | while (++jpiv < incl) { | |
1003 | int k = jbuf[jpiv] ; // select leftmost pivot | |
1004 | int kmin = k; | |
1005 | int jmin = jpiv; | |
1006 | for(int j=jpiv+1; j<incl; j++) if( jbuf[j]<kmin ) { kmin = jbuf[j]; jmin = j;} | |
1007 | // | |
1008 | // ------------------------------------ swap | |
1009 | if(jmin!=jpiv) { | |
1010 | jbuf[jpiv] = kmin; | |
1011 | jbuf[jmin] = k; | |
1012 | iw[kmin] = jpiv; | |
1013 | iw[k] = jmin; | |
1014 | int tj = levls[jpiv] ; | |
1015 | levls[jpiv] = levls[jmin]; | |
1016 | levls[jmin] = tj; | |
1017 | k = kmin; | |
1018 | } | |
1019 | // ------------------------------------ symbolic linear combinaiton of rows | |
1020 | AliVectorSparse& rowU = *fMatU->GetRow(k); | |
1021 | for(int j=0; j<rowU.GetNElems(); j++ ) { | |
1022 | int col = rowU.GetIndex(j); | |
1023 | int it = ulvl[k][j]+levls[jpiv]+1; | |
1024 | if( it > lofM ) continue; | |
1025 | int ip = iw[col]; | |
1026 | if( ip == -1 ) { | |
1027 | if( col < i) { | |
1028 | jbuf[incl] = col; | |
1029 | levls[incl] = it; | |
1030 | iw[col] = incl++; | |
1031 | } | |
1032 | else if( col > i ) { | |
1033 | jbuf[incu] = col; | |
1034 | levls[incu] = it; | |
1035 | iw[col] = incu++; | |
1036 | } | |
1037 | } | |
1038 | else levls[ip] = TMath::Min(levls[ip], it); | |
1039 | } | |
1040 | // | |
1041 | } // end - while loop | |
1042 | // | |
1043 | // reset iw | |
1044 | for (int j=0;j<incl;j++) iw[jbuf[j]] = -1; | |
1045 | for (int j=i;j<incu;j++) iw[jbuf[j]] = -1; | |
1046 | // | |
1047 | // copy L-part | |
1048 | AliVectorSparse& rowLi = *fMatL->GetRow(i); | |
1049 | rowLi.ReSize(incl); | |
1050 | if(incl>0) memcpy(rowLi.GetIndices(), jbuf, sizeof(UShort_t)*incl); | |
1051 | // copy U-part | |
1052 | int k = incu-i; | |
1053 | AliVectorSparse& rowUi = *fMatU->GetRow(i); | |
1054 | rowUi.ReSize(k); | |
1055 | if( k > 0 ) { | |
1056 | memcpy(rowUi.GetIndices(), jbuf+i, sizeof(UShort_t)*k); | |
1057 | ulvl[i] = new UChar_t[k]; // update matrix of levels | |
1058 | memcpy( ulvl[i], levls+i, k*sizeof(UChar_t) ); | |
1059 | } | |
1060 | } | |
1061 | // | |
1062 | // free temp space and leave | |
1063 | delete[] levls; | |
1064 | delete[] jbuf; | |
1065 | for(int i=fSize; i--;) if (fMatU->GetRow(i)->GetNElems()) delete[] ulvl[i]; | |
1066 | delete[] ulvl; | |
f8b1c575 | 1067 | delete[] iw; |
8a9ab0eb | 1068 | // |
1069 | fMatL->SortIndices(); | |
1070 | fMatU->SortIndices(); | |
1071 | return 0; | |
1072 | } |