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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1639 - (hide annotations)
Mon Jul 14 08:55:25 2008 UTC (12 years, 3 months ago) by gross
File MIME type: text/plain
File size: 4815 byte(s)


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

  ViewVC Help
Powered by ViewVC 1.1.26