paso/src/NewtonGMRES.c


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