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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1637 - (hide annotations)
Mon Jul 14 05:34:59 2008 UTC (12 years, 3 months ago) by ksteube
File MIME type: text/plain
File size: 3895 byte(s)
Resolved some compiler warnings
Changed blocktimer to not use strdup, intead malloc and strcpy

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 artak 1480 #include "Util.h"
33 gross 1476 #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 ksteube 1637 if (debug) printf("iteration step %d: norm of F =%e\n",itc,fnrm);
95 gross 1476 /*
96     * call GMRES to get increment
97     */
98     gmres_iter=lmaxit;
99     gmres_tol=eta;
100 ksteube 1637 if (debug) printf("GMRES called with tolerance = %e\n",eta);
101 gross 1476 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 gross 1511 Status=SOLVER_NO_ERROR;
105 gross 1476 /*
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