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	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	artak	1480	#include "Util.h"
33	gross	1476	#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			}