Solver.hpp

00001 //
00002 // Solvers
00003 //
00004 
00005 #ifndef _SOLVERS_
00006 # define _SOLVERS_
00007 
00008 # include "Operators.hpp"
00009 # include "SolverResidualCriteria.hpp"
00010 
00017 template<class I , int DIM>
00018 struct Solver {
00020             Solver() ;
00022   int  BiCGSTAB2(I *X , I *RHS) ;
00023   int  BiCGSTAB (I *X , I *RHS) ;
00024   int  CG       (I *X , I *RHS) ;
00025   int  GMRES    (I *X , I *RHS , const int m) ; 
00026 
00027   enum {NUMTMP_          = 100} ;
00029   enum {NUMTMP_BICGSTAB2 = 11 } ;
00031   enum {NUMTMP_CG        = 4  } ;
00032 
00034   Operator<I,DIM> *Op ;
00035 
00040   bool    (*StopCriterion )(double currentresidual , double startresidual , double stopvalue) ;
00041 
00043   double    eps     ;
00045   int       maxiter ;
00047   int       restart ;
00049   char     *print   ;
00051   int       prstep  ;
00053   I        *Tmp[NUMTMP_] ;                    //   , default NULL
00054   //I        *Residual ;
00055 
00057   double startresidual  ;
00059   double  finalresidual ;
00061   int     iterations    ;
00062   
00063   int     debugmode ;
00064 
00065   // for reference output of error |current - reference solution|_0
00066   I *RS;  // reference solution, if NULL no error is calculated
00067   I *tmp; // aux. memory
00068 } ;
00069 
00071 //                         //
00072 // Solver: Implementations //
00073 //                         //
00075 
00076 
00077 template<class I  , int DIM>
00078 Solver<I,DIM>::Solver() {
00079  eps    =1e-15  ;
00080  maxiter=1000   ;
00081  restart=MAXINT ;
00082  print  =NULL   ;
00083  prstep =10     ;
00084 
00085  Op=NULL  ;
00086  StopCriterion=NULL ;
00087 
00088  for (int i=0;i<NUMTMP_; i++) Tmp[i]=NULL ;
00089 
00090  debugmode=0 ;
00091  RS=tmp=NULL ;
00092 }
00093 
00094 //
00095 // preconditioned CG
00096 //
00097 
00098 template<class I  , int DIM>
00099 int Solver<I,DIM>::CG(I *X , I *RHS) 
00100 {
00101  int    i,it=0,lit ;
00102  
00103  double rho,rho2,alpha,beta , rr , ss, startres ;
00104 
00105  I *p,*q,*r,*z ;
00106 
00107  //Distribution and resize of aux. vectors
00108 
00109  for (i=0;i<NUMTMP_CG;i++) {
00110    Op->NotUsing(Tmp[i]) ;
00111 
00112    assert(Tmp[i]) ;
00113    assert(Tmp[i]!=X) ;
00114    assert(Tmp[i]!=RHS) ;
00115  }
00116 
00117  Op->NotUsing(X) ;
00118  Op->NotUsing(RHS) ;
00119  
00120  Op->InversePreconditioner(X,X) ;
00121 
00122  p =Tmp[0] ;
00123  q =Tmp[1] ;
00124  r =Tmp[2] ;
00125  z =Tmp[3] ;
00126 
00127  //Residual=r ;
00128 
00129  r->SetBoundaryConditions(RHS) ;
00130  z->SetBoundaryConditions(RHS) ;
00131 
00132  p->SetBoundaryConditions(X) ;
00133  q->SetBoundaryConditions(X) ; 
00134 
00135  // calc. start residual
00136 
00137 start:;
00138 
00139  for (i=0;i<NUMTMP_CG;i++) Tmp[i]->Set(0.0) ;
00140 
00141  lit=0 ;
00142 
00143  Op->Preconditioner(X,X) ;
00144  Op->Apply(X,r) ;
00145  r->Sub(RHS,r)  ;
00146  
00147  startres=1.0 ; 
00148 
00149  rr=r->InnerProd(r) ;
00150  if (print) std::cout<<print<<" initial residual "<<scientific<<rr<<"\n";
00151 
00152  if (StopCriterion==NULL) {std::cout<<"Solver: you must set StopCriterion \n"; exit(-1) ;} 
00153 
00154  while ( !StopCriterion(rr,startres,eps) && (it<maxiter)) {
00155    
00156    Op->Preconditioner(r,z) ;
00157    rho=r->InnerProd(z) ;
00158 
00159    if (it==0) p->Copy(z) ;
00160    else {
00161      beta=rho/rho2 ;
00162      p->Add(1.0,z , beta , p) ; 
00163    }
00164 
00165    Op->Apply(p,q) ;
00166    
00167    alpha=rho/p->InnerProd(q) ;
00168    
00169    X->PPlus(alpha,p) ;
00170 
00171    r->PPlus(-alpha,q) ;
00172 
00173    rr=r->InnerProd(r) ;
00174 
00175    rho2=rho ;
00176    it++ ;
00177 
00178    // Output
00179    if (print && (it%prstep==0)) {
00180      ss=r->MaxAbs() ;
00181      std::cout<<print<<' '<<setw(3)<<it<<' '<<scientific<<rr<<' '<<scientific<<ss<<"0.0"<<scientific<<fabs(rho)<<'\n' ;
00182    }
00183 
00184  }
00185 
00186  if (print) {
00187    ss=r->Max() ;
00188    std::cout<< "final  "<<it<<setw(3)<<' '<<scientific<<rr<<' '<<scientific<<ss<<'\n';
00189   }
00190 
00191  return it ;
00192 }
00193 
00194 //extern UniformData<2> K1,K2,K3,K4 ;
00195 
00196 //
00197 // preconditioned BiCGStab-Solver:
00198 //
00199 
00200 template<class I  , int DIM>
00201 int Solver<I,DIM>::BiCGSTAB(I *X ,I *RHS) 
00202 {
00203  int    i,it,lit ;
00204  
00205  double rhoi,rho1,alpha,beta,omega ;
00206  double rr,rr0,r0v,ts,tt,ss,startres,bestres,bestrss ;
00207 
00208  I *r,*p,*v,*s,*t,*r0,*y,*z,*x0,*best ;
00209 
00210  if (!X->Ext.IsSame(&RHS->Ext)) {
00211      X->Ext.Print() ;
00212    RHS->Ext.Print() ;
00213    exit(-1) ;
00214  }
00215 
00216  //Distribution and resize of aux. vectors
00217  for (i=0;i<NUMTMP_BICGSTAB2;i++) {
00218    Op->NotUsing(Tmp[i]) ;
00219    assert(Tmp[i]) ;
00220    assert(Tmp[i]!=X) ;
00221    assert(Tmp[i]!=RHS) ;
00222  }
00223 
00224  Op->NotUsing(X)   ;
00225  Op->NotUsing(RHS) ;
00226  for (i=0; i<=9; i++) Op->NotUsing(Tmp[i]) ;
00227 
00228  r =Tmp[0] ;
00229  p =Tmp[1] ;
00230  v =Tmp[2] ;
00231  s =Tmp[3] ;
00232  t =Tmp[4] ;
00233  r0=Tmp[5] ;
00234  y =Tmp[6] ;
00235  z =Tmp[7] ;
00236  x0=Tmp[8] ; 
00237  best=Tmp[9] ;
00238 
00239  // Residual=r ;
00240 
00241  r->CopyExtensions(RHS) ;
00242 r0->CopyExtensions(RHS) ;
00243  p->CopyExtensions(RHS) ;
00244  s->CopyExtensions(RHS) ;
00245  t->CopyExtensions(RHS) ;
00246  v->CopyExtensions(RHS) ;
00247 
00248  y->CopyExtensions(X) ;
00249  z->CopyExtensions(X) ;
00250 x0->CopyExtensions(X) ;
00251 
00252  it=0 ;
00253 start:;
00254 
00255  lit=0 ;
00256  Op->Apply(X,r) ;
00257  r ->Sub(RHS,r) ;
00258  r0->Copy(r)   ;
00259 
00260  // initialize with awful values
00261  rho1  = omega = 1e-200 ;
00262  alpha = bestres = bestrss = startres = 1e+100 ;
00263  
00264  rr = r->InnerProd(r) ;
00265  rr0=rr ;
00266  if (it==0) startres=rr ;
00267 
00268  ss=r->MaxAbs() ;
00269  if (print) std::cout<<print<<" initial residual "<<scientific<<rr<<' '<<scientific<<ss<<'\n' ;
00270 
00271  if (StopCriterion==NULL) {std::cout<<"Solver: you must set StopCriterion \n"; exit(-1) ;} 
00272 
00273  while ( !StopCriterion(rr,startres,eps) && (it<maxiter)) {
00274    
00275    rhoi=r0->InnerProd(r) ;
00276    
00277    if ( fabs(rhoi)<1e-15 ) {
00278      if (print) std::cout<<print<<" restart: "<<scientific<<rr/rr0<<" Forced:"<<(lit>restart)<<'\n' ;
00279      X->Copy(best) ;
00280      goto start ;
00281    }
00282 
00283    if (print && (it%prstep == 0)) std::cout<<print<<' '<<setw(3)<<it<<' '<<scientific<<rr<<' '<<sscientific<<s<<' '<<scientific<<fabs(rhoi)<<'\n' ;
00284    
00285    if (lit==0) { 
00286      p->Copy(r) ;
00287    } else {
00288      beta=(rhoi/rho1)*(alpha/omega) ;
00289      p->Add(1.0,r , beta,p , -beta*omega,v) ;
00290    }
00291 
00292    rho1=rhoi ;
00293    Op->Preconditioner(p,y) ;
00294    Op->Apply(y,v) ;
00295  
00296    r0v=r0->InnerProd(v) ;
00297    alpha=(rhoi + 1e-15)/(r0v + 1e-15) ;
00298    s->Add  (1.0,r , -alpha,v) ;
00299 
00300    ss = s->InnerProd(s) ;
00301    if (StopCriterion(ss,startres,eps)) {
00302      X->PPlus(alpha,y) ;
00303      rr=ss ;
00304      ss=s->MaxAbs() ; 
00305      break ;
00306    }
00307 
00308    Op->Preconditioner(s,z) ;
00309    Op->Apply(z,t) ;
00310    
00311    tt=t->InnerProd(t) ;
00312    ts=t->InnerProd(s) ;
00313    omega=(ts + 1e-15)/(tt + 1e-15); 
00314 
00315    X->Add(1.0,X, alpha,y, omega,z) ;
00316    r->Add(1.0,s , -omega,t) ;
00317 
00318    rr = r->InnerProd(r)  ;
00319    ss = r->MaxAbs() ;
00320    if (rr<bestres) {
00321      best->Copy(X)   ;
00322      bestres=rr ;
00323      bestrss=ss ;
00324    }
00325 
00326    it++;
00327    lit++ ;
00328  }
00329 
00330  // get best previous solution if necessary
00331  if (rr>bestres) { 
00332    X->Copy(best)   ; 
00333    rr=bestres ;
00334    ss=bestrss ;
00335   }
00336 
00337  // exact evaluation of last residual
00338   Op->Apply(X,r) ;
00339   r ->Sub(RHS,r) ;
00340  rr = r->InnerProd(r) ;
00341  ss = r->MaxAbs() ;
00342 
00343  if (print) {
00344    if (X->Ext.dimext) std::cout<<"It: "<<setw(3)<<it<<" Max: "<<scientific<<ss<<" L2: "<<scientific<<rr<<" defect: "<<scientific<<X->Ext.ext[0]<<'\n' ;
00345                  else std::cout<<"It: "<<setw(3)<<it<<" Max: "<<scientific<<ss<<" L2: "<<scientific<<rr<<"\n" ;
00346   } 
00347 
00348  for (i=0;i<NUMTMP_BICGSTAB2;i++) Tmp[i]->Ext.dimext=0 ;
00349 
00350  return it ;
00351 }
00352 
00353 //
00354 // Preconditioned BiCGStab(2)
00355 //
00356 
00357 template<class I  , int DIM>
00358 int Solver<I,DIM>::BiCGSTAB2(I *X , I *RHS) 
00359 {
00360  startresidual=finalresidual=-1 ;
00361  iterations=0 ;
00362 
00363  int    i,lit ;
00364  
00365  double rho0,rho1, omega1,omega2 , alpha,beta,gamma , mu,nu,tau;
00366  double rr,ss,bestresidual=1e+100,realerror=1e+100 ;
00367 
00368  I *r,*q,*v,*s,*t,*r0,*u,*y,*x0,*w,*best ;
00369 
00370  if (!X->Ext.IsSame(&RHS->Ext)) {
00371      X->Ext.Print() ;
00372    RHS->Ext.Print() ;
00373    exit(-1) ;
00374  }
00375  //Distribution and resize of aux. vectors
00376  
00377  for (i=0;i<NUMTMP_BICGSTAB2;i++) {
00378    Op->NotUsing(Tmp[i]) ;
00379 
00380    assert(Tmp[i]) ;
00381    assert(Tmp[i]!=X) ;
00382    assert(Tmp[i]!=RHS) ;
00383   
00384    //Tmp[i]->Ext.dimext=X->Ext.dimext ;
00385  }
00386 
00387  Op->NotUsing(X) ;
00388  Op->NotUsing(RHS) ;
00389  
00390  Op->InversePreconditioner(X,X) ;
00391 
00392  r =Tmp[0] ;
00393  r0=Tmp[1] ;
00394  u =Tmp[2] ;
00395  v =Tmp[3] ;
00396  q =Tmp[4] ;
00397  s =Tmp[5] ;
00398  t =Tmp[6] ;
00399  y =Tmp[7] ;
00400  w =Tmp[8] ;
00401 
00402  x0  =Tmp[9] ;
00403  best=Tmp[10] ;
00404 
00405  //Residual=r ;
00406 
00407  // Set Flags
00408 
00409  r->CopyExtensions(RHS) ;
00410 r0->CopyExtensions(RHS) ;
00411  u->CopyExtensions(RHS) ;
00412  v->CopyExtensions(RHS) ;
00413  q->CopyExtensions(RHS) ;
00414  s->CopyExtensions(RHS) ;
00415 
00416  t->CopyExtensions(X) ;
00417  y->CopyExtensions(X) ;
00418  w->CopyExtensions(X) ;
00419 
00420   x0->CopyExtensions(X) ;
00421 best->CopyExtensions(X) ;
00422 
00423  // calc. start residual
00424 
00425 start:;
00426 
00427  for (i=0;i<NUMTMP_BICGSTAB2;i++) Tmp[i]->Set(0.0) ;
00428 
00429  lit=0 ;
00430 
00431  Op->Preconditioner(X,x0) ;
00432  Op->Apply(x0,r) ;
00433  r->Sub(RHS,r)   ;
00434  r0->Copy(r)     ;
00435  
00436  rho0=omega2=1.0   ;
00437  alpha=0.0   ;
00438  u->Set(0.0) ;
00439 
00440  startresidual=rr=r->InnerProd(r) ;
00441  if (print) std::cout<<print<<" start residual "<<scientific<<startresidual<<"\n" ;
00442  best->Copy(X) ;
00443  bestresidual=rr ;
00444 
00445  if (StopCriterion==NULL) {std::cout<<"Solver: you must set StopCriterion \n"; exit(-1) ;} 
00446 
00447  while ( !StopCriterion(rr,startresidual,eps) && (iterations < maxiter)) {
00448    
00449    rho0 = -omega2*rho0 ;
00450    rho1=r0->InnerProd(r) ;
00451 
00452    // Output
00453    if (print && (iterations%prstep==0)) {
00454      ss=r->MaxAbs() ;
00455 
00456      if (RS) {
00457        tmp->CopyExtensions(X)     ;
00458        Op ->Preconditioner(X,tmp) ;
00459        tmp->MMinus(RS)            ;
00460        realerror=tmp->L2Norm(tmp) ;
00461      }
00462 
00463      std::cout<<print<<' '<<setw(3)<<iterations<<' '<<scientific<<rr<<' '<<scientific<<ss<<' '<<scientific<<realerror<<' '<<scientific<<fabs(rho1)<<'\n' ; 
00464    }
00465 
00466    if (lit>restart) {
00467      std::cout<<"Forced restart\n" ;
00468      X->Copy(x0) ;
00469      goto start ;
00470    }
00471 
00472 
00473    beta=alpha*rho1/rho0 ;
00474    rho0=rho1 ;
00475    
00476    u->Add(1.0,r, -beta,u) ;
00477 
00478    Op->Preconditioner(u,x0) ;
00479    Op->Apply(x0,v) ;
00480 
00481    gamma=r0->InnerProd(v) ;
00482    alpha=rho0/gamma;
00483 
00484    q->Add(1.0,r , -alpha,v) ;
00485 
00486    Op->Preconditioner(q,x0) ;
00487    Op->Apply(x0,s) ;
00488 
00489    y->Add(1.0,X, alpha,u) ;
00490    
00491 // odd BiCG
00492    rho1=r0->InnerProd(s) ;
00493    beta=alpha*rho1/rho0 ;
00494    rho0=rho1 ;
00495 
00496    v->Add(1.0,s , -beta,v) ;
00497 
00498    Op->Preconditioner(v,x0) ;
00499    Op->Apply(x0,w) ;
00500 
00501    gamma=r0->InnerProd(w) ;
00502    alpha=rho0/gamma ;
00503 
00504    u->Add  (1.0,q , -beta ,u) ;
00505    q->PPlus(        -alpha,v) ;
00506    s->PPlus(        -alpha,w) ;
00507     
00508    Op->Preconditioner(s,x0) ;
00509    Op->Apply(x0,t) ;
00510 
00511 //GCR(2)
00512    omega1=q->InnerProd(s) ;
00513    mu    =s->InnerProd(s) ;
00514    nu    =t->InnerProd(s) ;
00515    tau   =t->InnerProd(t) ;
00516    omega2=q->InnerProd(t) ;
00517    tau   =tau-nu*nu/mu ;
00518    omega2=(omega2-nu*omega1/mu)/tau ;
00519    omega1=(omega1-nu*omega2)/mu ;
00520 
00521    x0->Copy(X) ;
00522  
00523    X->Add(1.0,y ,  omega1,q ,  omega2,s , alpha,u) ;
00524    r->Add(1.0,q , -omega1,s , -omega2,t ) ;
00525 
00526    rr=r->InnerProd(r) ;
00527 
00528    if (rr<bestresidual) {
00529      best->Copy(X)   ;
00530      bestresidual=rr ;
00531    }
00532 
00533    u->Add(1.0,u , -omega1,v , -omega2,w ) ;
00534    
00535    iterations++;
00536    lit++ ;
00537  }
00538 
00539  // get best previous solution if necessary
00540  if (rr > bestresidual) { X->Copy(best) ; rr=bestresidual ; }
00541  finalresidual=rr ;
00542  
00543  Op->Preconditioner(X,X) ;
00544 
00545  if (print) {
00546    ss=r->Max() ;
00547    if (X->Ext.dimext) std::cout<<"It: "<<setw(3)<<iterations<<" Max: "<<scientific<<ss<<" L2: "<<scientific<<rr<<" defect: "<<scientific<<X->Ext.ext[0]<<'\n' ;
00548                  else std::cout<<"It: "<<setw(3)<<iterations<<" Max: "<<scientific<<ss<<" L2: "<<scientific<<rr<<"\n" ;
00549   }
00550 
00551  for (i=0;i<NUMTMP_BICGSTAB2;i++) Tmp[i]->Ext.dimext=0 ;
00552 
00553  return iterations ;
00554 }
00555 
00556 
00557 
00558 //
00559 // GMRES 
00560 //
00561 template <class I, int DIM>
00562 int Solver<I,DIM>::GMRES(I *X , I *RHS , const int m) {
00563  int i,k,l,j,it=0 ;
00564  // Checks
00565  assert(m < min(NUMTMP_ , 64-1) ) ;
00566  for (l=0; l<m+1; l++) assert (Tmp[l]) ;
00567  
00568  // Allocation
00569  double h[64-1][64] , y[64-1] , s[64] , c1[64-1] , c2[64-1] ,error, startres ;
00570  I *v[64],*z ;
00571  
00572  z=Tmp[0] ;
00573  z->SetBoundaryConditions(X->BC) ;
00574  for(l=0; l<m; l++) { 
00575    v[l]=Tmp[l+1] ;
00576    v[l]->SetBoundaryConditions(RHS->BC) ;
00577  }
00578  
00579  // Initial solution/residual
00580   
00581  Op->Apply(X,v[0]) ;
00582  v[0]->Add(1.0,RHS , -1.0,v[0]) ;
00583 
00584  error=v[0]->InnerProd(v[0]) ;
00585  startres=1.0 ; 
00586 
00587  if (StopCriterion==NULL) {std::cout<<"Solver: you must set StopCriterion \n"; exit(-1) ;} 
00588 
00589  if (StopCriterion(error,startres,eps)) return 0 ;
00590 
00591  // really iterate
00592  startres=error ;
00593  Op->InversePreconditioner(X,X) ;
00594  
00595  //
00596  // outer loop
00597  while ((!StopCriterion(error,startres,eps)) && (it<maxiter)) {
00598    double vv;
00599    vv=v[0]->InnerProd(v[0]) ;
00600    s[0]=vv ;
00601    v[0]->Mul(1/vv) ;
00602    
00603    // inner loop
00604    i=0 ;
00605    while ( (fabs(s[i])>eps) && (i<m-1) && (it<maxiter)) {
00606      I *w=v[i+1] ;
00607      Op->Preconditioner(v[i] , z) ;
00608      Op->Apply(z,w) ;
00609      
00610      // Gram Schmidt
00611      for (k=0; k<i; k++) {
00612        h[i][k]=w->InnerProd(v[k]) ;
00613        w->Add(1.0,w , -h[i][k],v[k]) ;
00614      }
00615      h[i][i+1]=w->InnerProd(w) ;
00616      w->Mul(1/h[i][i+1]) ; // v[i+1]==w !!
00617 
00618      // QR for Hessenberg matrix
00619      for (k=0; k<i; k++) {
00620        double h1 =  c1[k]*h[i][k] + c2[k]*h[i][k+1] ;
00621        double h2 = -c2[k]*h[i][k] + c1[k]*h[i][k+1] ;
00622     
00623        h[i][k  ]=h1 ;
00624        h[i][k+1]=h2 ;
00625      }
00626       
00627      double r=sqrt(h[i][i]*h[i][i] + h[i][i+1]*h[i][i+1]) ;
00628      c1[i]    =h[i][i  ]/r ;
00629      c2[i]    =h[i][i+1]/r ;
00630      h[i][i  ]=r ;
00631      h[i][i+1]=0 ;
00632      s[i+1]   =-c2[i]*s[i] ;
00633      s[i  ]   = c1[i]*s[i] ;
00634      
00635      error=fabs(s[i]) ;
00636 
00637      if (print && (it%prstep==0)) std::cout<<print<<' '<<it<<' '<<error<<'\n' ;
00638      i++  ;
00639      it++ ;
00640    }
00641 
00642    //Update:  Solve linear system 
00643    i=i-1 ;
00644    for (j=i; j>=0; j--) {
00645      y[j]=s[j]/h[j][j] ;
00646      for (k=j-1; k>=0; k--) s[k] -= h[j][k]*y[j] ;
00647    }
00648 
00649    for (j=i; j>=0; j--) X->Add(1.0,X ,  y[j],v[j]) ;
00650 
00651    Op->Preconditioner(X , z) ;
00652    Op->Apply(z,v[0]) ;
00653    v[0]->Add(1.0,RHS , -1.0,v[0]) ;
00654    if (print) std::cout<< "GMRES restart\n" ;
00655  }
00656  
00657  if (print) std::cout<< "final "<<error<<"\n" ; 
00658  
00659  Op->Preconditioner(X , X) ;
00660  return it ;
00661 }
00662 
00663 #endif
Generated at Mon Aug 19 10:02:32 2002 for AWFD by