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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1801 - (show annotations)
Fri Sep 19 01:37:09 2008 UTC (12 years, 1 month ago) by ksteube
File MIME type: text/plain
File size: 5258 byte(s)
Fixed serialization of I/O for MPI...code didn't compile without MPI

1 /* $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 #include "PasoUtil.h"
33 #ifdef _OPENMP
34 #include <omp.h>
35 #endif
36 err_t Paso_Solver_NewtonGMRES(
37 Paso_Function *F, /* function evaluation */
38 double *x, /* in: initial guess, out: new approximation */
39 Paso_Options* options,
40 Paso_Performance* pp)
41
42 {
43 const double inner_tolerance_safety=.9;
44 dim_t gmres_iter;
45 double stop_tol, norm_f,norm_fo, reduction_f,old_inner_tolerance, gmres_tol, rtmp;
46 bool_t convergeFlag=FALSE, maxIterFlag=FALSE, breakFlag=FALSE;
47 double *f=NULL, *step=NULL;
48 err_t Status=SOLVER_NO_ERROR;
49 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 /*
60 * max_inner_tolerance = Maximum error tolerance for residual in inner
61 * iteration. The inner iteration terminates
62 * 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 */
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 norm_f=Paso_l2(n,f,F->mpi_info);
77 /*
78 * stoping criterium:
79 */
80 stop_tol=atol + rtol*norm_f;
81 iteration_count=1;
82 if (debug) printf("Start Jacobi-free Newton scheme:\n\ttolerance = %e\n\tstopping tolerance = %e\n",rtol,stop_tol);
83 /*
84 * main iteration loop
85 */
86 while (! (convergeFlag || maxIterFlag || breakFlag)) {
87 /*
88 * keep track of the ratio (reduction_f = norm_f/frnmo) of successive residual norms and
89 * the iteration counter (iteration_count)
90 */
91 if (debug) printf("iteration step %d: norm of F =%lg\n",iteration_count,norm_f);
92 /*
93 * call GMRES to get increment
94 */
95 gmres_iter=lmaxit;
96 gmres_tol=inner_tolerance;
97 if (debug) printf("GMRES called with tolerance = %e (max iter=%d)\n",inner_tolerance,gmres_iter);
98 Status=Paso_Solver_GMRES2(F,f,x,step,&gmres_iter,&gmres_tol,pp);
99 if (debug) printf("GMRES finalized after %d steps (residual = %e)\n",gmres_iter,gmres_tol);
100 iteration_count+=gmres_iter;
101 if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) {
102 Status=SOLVER_NO_ERROR;
103 /*
104 * update x:
105 */
106 norm_fo=norm_f;
107 Paso_Update(n,1.,x,1.,step);
108 Paso_FunctionCall(F,f,x);
109 iteration_count++;
110 norm_f=Paso_l2(n,f,F->mpi_info);
111 reduction_f=norm_f/norm_fo;
112 /*
113 * adjust inner_tolerance
114 */
115 if (adapt_inner_tolerance) {
116 old_inner_tolerance=inner_tolerance;
117 inner_tolerance=inner_tolerance_safety * reduction_f * reduction_f;
118 rtmp=inner_tolerance_safety * old_inner_tolerance * old_inner_tolerance;
119 if (rtmp>.1) inner_tolerance=MAX(inner_tolerance,rtmp);
120 inner_tolerance=MAX(MIN(inner_tolerance,max_inner_tolerance), .5*stop_tol/norm_f);
121 }
122 convergeFlag = (norm_f <= stop_tol);
123 } else {
124 breakFlag=TRUE;
125 }
126 maxIterFlag = (iteration_count > maxit);
127 }
128 if (debug) {
129 if (convergeFlag) printf("convergence reached after %d steps with residual %e.\n",iteration_count,norm_f);
130 if (breakFlag) printf("iteration break down after %d steps.\n",iteration_count);
131 if (maxIterFlag) printf("maximum number of iteration step %d is reached.\n",maxit);
132 }
133 if (breakFlag) Status=SOLVER_BREAKDOWN;
134 if (maxIterFlag) Status=SOLVER_MAXITER_REACHED;
135 }
136 TMPMEMFREE(f);
137 TMPMEMFREE(step);
138 printf("STATUS return = %d\n",Status);
139 return Status;
140 }

  ViewVC Help
Powered by ViewVC 1.1.26