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