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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1480 - (show annotations)
Tue Apr 8 02:02:48 2008 UTC (11 years, 3 months ago) by artak
File MIME type: text/plain
File size: 3896 byte(s)
Unil.h header is added
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 "Util.h"
33 #ifdef _OPENMP
34 #include <omp.h>
35 #endif
36 err_t Paso_Solver_NewtonGMRES(
37 Paso_Function *F,
38 double *x,
39 Paso_Options* options,
40 Paso_Performance* pp)
41
42 {
43 const double gamma=.9;
44 dim_t gmres_iter;
45 double stop_tol, fnrm,fnrmo, rat,etaold,etanew, gmres_tol;
46 bool_t convergeFlag=FALSE, maxIterFlag=FALSE, breakFlag=FALSE;
47 double *f=NULL, *step=NULL;
48 err_t Status=SOLVER_NO_ERROR;
49 bool_t debug=options->verbose;
50 dim_t n = n=F->local_n;
51 dim_t itc=0;
52 double atol=options->absolute_tolerance;
53 double rtol=options->tolerance;
54 dim_t maxit=options->iter_max;
55 dim_t lmaxit=options->inner_iter_max;
56 bool_t adapt_eta=options->adapt_inner_tolerance;
57 double etamax=options->inner_tolerance;
58 double eta=etamax;
59 /*
60 * etamax = Maximum error tolerance for residual in inner
61 * iteration. The inner iteration terminates
62 * when the relative linear residual is
63 * smaller than eta*| F(x_c) |. eta is determined
64 * by the modified Eisenstat-Walker formula if etamax > 0.
65 * If etamax < 0, then eta = |etamax| for the entire
66 * iteration.
67 */
68
69 f=TMPMEMALLOC(n,double);
70 step=TMPMEMALLOC(n,double);
71 if (f==NULL || step ==NULL) {
72 Status=SOLVER_MEMORY_ERROR;
73 } else {
74 /*
75 * initial evaluation of F
76 */
77 Paso_FunctionCall(F,f,x);
78 fnrm=Paso_l2(n,f,F->mpi_info);
79 /*
80 * stoping criterium:
81 */
82 stop_tol=atol + rtol*fnrm;
83 itc=0;
84 if (debug) printf("Start Jacobi-free Newton scheme:\n\ttolerance = %e\n\tstoping tolerance = %e\n",rtol,stop_tol);
85 /*
86 * main iteration loop
87 */
88 while (! (convergeFlag || maxIterFlag || breakFlag)) {
89 /*
90 * keep track of the ratio (rat = fnrm/frnmo) of successive residual norms and
91 * the iteration counter (itc)
92 */
93 itc++;
94 if (debug) printf("iteration step %d: norm of F =%d\n",itc,fnrm);
95 /*
96 * call GMRES to get increment
97 */
98 gmres_iter=lmaxit;
99 gmres_tol=eta;
100 if (debug) printf("GMRES called with tolerance = %d\n",eta);
101 Status=Paso_Solver_NLGMRES(F,f,x,step,&gmres_iter,&gmres_tol,pp);
102 itc+=gmres_iter;
103 if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) {
104 Status==SOLVER_NO_ERROR;
105 /*
106 * update x:
107 */
108 Paso_Update(n,1.,x,1.,step);
109 Paso_FunctionCall(F,f,x);
110 fnrm=Paso_l2(n,f,F->mpi_info);
111 fnrmo=fnrm;
112 rat=fnrm/fnrmo;
113 /*
114 * adjust eta
115 */
116 if (adapt_eta) {
117 etaold=eta;
118 etanew=gamma*rat*rat;
119 if (gamma*etaold*etaold >1.) etanew=MAX(etanew,gamma*etaold*etaold);
120 eta=MAX(MIN(etanew,etamax), .5*stop_tol/fnrm);
121 }
122 } else {
123 breakFlag=TRUE;
124 }
125 maxIterFlag = (itc > maxit);
126 convergeFlag = (fnrm <= stop_tol);
127 }
128 }
129 TMPMEMFREE(f);
130 TMPMEMFREE(step);
131 return Status;
132 }

  ViewVC Help
Powered by ViewVC 1.1.26