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 |
} |