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

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

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

revision 4855 by caltinay, Wed Apr 9 03:58:08 2014 UTC revision 4856 by caltinay, Wed Apr 9 06:46:55 2014 UTC
# Line 34  Line 34 
34  #include "Common.h"  #include "Common.h"
35  #include "Solver.h"  #include "Solver.h"
36  #include "PasoUtil.h"  #include "PasoUtil.h"
 #ifdef _OPENMP  
 #include <omp.h>  
 #endif  
 err_t Paso_Solver_NewtonGMRES(  
     Paso_Function* F,    /* function evaluation */  
     double *x,           /* in: initial guess, out: new approximation */  
     paso::Options* options,  
     Paso_Performance* pp)  
37    
38    namespace paso {
39    
40    err_t Solver_NewtonGMRES(Function* F, double* x, Options* options,
41                             Paso_Performance* pp)
42  {  {
43     const double inner_tolerance_safety=.9;      const double inner_tolerance_safety=.9;
44     dim_t gmres_iter;      dim_t gmres_iter;
45     double stop_tol, norm2_f,norm2_fo, normsup_f,reduction_f, gmres_tol, rtmp, quad_tolerance;      double stop_tol, norm2_f,norm2_fo, normsup_f,reduction_f, gmres_tol, rtmp, quad_tolerance;
46     bool convergeFlag=FALSE, maxIterFlag=FALSE, breakFlag=FALSE;      bool 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     const bool debug=options->verbose;      const bool debug=options->verbose;
50     const dim_t n = F->n;      const dim_t n = F->getLen();
51     dim_t iteration_count=0;      dim_t iteration_count=0;
52     const double atol=options->absolute_tolerance;  /* absolute tolerance */      const double atol=options->absolute_tolerance;  /* absolute tolerance */
53     const double rtol=options->tolerance;           /* relative tolerance */      const double rtol=options->tolerance;           /* relative tolerance */
54     const dim_t maxit=options->iter_max;            /* max iteration counter */      const dim_t maxit=options->iter_max;            /* max iteration counter */
55     const dim_t lmaxit=options->inner_iter_max*10;  /* max inner iteration counter */      const dim_t lmaxit=options->inner_iter_max*10;  /* max inner iteration counter */
56     const bool adapt_inner_tolerance=options->adapt_inner_tolerance;      const bool adapt_inner_tolerance=options->adapt_inner_tolerance;
57     const double max_inner_tolerance=options->inner_tolerance *1.e-11;        const double max_inner_tolerance=options->inner_tolerance *1.e-11;  
58     double inner_tolerance=max_inner_tolerance;      double inner_tolerance=max_inner_tolerance;
59    /*      /*
60     * max_inner_tolerance = 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 smaller than     *              when the relative linear residual is smaller than
# Line 68  err_t Paso_Solver_NewtonGMRES( Line 64  err_t Paso_Solver_NewtonGMRES(
64     *              by the modified Eisenstat-Walker formula, if     *              by the modified Eisenstat-Walker formula, if
65     *              adapt_inner_tolerance is set, otherwise     *              adapt_inner_tolerance is set, otherwise
66     *              inner_tolerance = max_inner_tolerance iteration.     *              inner_tolerance = max_inner_tolerance iteration.
67     */      */
68        
69    f=new double[n];      f=new double[n];
70    step=new double[n];      step=new double[n];
71    if (f==NULL || step ==NULL) {      /*
72        Status=SOLVER_MEMORY_ERROR;       * initial evaluation of F
73    } else {       */
74       /*      F->call(f, x, pp);
75        * initial evaluation of F      iteration_count++;
76        */      norm2_f=paso::util::l2(n,f,F->mpi_info);
77       Paso_FunctionCall(F,f,x,pp);      normsup_f=paso::util::lsup(n,f,F->mpi_info);
78       iteration_count++;      /*
79       norm2_f=paso::util::l2(n,f,F->mpi_info);       * stopping criterion:
80       normsup_f=paso::util::lsup(n,f,F->mpi_info);       */
81       /*      stop_tol=atol + rtol*normsup_f;
82        * stopping criterion:      if (stop_tol<=0) {
83        */          Status=SOLVER_INPUT_ERROR;
84       stop_tol=atol + rtol*normsup_f;          if (debug) printf("NewtonGMRES: zero tolerance given.\n");
85       if (stop_tol<=0) {      } else {
86         Status=SOLVER_INPUT_ERROR;          iteration_count=1;
87         if (debug) printf("NewtonGMRES: zero tolerance given.\n");          if (debug) {
88       } else {              printf("NewtonGMRES: Start Jacobi-free Newton scheme\n");
89         iteration_count=1;              printf("NewtonGMRES: lsup tolerance rel/abs= %e/%e\n",rtol,atol);
90         if (debug) {              printf("NewtonGMRES: lsup stopping tolerance = %e\n",stop_tol);
91          printf("NewtonGMRES: Start Jacobi-free Newton scheme\n");              printf("NewtonGMRES: max. inner iterations (GMRES) = %d\n",lmaxit);
         printf("NewtonGMRES: lsup tolerance rel/abs= %e/%e\n",rtol,atol);  
         printf("NewtonGMRES: lsup stopping tolerance = %e\n",stop_tol);  
           
         printf("NewtonGMRES: max. inner iterations (GMRES) = %d\n",lmaxit);  
         if (adapt_inner_tolerance) {  
             printf("NewtonGMRES: inner tolerance is adapted.\n");  
         printf("NewtonGMRES: max. inner l2 tolerance (GMRES) = %e\n",max_inner_tolerance);  
         } else {  
             printf("NewtonGMRES: inner l2 tolerance (GMRES) = %e\n",inner_tolerance);  
         }  
        }  
        /*  
         *  main iteration loop  
         */  
        while (! (convergeFlag || maxIterFlag || breakFlag)) {  
          /*  
           * keep track of the ratio (reduction_f = norm2_f/norm2_fo) of  
           * successive residual norms and the iteration counter (iteration_count)  
           */  
          if (debug) printf("NewtonGMRES: iteration step %d: lsup-norm of F =%e\n",iteration_count,normsup_f);  
          /*  
           * call GMRES to get increment  
           */  
          gmres_iter=lmaxit;  
          gmres_tol=inner_tolerance;  
          Status=Paso_Solver_GMRES2(F,f,x,step,&gmres_iter,&gmres_tol,pp);  
      inner_tolerance=MAX(inner_tolerance, gmres_tol/norm2_f);  
      printf("NewtonGMRES: actual rel. inner tolerance = %e\n",inner_tolerance);  
          iteration_count+=gmres_iter;  
          if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) {  
             Status=SOLVER_NO_ERROR;  
             /*  
              * update x:  
              */  
             norm2_fo=norm2_f;  
             paso::util::update(n,1.,x,1.,step);  
             Paso_FunctionCall(F,f,x,pp);  
             iteration_count++;  
             norm2_f=paso::util::l2(n,f,F->mpi_info);  
             normsup_f=paso::util::lsup(n,f,F->mpi_info);  
         reduction_f=norm2_f/norm2_fo;  
             /*  
              *   adjust inner_tolerance  
              */  
92              if (adapt_inner_tolerance) {              if (adapt_inner_tolerance) {
93           quad_tolerance = inner_tolerance_safety * reduction_f * reduction_f;                  printf("NewtonGMRES: inner tolerance is adapted.\n");
94           rtmp=inner_tolerance_safety * inner_tolerance * inner_tolerance;                  printf("NewtonGMRES: max. inner l2 tolerance (GMRES) = %e\n",max_inner_tolerance);
95                   if (rtmp>.1)  {              } else {
96                inner_tolerance=MIN(max_inner_tolerance, MAX(quad_tolerance,rtmp));                  printf("NewtonGMRES: inner l2 tolerance (GMRES) = %e\n",inner_tolerance);
97           } else {              }
98                inner_tolerance=MIN(max_inner_tolerance, quad_tolerance);          }
99           }          /*
100                   inner_tolerance=MIN(max_inner_tolerance, MAX(inner_tolerance, .5*stop_tol/normsup_f));           *  main iteration loop
101             */
102            while (! (convergeFlag || maxIterFlag || breakFlag)) {
103                // keep track of the ratio (reduction_f = norm2_f/norm2_fo) of
104                // successive residual norms and the iteration counter
105                // (iteration_count)
106                if (debug)
107                    printf("NewtonGMRES: iteration step %d: lsup-norm of F =%e\n",iteration_count,normsup_f);
108    
109                // call GMRES to get increment
110                gmres_iter=lmaxit;
111                gmres_tol=inner_tolerance;
112                Status = Solver_GMRES2(F,f,x,step,&gmres_iter,&gmres_tol,pp);
113                inner_tolerance=MAX(inner_tolerance, gmres_tol/norm2_f);
114                printf("NewtonGMRES: actual rel. inner tolerance = %e\n",inner_tolerance);
115                iteration_count+=gmres_iter;
116                if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) {
117                    Status=SOLVER_NO_ERROR;
118                    // update x
119                    norm2_fo=norm2_f;
120                    paso::util::update(n,1.,x,1.,step);
121                    F->call(f, x, pp);
122                    iteration_count++;
123                    norm2_f=paso::util::l2(n,f,F->mpi_info);
124                    normsup_f=paso::util::lsup(n,f,F->mpi_info);
125                    reduction_f=norm2_f/norm2_fo;
126                    // adjust inner_tolerance
127                    if (adapt_inner_tolerance) {
128                        quad_tolerance = inner_tolerance_safety * reduction_f * reduction_f;
129                        rtmp=inner_tolerance_safety * inner_tolerance * inner_tolerance;
130                        if (rtmp>.1)  {
131                            inner_tolerance=MIN(max_inner_tolerance, MAX(quad_tolerance,rtmp));
132                        } else {
133                            inner_tolerance=MIN(max_inner_tolerance, quad_tolerance);
134                        }
135                        inner_tolerance=MIN(max_inner_tolerance, MAX(inner_tolerance, .5*stop_tol/normsup_f));
136                    }
137                    convergeFlag = (normsup_f <= stop_tol);
138                } else {
139                    breakFlag = true;
140              }              }
141              convergeFlag = (normsup_f <= stop_tol);              maxIterFlag = (iteration_count > maxit);
142           } else {          }
143              breakFlag=TRUE;          if (debug) {
144           }              printf("NewtonGMRES: iteration step %d: lsup-norm of F =%e\n",iteration_count,normsup_f);
145           maxIterFlag = (iteration_count > maxit);  
146           }              if (convergeFlag) printf("NewtonGMRES: convergence reached after %d steps.\n",iteration_count);
147           if (debug) {              if (breakFlag) printf("NewtonGMRES: iteration break down after %d steps.\n",iteration_count);
148                printf("NewtonGMRES: iteration step %d: lsup-norm of F =%e\n",iteration_count,normsup_f);              if (maxIterFlag) printf("NewtonGMRES: maximum number of iteration steps %d reached.\n",maxit);
149            }
150                if (convergeFlag) printf("NewtonGMRES: convergence reached after %d steps.\n",iteration_count);          if (breakFlag) Status = SOLVER_BREAKDOWN;
151                if (breakFlag) printf("NewtonGMRES: iteration break down after %d steps.\n",iteration_count);          if (maxIterFlag) Status = SOLVER_MAXITER_REACHED;
152                if (maxIterFlag) printf("NewtonGMRES: maximum number of iteration steps %d reached.\n",maxit);      }
153           }      delete[] f;
154          if (breakFlag) Status=SOLVER_BREAKDOWN;      delete[] step;
155          if (maxIterFlag) Status=SOLVER_MAXITER_REACHED;      if (debug) printf("NewtonGMRES: STATUS return = %d\n",Status);
156        }      return Status;
   }  
   delete[] f;  
   delete[] step;  
   if (debug) printf("NewtonGMRES: STATUS return = %d\n",Status);  
   return Status;  
157  }  }
158    
159    } // namespace paso
160    

Legend:
Removed from v.4855  
changed lines
  Added in v.4856

  ViewVC Help
Powered by ViewVC 1.1.26