/[escript]/trunk/paso/src/NewtonGMRES.c
ViewVC logotype

Diff of /trunk/paso/src/NewtonGMRES.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1638 by ksteube, Mon Jul 14 05:34:59 2008 UTC revision 1639 by gross, Mon Jul 14 08:55:25 2008 UTC
# Line 29  Line 29 
29    
30  #include "Common.h"  #include "Common.h"
31  #include "Solver.h"  #include "Solver.h"
32  #include "Util.h"  #include "PasoUtil.h"
33  #ifdef _OPENMP  #ifdef _OPENMP
34  #include <omp.h>  #include <omp.h>
35  #endif  #endif
36  err_t Paso_Solver_NewtonGMRES(  err_t Paso_Solver_NewtonGMRES(
37      Paso_Function *F,      Paso_Function *F,    /* function evaluation */
38      double *x,      double *x,           /* in: initial guess, out: new approximation */
39      Paso_Options* options,      Paso_Options* options,
40      Paso_Performance* pp)      Paso_Performance* pp)
41    
42  {  {
43     const double gamma=.9;     const double inner_tolerance_safety=.9;
44     dim_t gmres_iter;     dim_t gmres_iter;
45     double stop_tol, fnrm,fnrmo, rat,etaold,etanew, gmres_tol;     double stop_tol, norm_f,norm_fo, reduction_f,old_inner_tolerance, gmres_tol, rtmp;
46     bool_t convergeFlag=FALSE, maxIterFlag=FALSE, breakFlag=FALSE;     bool_t convergeFlag=FALSE, maxIterFlag=FALSE, breakFlag=FALSE;
47     double *f=NULL, *step=NULL;     double *f=NULL, *step=NULL;
48     err_t Status=SOLVER_NO_ERROR;     err_t Status=SOLVER_NO_ERROR;
49     bool_t debug=options->verbose;     const bool_t debug=options->verbose;
50     dim_t n = n=F->local_n;     const dim_t n = F->n;
51     dim_t itc=0;     dim_t iteration_count=0;
52     double atol=options->absolute_tolerance;     const double atol=options->absolute_tolerance;  /* absolute tolerance */
53     double rtol=options->tolerance;     const double rtol=options->tolerance;           /* relative tolerance */
54     dim_t maxit=options->iter_max;     const dim_t maxit=options->iter_max;            /* max iteration counter */
55     dim_t lmaxit=options->inner_iter_max;     const dim_t lmaxit=options->inner_iter_max;     /* max inner iteration counter */
56     bool_t adapt_eta=options->adapt_inner_tolerance;     const bool_t adapt_inner_tolerance=options->adapt_inner_tolerance;
57     double etamax=options->inner_tolerance;     const double max_inner_tolerance=options->inner_tolerance;  /*inner tolerance counter */
58     double eta=etamax;     double inner_tolerance=max_inner_tolerance;
59    /*    /*
60     * etamax = Maximum error tolerance for residual in inner     * max_inner_tolerance = Maximum error tolerance for residual in inner
61     *              iteration. The inner iteration terminates     *              iteration. The inner iteration terminates
62     *              when the relative linear residual is     *              when the relative linear residual is smaller than inner_tolerance*| F(x_c) |. inner_tolerance is determined
63     *              smaller than eta*| F(x_c) |. eta is determined     *              by the modified Eisenstat-Walker formula, if adapt_inner_tolerance is set, otherwise
64     *              by the modified Eisenstat-Walker formula if etamax > 0.     *              inner_tolerance = max_inner_tolerance iteration.
    *              If etamax < 0, then eta = |etamax| for the entire  
    *              iteration.  
65     */     */
66        
67    f=TMPMEMALLOC(n,double);    f=TMPMEMALLOC(n,double);
# Line 75  err_t Paso_Solver_NewtonGMRES( Line 73  err_t Paso_Solver_NewtonGMRES(
73        * initial evaluation of F        * initial evaluation of F
74        */        */
75       Paso_FunctionCall(F,f,x);       Paso_FunctionCall(F,f,x);
76       fnrm=Paso_l2(n,f,F->mpi_info);       norm_f=Paso_l2(n,f,F->mpi_info);
77       /*       /*
78        * stoping criterium:        * stoping criterium:
79        */        */
80       stop_tol=atol + rtol*fnrm;       stop_tol=atol + rtol*norm_f;
81       itc=0;       iteration_count=0;
82       if (debug) printf("Start Jacobi-free Newton scheme:\n\ttolerance = %e\n\tstoping tolerance = %e\n",rtol,stop_tol);       if (debug) printf("Start Jacobi-free Newton scheme:\n\ttolerance = %e\n\tstoping tolerance = %e\n",rtol,stop_tol);
83       /*       /*
84        *  main iteration loop        *  main iteration loop
85        */        */
86       while (! (convergeFlag || maxIterFlag || breakFlag)) {       while (! (convergeFlag || maxIterFlag || breakFlag)) {
87           /*           /*
88            * keep track of the ratio (rat = fnrm/frnmo) of successive residual norms and            * keep track of the ratio (reduction_f = norm_f/frnmo) of successive residual norms and
89            * the iteration counter (itc)            * the iteration counter (iteration_count)
90            */            */
91           itc++;           iteration_count++;
92           if (debug) printf("iteration step %d: norm of F =%e\n",itc,fnrm);           if (debug) printf("iteration step %d: norm of F =%d\n",iteration_count,norm_f);
93           /*           /*
94            * call GMRES to get increment            * call GMRES to get increment
95            */            */
96           gmres_iter=lmaxit;           gmres_iter=lmaxit;
97           gmres_tol=eta;           gmres_tol=inner_tolerance;
98           if (debug) printf("GMRES called with tolerance = %e\n",eta);           if (debug) printf("GMRES called with tolerance = %e (max iter=%d)\n",inner_tolerance,gmres_iter);
99           Status=Paso_Solver_NLGMRES(F,f,x,step,&gmres_iter,&gmres_tol,pp);           Status=Paso_Solver_GMRES2(F,f,x,step,&gmres_iter,&gmres_tol,pp);
100           itc+=gmres_iter;           if (debug) printf("GMRES finalized after %d steps (residual = %e)\n",gmres_iter,gmres_tol);
101             iteration_count+=gmres_iter;
102           if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) {           if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) {
103              Status=SOLVER_NO_ERROR;              Status=SOLVER_NO_ERROR;
104              /*              /*
105               * update x:               * update x:
106               */               */
107                norm_fo=norm_f;
108              Paso_Update(n,1.,x,1.,step);              Paso_Update(n,1.,x,1.,step);
109              Paso_FunctionCall(F,f,x);              Paso_FunctionCall(F,f,x);
110              fnrm=Paso_l2(n,f,F->mpi_info);              iteration_count++;
111              fnrmo=fnrm;              norm_f=Paso_l2(n,f,F->mpi_info);
112              rat=fnrm/fnrmo;              reduction_f=norm_f/norm_fo;
113              /*              /*
114               *   adjust eta               *   adjust inner_tolerance
115               */               */
116               if (adapt_eta) {               if (adapt_inner_tolerance) {
117                   etaold=eta;                   old_inner_tolerance=inner_tolerance;
118                   etanew=gamma*rat*rat;                   inner_tolerance=inner_tolerance_safety * reduction_f * reduction_f;
119                   if (gamma*etaold*etaold >1.) etanew=MAX(etanew,gamma*etaold*etaold);                   rtmp=inner_tolerance_safety * old_inner_tolerance * old_inner_tolerance;
120                   eta=MAX(MIN(etanew,etamax), .5*stop_tol/fnrm);                   if (rtmp>.1) inner_tolerance=MAX(inner_tolerance,rtmp);
121                  }                   inner_tolerance=MAX(MIN(inner_tolerance,max_inner_tolerance), .5*stop_tol/norm_f);
122                  }
123                  convergeFlag = (norm_f <= stop_tol);
124           } else {           } else {
125             breakFlag=TRUE;             breakFlag=TRUE;
126           }           }
127           maxIterFlag = (itc > maxit);           maxIterFlag = (iteration_count > maxit);
          convergeFlag = (fnrm <= stop_tol);  
128        }        }
129    }    }
130    TMPMEMFREE(f);    TMPMEMFREE(f);

Legend:
Removed from v.1638  
changed lines
  Added in v.1639

  ViewVC Help
Powered by ViewVC 1.1.26