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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4873 - (show annotations)
Wed Apr 16 06:38:51 2014 UTC (5 years, 5 months ago) by caltinay
File size: 6409 byte(s)
whitespace only changes.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /*
19 * Purpose
20 * =======
21 *
22 * NewtonGMRES solves the nonlinear system F(x)=0 using the
23 * restarted GMRES method.
24 *
25 * Convergence test: norm(F(x))< rtol*norm(F(x0))+atol
26 *
27 * x input initial guess
28 * output solution approximation
29 *
30 * Based on code by C. T. Kelley, July 1, 1994.
31 *
32 */
33
34 #include "Solver.h"
35 #include "PasoUtil.h"
36
37 namespace paso {
38
39 err_t Solver_NewtonGMRES(Function* F, double* x, Options* options,
40 Performance* pp)
41 {
42 const double inner_tolerance_safety=.9;
43 dim_t gmres_iter;
44 double stop_tol, norm2_f,norm2_fo, normsup_f,reduction_f, gmres_tol, rtmp, quad_tolerance;
45 bool convergeFlag=false, maxIterFlag=false, breakFlag=false;
46 double *f=NULL, *step=NULL;
47 err_t Status=SOLVER_NO_ERROR;
48 const bool debug=options->verbose;
49 const dim_t n = F->getLen();
50 dim_t iteration_count=0;
51 const double atol=options->absolute_tolerance; /* absolute tolerance */
52 const double rtol=options->tolerance; /* relative tolerance */
53 const dim_t maxit=options->iter_max; /* max iteration counter */
54 const dim_t lmaxit=options->inner_iter_max*10; /* max inner iteration counter */
55 const bool adapt_inner_tolerance=options->adapt_inner_tolerance;
56 const double max_inner_tolerance=options->inner_tolerance *1.e-11;
57 double inner_tolerance=max_inner_tolerance;
58 /*
59 * max_inner_tolerance = Maximum error tolerance for residual in inner
60 * iteration. The inner iteration terminates
61 * when the relative linear residual is smaller than
62 * inner_tolerance*| F(x_c) |. Inner_tolerance is determined
63 * by the modified Eisenstat-Walker formula, if
64 * adapt_inner_tolerance is set, otherwise
65 * inner_tolerance = max_inner_tolerance iteration.
66 */
67
68 f=new double[n];
69 step=new double[n];
70 /*
71 * initial evaluation of F
72 */
73 F->call(f, x, pp);
74 iteration_count++;
75 norm2_f=util::l2(n,f,F->mpi_info);
76 normsup_f=util::lsup(n,f,F->mpi_info);
77 /*
78 * stopping criterion:
79 */
80 stop_tol=atol + rtol*normsup_f;
81 if (stop_tol<=0) {
82 Status=SOLVER_INPUT_ERROR;
83 if (debug) printf("NewtonGMRES: zero tolerance given.\n");
84 } else {
85 iteration_count=1;
86 if (debug) {
87 printf("NewtonGMRES: Start Jacobi-free Newton scheme\n");
88 printf("NewtonGMRES: lsup tolerance rel/abs= %e/%e\n",rtol,atol);
89 printf("NewtonGMRES: lsup stopping tolerance = %e\n",stop_tol);
90 printf("NewtonGMRES: max. inner iterations (GMRES) = %d\n",lmaxit);
91 if (adapt_inner_tolerance) {
92 printf("NewtonGMRES: inner tolerance is adapted.\n");
93 printf("NewtonGMRES: max. inner l2 tolerance (GMRES) = %e\n",max_inner_tolerance);
94 } else {
95 printf("NewtonGMRES: inner l2 tolerance (GMRES) = %e\n",inner_tolerance);
96 }
97 }
98 /*
99 * main iteration loop
100 */
101 while (! (convergeFlag || maxIterFlag || breakFlag)) {
102 // keep track of the ratio (reduction_f = norm2_f/norm2_fo) of
103 // successive residual norms and the iteration counter
104 // (iteration_count)
105 if (debug)
106 printf("NewtonGMRES: iteration step %d: lsup-norm of F =%e\n",iteration_count,normsup_f);
107
108 // call GMRES to get increment
109 gmres_iter=lmaxit;
110 gmres_tol=inner_tolerance;
111 Status = Solver_GMRES2(F, f, x, step, &gmres_iter, &gmres_tol, pp);
112 inner_tolerance=std::max(inner_tolerance, gmres_tol/norm2_f);
113 printf("NewtonGMRES: actual rel. inner tolerance = %e\n",inner_tolerance);
114 iteration_count+=gmres_iter;
115 if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) {
116 Status=SOLVER_NO_ERROR;
117 // update x
118 norm2_fo=norm2_f;
119 util::update(n,1.,x,1.,step);
120 F->call(f, x, pp);
121 iteration_count++;
122 norm2_f=util::l2(n,f,F->mpi_info);
123 normsup_f=util::lsup(n,f,F->mpi_info);
124 reduction_f=norm2_f/norm2_fo;
125 // adjust inner_tolerance
126 if (adapt_inner_tolerance) {
127 quad_tolerance = inner_tolerance_safety * reduction_f * reduction_f;
128 rtmp=inner_tolerance_safety * inner_tolerance * inner_tolerance;
129 if (rtmp>.1) {
130 inner_tolerance=std::min(max_inner_tolerance, std::max(quad_tolerance,rtmp));
131 } else {
132 inner_tolerance=std::min(max_inner_tolerance, quad_tolerance);
133 }
134 inner_tolerance=std::min(max_inner_tolerance, std::max(inner_tolerance, .5*stop_tol/normsup_f));
135 }
136 convergeFlag = (normsup_f <= stop_tol);
137 } else {
138 breakFlag = true;
139 }
140 maxIterFlag = (iteration_count > maxit);
141 }
142 if (debug) {
143 printf("NewtonGMRES: iteration step %d: lsup-norm of F =%e\n",iteration_count,normsup_f);
144
145 if (convergeFlag) printf("NewtonGMRES: convergence reached after %d steps.\n",iteration_count);
146 if (breakFlag) printf("NewtonGMRES: iteration break down after %d steps.\n",iteration_count);
147 if (maxIterFlag) printf("NewtonGMRES: maximum number of iteration steps %d reached.\n",maxit);
148 }
149 if (breakFlag) Status = SOLVER_BREAKDOWN;
150 if (maxIterFlag) Status = SOLVER_MAXITER_REACHED;
151 }
152 delete[] f;
153 delete[] step;
154 if (debug) printf("NewtonGMRES: STATUS return = %d\n",Status);
155 return Status;
156 }
157
158 } // namespace paso
159

Properties

Name Value
svn:mergeinfo /branches/amg_from_3530/paso/src/NewtonGMRES.cpp:3531-3826 /branches/lapack2681/paso/src/NewtonGMRES.cpp:2682-2741 /branches/pasowrap/paso/src/NewtonGMRES.cpp:3661-3674 /branches/py3_attempt2/paso/src/NewtonGMRES.cpp:3871-3891 /branches/restext/paso/src/NewtonGMRES.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/NewtonGMRES.cpp:3669-3791 /branches/stage3.0/paso/src/NewtonGMRES.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/NewtonGMRES.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/NewtonGMRES.cpp:3517-3974 /release/3.0/paso/src/NewtonGMRES.cpp:2591-2601 /trunk/paso/src/NewtonGMRES.cpp:4257-4344 /trunk/ripley/test/python/paso/src/NewtonGMRES.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26