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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26