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

Revision 1798 - (show annotations)
Wed Sep 17 06:21:12 2008 UTC (12 years, 1 month ago) by gross
File MIME type: text/plain
File size: 5258 byte(s)
```Fixes for the JacobeanFreeNewton scheme. Still needs to be tested under OPENMP but runs under 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 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 %s 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 }