54d47f07de5ab16bc6209165d6a82873e5554b06
[u/mrichter/AliRoot.git] / STEER / AliMinResSolve.cxx
1 #include "AliMinResSolve.h"
2 #include "AliLog.h"
3 #include <TStopwatch.h>
4
5 using namespace std;
6
7 ClassImp(AliMinResSolve)
8
9 //______________________________________________________
10 AliMinResSolve::AliMinResSolve() : 
11 fSize(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 //______________________________________________________
18 AliMinResSolve::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 //______________________________________________________
27 AliMinResSolve::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 //______________________________________________________
35 AliMinResSolve::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 //______________________________________________________
43 AliMinResSolve::~AliMinResSolve()
44 {
45   ClearAux();
46 }
47
48 //______________________________________________________
49 AliMinResSolve& 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 //_______________________________________________________________
61 Int_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
73     for (int i=0;i<fSize;i++) fPrecDiag[i] = fMatrix->QueryDiag(i);
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++) {
78         double vl = fMatrix->Query(j,i);
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 ________________________________
95 Bool_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 //________________________________________________________________________________
101 Bool_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 ________________________________
256 Bool_t AliMinResSolve::SolveMinRes(TVectorD& VecSol,Int_t precon,int itnlim,double rtol)
257 {
258   return SolveMinRes(VecSol.GetMatrixArray(), precon, itnlim, rtol);
259 }
260
261 //________________________________________________________________________________
262 Bool_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   //
314   for (int i=fSize;i--;) fPVecY[i]  = fPVecR1[i] = fRHS[i];
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 //______________________________________________________________
498 void AliMinResSolve::ApplyPrecon(const TVectorD& vecRHS, TVectorD& vecOut) const
499 {
500   ApplyPrecon(vecRHS.GetMatrixArray(), vecOut.GetMatrixArray());
501 }
502
503 //______________________________________________________________
504 void 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--;) {
511       double d = fMatrix->QueryDiag(i);
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;
520       for (int j=i;j--;) {double vl = fMatrix->Query(i,j);if (vl!=0) el += fPrecAux[j]*vl;}
521       fPrecAux[i] = fPrecDiag[i]*(vecRHS[i]-el); 
522     }
523     for (int i=fSize;i--;) {
524       double el = 0;
525       for (int j=i+1;j<fSize;j++) {double vl = fMatrix->Query(i,j);if (vl!=0) el += vl*vecOut[j];}
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 //___________________________________________________________
554 Bool_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 //___________________________________________________________
578 Bool_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 //___________________________________________________________
605 void 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 //___________________________________________________________
621 Int_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
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     }
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 
682         double vl = fMatrix->Query(col,i);    // the lower part of the column I
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 //___________________________________________________________
727 Int_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 
773       double vl = fMatrix->Query(i,j);
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 //___________________________________________________________
819 Int_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   //
827   TStopwatch sw;
828   printf("PreconILUKsymb>>\n");
829   sw.Start();
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
874         if (fMatrix->Query(col,i)==0) continue;    // Due to the symmetry  == matrix(i,col)
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();
950   sw.Stop();
951   sw.Print();
952   printf("PreconILUKsymb<<\n");
953   return 0;
954 }
955
956
957 //___________________________________________________________
958 Int_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++) {
993       if (fMatrix->Query(i,j)==0) continue;
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 }