/[escript]/trunk/paso/src/NewtonGMRES.c
ViewVC logotype

Contents of /trunk/paso/src/NewtonGMRES.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3642 - (show annotations)
Thu Oct 27 03:41:51 2011 UTC (7 years, 5 months ago) by caltinay
File MIME type: text/plain
File size: 6205 byte(s)
Assorted spelling/comment fixes in paso.

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 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 nonlinear 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, norm2_f,norm2_fo, normsup_f,reduction_f, gmres_tol, rtmp, quad_tolerance;
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*10; /* max inner iteration counter */
57 const bool_t adapt_inner_tolerance=options->adapt_inner_tolerance;
58 const double max_inner_tolerance=options->inner_tolerance *1.e-11;
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
64 * inner_tolerance*| F(x_c) |. Inner_tolerance is determined
65 * by the modified Eisenstat-Walker formula, if
66 * adapt_inner_tolerance is set, otherwise
67 * inner_tolerance = max_inner_tolerance iteration.
68 */
69
70 f=TMPMEMALLOC(n,double);
71 step=TMPMEMALLOC(n,double);
72 if (f==NULL || step ==NULL) {
73 Status=SOLVER_MEMORY_ERROR;
74 } else {
75 /*
76 * initial evaluation of F
77 */
78 Paso_FunctionCall(F,f,x,pp);
79 iteration_count++;
80 norm2_f=Paso_l2(n,f,F->mpi_info);
81 normsup_f=Paso_lsup(n,f,F->mpi_info);
82 /*
83 * stopping criterion:
84 */
85 stop_tol=atol + rtol*normsup_f;
86 if (stop_tol<=0) {
87 Status=SOLVER_INPUT_ERROR;
88 if (debug) printf("NewtonGMRES: zero tolerance given.\n");
89 } else {
90 iteration_count=1;
91 if (debug) {
92 printf("NewtonGMRES: Start Jacobi-free Newton scheme\n");
93 printf("NewtonGMRES: lsup tolerance rel/abs= %e/%e\n",rtol,atol);
94 printf("NewtonGMRES: lsup stopping tolerance = %e\n",stop_tol);
95
96 printf("NewtonGMRES: max. inner iterations (GMRES) = %d\n",lmaxit);
97 if (adapt_inner_tolerance) {
98 printf("NewtonGMRES: inner tolerance is adapted.\n");
99 printf("NewtonGMRES: max. inner l2 tolerance (GMRES) = %e\n",max_inner_tolerance);
100 } else {
101 printf("NewtonGMRES: inner l2 tolerance (GMRES) = %e\n",inner_tolerance);
102 }
103 }
104 /*
105 * main iteration loop
106 */
107 while (! (convergeFlag || maxIterFlag || breakFlag)) {
108 /*
109 * keep track of the ratio (reduction_f = norm2_f/norm2_fo) of
110 * successive residual norms and the iteration counter (iteration_count)
111 */
112 if (debug) printf("NewtonGMRES: iteration step %d: lsup-norm of F =%e\n",iteration_count,normsup_f);
113 /*
114 * call GMRES to get increment
115 */
116 gmres_iter=lmaxit;
117 gmres_tol=inner_tolerance;
118 Status=Paso_Solver_GMRES2(F,f,x,step,&gmres_iter,&gmres_tol,pp);
119 inner_tolerance=MAX(inner_tolerance, gmres_tol/norm2_f);
120 printf("NewtonGMRES: actual rel. inner tolerance = %e\n",inner_tolerance);
121 iteration_count+=gmres_iter;
122 if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) {
123 Status=SOLVER_NO_ERROR;
124 /*
125 * update x:
126 */
127 norm2_fo=norm2_f;
128 Paso_Update(n,1.,x,1.,step);
129 Paso_FunctionCall(F,f,x,pp);
130 iteration_count++;
131 norm2_f=Paso_l2(n,f,F->mpi_info);
132 normsup_f=Paso_lsup(n,f,F->mpi_info);
133 reduction_f=norm2_f/norm2_fo;
134 /*
135 * adjust inner_tolerance
136 */
137 if (adapt_inner_tolerance) {
138 quad_tolerance = inner_tolerance_safety * reduction_f * reduction_f;
139 rtmp=inner_tolerance_safety * inner_tolerance * inner_tolerance;
140 if (rtmp>.1) {
141 inner_tolerance=MIN(max_inner_tolerance, MAX(quad_tolerance,rtmp));
142 } else {
143 inner_tolerance=MIN(max_inner_tolerance, quad_tolerance);
144 }
145 inner_tolerance=MIN(max_inner_tolerance, MAX(inner_tolerance, .5*stop_tol/normsup_f));
146 }
147 convergeFlag = (normsup_f <= stop_tol);
148 } else {
149 breakFlag=TRUE;
150 }
151 maxIterFlag = (iteration_count > maxit);
152 }
153 if (debug) {
154 printf("NewtonGMRES: iteration step %d: lsup-norm of F =%e\n",iteration_count,normsup_f);
155
156 if (convergeFlag) printf("NewtonGMRES: convergence reached after %d steps.\n",iteration_count);
157 if (breakFlag) printf("NewtonGMRES: iteration break down after %d steps.\n",iteration_count);
158 if (maxIterFlag) printf("NewtonGMRES: maximum number of iteration steps %d reached.\n",maxit);
159 }
160 if (breakFlag) Status=SOLVER_BREAKDOWN;
161 if (maxIterFlag) Status=SOLVER_MAXITER_REACHED;
162 }
163 }
164 TMPMEMFREE(f);
165 TMPMEMFREE(step);
166 if (debug) printf("NewtonGMRES: STATUS return = %d\n",Status);
167 return Status;
168 }
169

  ViewVC Help
Powered by ViewVC 1.1.26