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 "Util.h"
#ifdef _OPENMP
#include <omp.h>
#endif
err_t Paso_Solver_NewtonGMRES(
    Paso_Function *F,
    double *x,
    Paso_Options* options,
    Paso_Performance* pp) 

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