# Contents of /branches/arrexp_2137_win/paso/src/NewtonGMRES.c

Revision 2202 - (show annotations)
Fri Jan 9 01:28:32 2009 UTC (10 years, 4 months ago) by jfenwick
File MIME type: text/plain
File size: 5443 byte(s)
```Branching the array experiments from version 2137.
The idea is to make the changes required for the c++ changes to compile
on windows without bringing in the later python changes.

```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2008 by University of Queensland 5 * Earth Systems Science Computational Center (ESSCC) 6 * http://www.uq.edu.au/esscc 7 * 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 /* 16 * Purpose 17 * ======= 18 * 19 * NewtonGMRES solves the non-linear system F(x)=0 using the 20 * restarted GMRES method. 21 * 22 * Convergence test: norm(F(x))< rtol*norm(F(x0))+atol 23 * 24 * x input initial guess 25 * output solution approximation 26 * 27 * based on code by C. T. Kelley, July 1, 1994. 28 * 29 */ 30 31 #include "Common.h" 32 #include "Solver.h" 33 #include "PasoUtil.h" 34 #ifdef _OPENMP 35 #include 36 #endif 37 err_t Paso_Solver_NewtonGMRES( 38 Paso_Function *F, /* function evaluation */ 39 double *x, /* in: initial guess, out: new approximation */ 40 Paso_Options* options, 41 Paso_Performance* pp) 42 43 { 44 const double inner_tolerance_safety=.9; 45 dim_t gmres_iter; 46 double stop_tol, norm_f,norm_fo, reduction_f,old_inner_tolerance, gmres_tol, rtmp; 47 bool_t convergeFlag=FALSE, maxIterFlag=FALSE, breakFlag=FALSE; 48 double *f=NULL, *step=NULL; 49 err_t Status=SOLVER_NO_ERROR; 50 const bool_t debug=options->verbose; 51 const dim_t n = F->n; 52 dim_t iteration_count=0; 53 const double atol=options->absolute_tolerance; /* absolute tolerance */ 54 const double rtol=options->tolerance; /* relative tolerance */ 55 const dim_t maxit=options->iter_max; /* max iteration counter */ 56 const dim_t lmaxit=options->inner_iter_max; /* max inner iteration counter */ 57 const bool_t adapt_inner_tolerance=options->adapt_inner_tolerance; 58 const double max_inner_tolerance=options->inner_tolerance; /*inner tolerance counter */ 59 double inner_tolerance=max_inner_tolerance; 60 /* 61 * max_inner_tolerance = Maximum error tolerance for residual in inner 62 * iteration. The inner iteration terminates 63 * when the relative linear residual is smaller than inner_tolerance*| F(x_c) |. inner_tolerance is determined 64 * by the modified Eisenstat-Walker formula, if adapt_inner_tolerance is set, otherwise 65 * inner_tolerance = max_inner_tolerance 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,pp); 77 norm_f=Paso_l2(n,f,F->mpi_info); 78 /* 79 * stoping criterium: 80 */ 81 stop_tol=atol + rtol*norm_f; 82 if (stop_tol<=0) { 83 Status=SOLVER_INPUT_ERROR; 84 if (debug) printf("zero tolerance given.\n"); 85 } else { 86 iteration_count=1; 87 if (debug) printf("Start Jacobi-free Newton scheme:\n\ttolerance = %e\n\tstopping tolerance = %e\n",rtol,stop_tol); 88 /* 89 * main iteration loop 90 */ 91 while (! (convergeFlag || maxIterFlag || breakFlag)) { 92 /* 93 * keep track of the ratio (reduction_f = norm_f/frnmo) of successive residual norms and 94 * the iteration counter (iteration_count) 95 */ 96 if (debug) printf("iteration step %d: norm of F =%g\n",iteration_count,norm_f); 97 /* 98 * call GMRES to get increment 99 */ 100 gmres_iter=lmaxit; 101 gmres_tol=inner_tolerance; 102 if (debug) printf("GMRES called with tolerance = %e (max iter=%d)\n",inner_tolerance,gmres_iter); 103 Status=Paso_Solver_GMRES2(F,f,x,step,&gmres_iter,&gmres_tol,pp); 104 if (debug) printf("GMRES finalized after %d steps (residual = %e)\n",gmres_iter,gmres_tol); 105 iteration_count+=gmres_iter; 106 if ((Status==SOLVER_NO_ERROR) || (Status==SOLVER_MAXITER_REACHED)) { 107 Status=SOLVER_NO_ERROR; 108 /* 109 * update x: 110 */ 111 norm_fo=norm_f; 112 Paso_Update(n,1.,x,1.,step); 113 Paso_FunctionCall(F,f,x,pp); 114 iteration_count++; 115 norm_f=Paso_l2(n,f,F->mpi_info); 116 reduction_f=norm_f/norm_fo; 117 /* 118 * adjust inner_tolerance 119 */ 120 if (adapt_inner_tolerance) { 121 old_inner_tolerance=inner_tolerance; 122 inner_tolerance=inner_tolerance_safety * reduction_f * reduction_f; 123 rtmp=inner_tolerance_safety * old_inner_tolerance * old_inner_tolerance; 124 if (rtmp>.1) inner_tolerance=MAX(inner_tolerance,rtmp); 125 inner_tolerance=MAX(MIN(inner_tolerance,max_inner_tolerance), .5*stop_tol/norm_f); 126 } 127 convergeFlag = (norm_f <= stop_tol); 128 } else { 129 breakFlag=TRUE; 130 } 131 maxIterFlag = (iteration_count > maxit); 132 } 133 if (debug) { 134 if (convergeFlag) printf("convergence reached after %d steps with residual %e.\n",iteration_count,norm_f); 135 if (breakFlag) printf("iteration break down after %d steps.\n",iteration_count); 136 if (maxIterFlag) printf("maximum number of iteration step %d is reached.\n",maxit); 137 } 138 if (breakFlag) Status=SOLVER_BREAKDOWN; 139 if (maxIterFlag) Status=SOLVER_MAXITER_REACHED; 140 } 141 } 142 TMPMEMFREE(f); 143 TMPMEMFREE(step); 144 printf("STATUS return = %d\n",Status); 145 return Status; 146 }