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	/* $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	#include "PasoUtil.h"
33	#ifdef _OPENMP
34	#include <omp.h>
35	#endif
36	err_t Paso_Solver_NewtonGMRES(
37	Paso_Function F, / function evaluation */
38	double x, / in: initial guess, out: new approximation */
39	Paso_Options* options,
40	Paso_Performance* pp)
41
42	{
43	const double inner_tolerance_safety=.9;
44	dim_t gmres_iter;
45	double stop_tol, norm_f,norm_fo, reduction_f,old_inner_tolerance, gmres_tol, rtmp;
46	bool_t convergeFlag=FALSE, maxIterFlag=FALSE, breakFlag=FALSE;
47	double f=NULL, step=NULL;
48	err_t Status=SOLVER_NO_ERROR;
49	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	/*
60	* max_inner_tolerance = Maximum error tolerance for residual in inner
61	* iteration. The inner iteration terminates
62	* 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	*/
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	norm_f=Paso_l2(n,f,F->mpi_info);
77	/*
78	* stoping criterium:
79	*/
80	stop_tol=atol + rtol*norm_f;
81	iteration_count=0;
82	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	* keep track of the ratio (reduction_f = norm_f/frnmo) of successive residual norms and
89	* the iteration counter (iteration_count)
90	*/
91	iteration_count++;
92	if (debug) printf("iteration step %d: norm of F =%lg\n",iteration_count,norm_f);
93	/*
94	* call GMRES to get increment
95	*/
96	gmres_iter=lmaxit;
97	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	if ((Status==SOLVER_NO_ERROR) \|\| (Status==SOLVER_MAXITER_REACHED)) {
103	Status=SOLVER_NO_ERROR;
104	/*
105	* update x:
106	*/
107	norm_fo=norm_f;
108	Paso_Update(n,1.,x,1.,step);
109	Paso_FunctionCall(F,f,x);
110	iteration_count++;
111	norm_f=Paso_l2(n,f,F->mpi_info);
112	reduction_f=norm_f/norm_fo;
113	/*
114	* adjust inner_tolerance
115	*/
116	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	} else {
125	breakFlag=TRUE;
126	}
127	maxIterFlag = (iteration_count > maxit);
128	}
129	}
130	TMPMEMFREE(f);
131	TMPMEMFREE(step);
132	return Status;
133	}