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