paso/src/NewtonGMRES.c

/* $Id:$ */

/*******************************************************
 *
 *       Copyright 2008 by University of Queensland
 *
 *                http://esscc.uq.edu.au
 *        Primary Business: Queensland, Australia
 *  Licensed under the Open Software License version 3.0
 *     http://www.opensource.org/licenses/osl-3.0.php
 *
 *******************************************************/

/*
*  Purpose
*  =======
*
*  NewtonGMRES solves the non-linear system F(x)=0 using the
*  restarted GMRES method.
*
*  Convergence test: norm(F(x))< rtol*norm(F(x0))+atol
*
*  x input initial guess
*    output solution approximation
*
*  based on code by C. T. Kelley, July 1, 1994.
*
*/

#include "Common.h"
#include "Solver.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) 

{
   const double inner_tolerance_safety=.9;
   dim_t gmres_iter;
   double stop_tol, norm_f,norm_fo, reduction_f,old_inner_tolerance, gmres_tol, rtmp;
   bool_t convergeFlag=FALSE, maxIterFlag=FALSE, breakFlag=FALSE;
   double *f=NULL, *step=NULL;
   err_t Status=SOLVER_NO_ERROR;
   const bool_t debug=options->verbose;
   const dim_t n = F->n; 
   dim_t iteration_count=0;
   const double atol=options->absolute_tolerance;  /* absolute tolerance */
   const double rtol=options->tolerance;           /* relative tolerance */
   const dim_t maxit=options->iter_max;            /* max iteration counter */
   const dim_t lmaxit=options->inner_iter_max;     /* max inner iteration counter */
   const bool_t adapt_inner_tolerance=options->adapt_inner_tolerance;
   const double max_inner_tolerance=options->inner_tolerance;  /*inner tolerance counter */
   double inner_tolerance=max_inner_tolerance;
  /*
   * max_inner_tolerance = Maximum error tolerance for residual in inner
   *              iteration. The inner iteration terminates
   *              when the relative linear residual is smaller than inner_tolerance*| F(x_c) |. inner_tolerance is determined
   *              by the modified Eisenstat-Walker formula, if adapt_inner_tolerance is set, otherwise 
   *              inner_tolerance = max_inner_tolerance iteration.
   */
  
  f=TMPMEMALLOC(n,double);
  step=TMPMEMALLOC(n,double);
  if (f==NULL || step ==NULL) {
      Status=SOLVER_MEMORY_ERROR;
  } else {
     /* 
      * initial evaluation of F
      */
     Paso_FunctionCall(F,f,x);
     norm_f=Paso_l2(n,f,F->mpi_info);
     /*
      * stoping criterium:
      */
     stop_tol=atol + rtol*norm_f;
     iteration_count=0;
     if (debug) printf("Start Jacobi-free Newton scheme:\n\ttolerance = %e\n\tstoping tolerance = %e\n",rtol,stop_tol);
     /* 
      *  main iteration loop
      */
     while (! (convergeFlag || maxIterFlag || breakFlag)) {
         /* 
          * keep track of the ratio (reduction_f = norm_f/frnmo) of successive residual norms and 
          * the iteration counter (iteration_count)
          */
         iteration_count++;
         if (debug) printf("iteration step %d: norm of F =%lg\n",iteration_count,norm_f);
         /*
          * call GMRES to get increment
          */
         gmres_iter=lmaxit;
         gmres_tol=inner_tolerance;
         if (debug) printf("GMRES called with tolerance = %e (max iter=%d)\n",inner_tolerance,gmres_iter);
         Status=Paso_Solver_GMRES2(F,f,x,step,&gmres_iter,&gmres_tol,pp);
         if (debug) printf("GMRES finalized after %d steps (residual = %e)\n",gmres_iter,gmres_tol);
         iteration_count+=gmres_iter;
         if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) {
            Status=SOLVER_NO_ERROR;
            /* 
             * update x:
             */
            norm_fo=norm_f; 
            Paso_Update(n,1.,x,1.,step);
            Paso_FunctionCall(F,f,x);
            iteration_count++;
            norm_f=Paso_l2(n,f,F->mpi_info);
            reduction_f=norm_f/norm_fo;
            /*
             *   adjust inner_tolerance 
             */
             if (adapt_inner_tolerance) {
                 old_inner_tolerance=inner_tolerance;
                 inner_tolerance=inner_tolerance_safety * reduction_f * reduction_f;
                 rtmp=inner_tolerance_safety * old_inner_tolerance * old_inner_tolerance;
                 if (rtmp>.1) inner_tolerance=MAX(inner_tolerance,rtmp);
                 inner_tolerance=MAX(MIN(inner_tolerance,max_inner_tolerance), .5*stop_tol/norm_f);
              }
              convergeFlag = (norm_f <= stop_tol);
         } else {
           breakFlag=TRUE;
         }
         maxIterFlag = (iteration_count > maxit);
      }
  }
  TMPMEMFREE(f);
  TMPMEMFREE(step);
  return Status;
}
1	gross	1476	/* $Id:$ */
2
3			/*******************************************************
4			*
5			* Copyright 2008 by University of Queensland
6			*
7			* http://esscc.uq.edu.au
8			* Primary Business: Queensland, Australia
9			* Licensed under the Open Software License version 3.0
10			* http://www.opensource.org/licenses/osl-3.0.php
11			*
12			*******************************************************/
13
14			/*
15			* Purpose
16			* =======
17			*
18			* NewtonGMRES solves the non-linear system F(x)=0 using the
19			* restarted GMRES method.
20			*
21			* Convergence test: norm(F(x))< rtol*norm(F(x0))+atol
22			*
23			* x input initial guess
24			* output solution approximation
25			*
26			* based on code by C. T. Kelley, July 1, 1994.
27			*
28			*/
29
30			#include "Common.h"
31			#include "Solver.h"
32	gross	1639	#include "PasoUtil.h"
33	gross	1476	#ifdef _OPENMP
34			#include <omp.h>
35			#endif
36			err_t Paso_Solver_NewtonGMRES(
37	gross	1639	Paso_Function F, / function evaluation */
38			double x, / in: initial guess, out: new approximation */
39	gross	1476	Paso_Options* options,
40			Paso_Performance* pp)
41
42			{
43	gross	1639	const double inner_tolerance_safety=.9;
44	gross	1476	dim_t gmres_iter;
45	gross	1639	double stop_tol, norm_f,norm_fo, reduction_f,old_inner_tolerance, gmres_tol, rtmp;
46	gross	1476	bool_t convergeFlag=FALSE, maxIterFlag=FALSE, breakFlag=FALSE;
47			double f=NULL, step=NULL;
48			err_t Status=SOLVER_NO_ERROR;
49	gross	1639	const bool_t debug=options->verbose;
50			const dim_t n = F->n;
51			dim_t iteration_count=0;
52			const double atol=options->absolute_tolerance; /* absolute tolerance */
53			const double rtol=options->tolerance; /* relative tolerance */
54			const dim_t maxit=options->iter_max; /* max iteration counter */
55			const dim_t lmaxit=options->inner_iter_max; /* max inner iteration counter */
56			const bool_t adapt_inner_tolerance=options->adapt_inner_tolerance;
57			const double max_inner_tolerance=options->inner_tolerance; /inner tolerance counter /
58			double inner_tolerance=max_inner_tolerance;
59	gross	1476	/*
60	gross	1639	* max_inner_tolerance = Maximum error tolerance for residual in inner
61	gross	1476	* iteration. The inner iteration terminates
62	gross	1639	* when the relative linear residual is smaller than inner_tolerance*\| F(x_c) \|. inner_tolerance is determined
63			* by the modified Eisenstat-Walker formula, if adapt_inner_tolerance is set, otherwise
64			* inner_tolerance = max_inner_tolerance iteration.
65	gross	1476	*/
66
67			f=TMPMEMALLOC(n,double);
68			step=TMPMEMALLOC(n,double);
69			if (f==NULL \|\| step ==NULL) {
70			Status=SOLVER_MEMORY_ERROR;
71			} else {
72			/*
73			* initial evaluation of F
74			*/
75			Paso_FunctionCall(F,f,x);
76	gross	1639	norm_f=Paso_l2(n,f,F->mpi_info);
77	gross	1476	/*
78			* stoping criterium:
79			*/
80	gross	1639	stop_tol=atol + rtol*norm_f;
81			iteration_count=0;
82	gross	1476	if (debug) printf("Start Jacobi-free Newton scheme:\n\ttolerance = %e\n\tstoping tolerance = %e\n",rtol,stop_tol);
83			/*
84			* main iteration loop
85			*/
86			while (! (convergeFlag \|\| maxIterFlag \|\| breakFlag)) {
87			/*
88	gross	1639	* keep track of the ratio (reduction_f = norm_f/frnmo) of successive residual norms and
89			* the iteration counter (iteration_count)
90	gross	1476	*/
91	gross	1639	iteration_count++;
92	phornby	1651	if (debug) printf("iteration step %d: norm of F =%lg\n",iteration_count,norm_f);
93	gross	1476	/*
94			* call GMRES to get increment
95			*/
96			gmres_iter=lmaxit;
97	gross	1639	gmres_tol=inner_tolerance;
98			if (debug) printf("GMRES called with tolerance = %e (max iter=%d)\n",inner_tolerance,gmres_iter);
99			Status=Paso_Solver_GMRES2(F,f,x,step,&gmres_iter,&gmres_tol,pp);
100			if (debug) printf("GMRES finalized after %d steps (residual = %e)\n",gmres_iter,gmres_tol);
101			iteration_count+=gmres_iter;
102	gross	1476	if ((Status==SOLVER_NO_ERROR) \|\| (Status==SOLVER_MAXITER_REACHED)) {
103	gross	1511	Status=SOLVER_NO_ERROR;
104	gross	1476	/*
105			* update x:
106			*/
107	gross	1639	norm_fo=norm_f;
108	gross	1476	Paso_Update(n,1.,x,1.,step);
109			Paso_FunctionCall(F,f,x);
110	gross	1639	iteration_count++;
111			norm_f=Paso_l2(n,f,F->mpi_info);
112			reduction_f=norm_f/norm_fo;
113	gross	1476	/*
114	gross	1639	* adjust inner_tolerance
115	gross	1476	*/
116	gross	1639	if (adapt_inner_tolerance) {
117			old_inner_tolerance=inner_tolerance;
118			inner_tolerance=inner_tolerance_safety * reduction_f * reduction_f;
119			rtmp=inner_tolerance_safety * old_inner_tolerance * old_inner_tolerance;
120			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	gross	1476	} else {
125			breakFlag=TRUE;
126			}
127	gross	1639	maxIterFlag = (iteration_count > maxit);
128	gross	1476	}
129			}
130			TMPMEMFREE(f);
131			TMPMEMFREE(step);
132			return Status;
133			}