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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1887 - (hide annotations)
Wed Oct 15 03:26:25 2008 UTC (10 years, 9 months ago) by ksteube
File MIME type: text/plain
File size: 5443 byte(s)
Fixed two typos that stopped the test suite from running.

Also, gcc 4.3.2 issued several warnings not seen before:
ignoring the return value of fscanf and using the wrong format
with printf.

1 gross 1476
2     /*******************************************************
3 ksteube 1811 *
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 gross 1476
14 ksteube 1811
15 gross 1476 /*
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 gross 1639 #include "PasoUtil.h"
34 gross 1476 #ifdef _OPENMP
35     #include <omp.h>
36     #endif
37     err_t Paso_Solver_NewtonGMRES(
38 gross 1639 Paso_Function *F, /* function evaluation */
39     double *x, /* in: initial guess, out: new approximation */
40 gross 1476 Paso_Options* options,
41     Paso_Performance* pp)
42    
43     {
44 gross 1639 const double inner_tolerance_safety=.9;
45 gross 1476 dim_t gmres_iter;
46 gross 1639 double stop_tol, norm_f,norm_fo, reduction_f,old_inner_tolerance, gmres_tol, rtmp;
47 gross 1476 bool_t convergeFlag=FALSE, maxIterFlag=FALSE, breakFlag=FALSE;
48     double *f=NULL, *step=NULL;
49     err_t Status=SOLVER_NO_ERROR;
50 gross 1639 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 gross 1476 /*
61 gross 1639 * max_inner_tolerance = Maximum error tolerance for residual in inner
62 gross 1476 * iteration. The inner iteration terminates
63 gross 1639 * 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 gross 1476 */
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 gross 1804 Paso_FunctionCall(F,f,x,pp);
77 gross 1639 norm_f=Paso_l2(n,f,F->mpi_info);
78 gross 1476 /*
79     * stoping criterium:
80     */
81 gross 1639 stop_tol=atol + rtol*norm_f;
82 gross 1804 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 gross 1476 /*
93 gross 1639 * keep track of the ratio (reduction_f = norm_f/frnmo) of successive residual norms and
94     * the iteration counter (iteration_count)
95 gross 1476 */
96 ksteube 1887 if (debug) printf("iteration step %d: norm of F =%g\n",iteration_count,norm_f);
97 gross 1476 /*
98     * call GMRES to get increment
99     */
100     gmres_iter=lmaxit;
101 gross 1639 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 gross 1476 if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) {
107 gross 1511 Status=SOLVER_NO_ERROR;
108 gross 1476 /*
109     * update x:
110     */
111 gross 1639 norm_fo=norm_f;
112 gross 1476 Paso_Update(n,1.,x,1.,step);
113 gross 1804 Paso_FunctionCall(F,f,x,pp);
114 gross 1639 iteration_count++;
115     norm_f=Paso_l2(n,f,F->mpi_info);
116     reduction_f=norm_f/norm_fo;
117 gross 1476 /*
118 gross 1639 * adjust inner_tolerance
119 gross 1476 */
120 gross 1798 if (adapt_inner_tolerance) {
121 gross 1639 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 gross 1798 }
127     convergeFlag = (norm_f <= stop_tol);
128 gross 1476 } else {
129 gross 1798 breakFlag=TRUE;
130 gross 1476 }
131 gross 1639 maxIterFlag = (iteration_count > maxit);
132 gross 1804 }
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 gross 1476 }
141     }
142     TMPMEMFREE(f);
143     TMPMEMFREE(step);
144 gross 1798 printf("STATUS return = %d\n",Status);
145 gross 1476 return Status;
146     }

  ViewVC Help
Powered by ViewVC 1.1.26