1e7e93795SLois Curfman McInnes #ifndef lint 2*5d2e0e51SBarry Smith static char vcid[] = "$Id: snesut.c,v 1.7 1996/01/23 00:19:51 bsmith Exp bsmith $"; 3e7e93795SLois Curfman McInnes #endif 4e7e93795SLois Curfman McInnes 5e7e93795SLois Curfman McInnes #include <math.h> 6e7e93795SLois Curfman McInnes #include "snesimpl.h" /*I "snes.h" I*/ 7e7e93795SLois Curfman McInnes 84b828684SBarry Smith /*@C 9e7e93795SLois Curfman McInnes SNESDefaultMonitor - Default SNES monitoring routine. 10e7e93795SLois Curfman McInnes 11e7e93795SLois Curfman McInnes Input Parameters: 12e7e93795SLois Curfman McInnes . snes - the SNES context 13e7e93795SLois Curfman McInnes . its - iteration number 14e7e93795SLois Curfman McInnes . fgnorm - 2-norm of residual (or gradient) 15e7e93795SLois Curfman McInnes . dummy - unused context 16e7e93795SLois Curfman McInnes 17e7e93795SLois Curfman McInnes Notes: 18e7e93795SLois Curfman McInnes For SNES_NONLINEAR_EQUATIONS methods the routine prints the 19e7e93795SLois Curfman McInnes residual norm at each iteration. 20e7e93795SLois Curfman McInnes 21e7e93795SLois Curfman McInnes For SNES_UNCONSTRAINED_MINIMIZATION methods the routine prints the 22e7e93795SLois Curfman McInnes function value and gradient norm at each iteration. 23e7e93795SLois Curfman McInnes 24e7e93795SLois Curfman McInnes .keywords: SNES, nonlinear, default, monitor, norm 25e7e93795SLois Curfman McInnes 26e7e93795SLois Curfman McInnes .seealso: SNESSetMonitor() 27e7e93795SLois Curfman McInnes @*/ 28e7e93795SLois Curfman McInnes int SNESDefaultMonitor(SNES snes,int its,double fgnorm,void *dummy) 29e7e93795SLois Curfman McInnes { 30e7e93795SLois Curfman McInnes if (snes->method_class == SNES_NONLINEAR_EQUATIONS) 31052efed2SBarry Smith MPIU_printf(snes->comm, "iter = %d, SNES Function norm %g \n",its,fgnorm); 32e7e93795SLois Curfman McInnes else if (snes->method_class == SNES_UNCONSTRAINED_MINIMIZATION) 33e7e93795SLois Curfman McInnes MPIU_printf(snes->comm, 34e7e93795SLois Curfman McInnes "iter = %d, Function value %g, Gradient norm %g \n",its,snes->fc,fgnorm); 35e7e93795SLois Curfman McInnes else SETERRQ(1,"SNESDefaultMonitor:Unknown method class"); 36e7e93795SLois Curfman McInnes return 0; 37e7e93795SLois Curfman McInnes } 38e7e93795SLois Curfman McInnes /* ---------------------------------------------------------------- */ 39e7e93795SLois Curfman McInnes int SNESDefaultSMonitor(SNES snes,int its, double fgnorm,void *dummy) 40e7e93795SLois Curfman McInnes { 41e7e93795SLois Curfman McInnes if (snes->method_class == SNES_NONLINEAR_EQUATIONS) { 42e7e93795SLois Curfman McInnes if (fgnorm > 1.e-9 || fgnorm == 0.0) { 43e7e93795SLois Curfman McInnes MPIU_printf(snes->comm, "iter = %d, Function norm %g \n",its,fgnorm); 44e7e93795SLois Curfman McInnes } 45e7e93795SLois Curfman McInnes else if (fgnorm > 1.e-11){ 46e7e93795SLois Curfman McInnes MPIU_printf(snes->comm, "iter = %d, Function norm %5.3e \n",its,fgnorm); 47e7e93795SLois Curfman McInnes } 48e7e93795SLois Curfman McInnes else { 49e7e93795SLois Curfman McInnes MPIU_printf(snes->comm, "iter = %d, Function norm < 1.e-11\n",its); 50e7e93795SLois Curfman McInnes } 51e7e93795SLois Curfman McInnes } else if (snes->method_class == SNES_UNCONSTRAINED_MINIMIZATION) { 52e7e93795SLois Curfman McInnes if (fgnorm > 1.e-9 || fgnorm == 0.0) { 53e7e93795SLois Curfman McInnes MPIU_printf(snes->comm, 54e7e93795SLois Curfman McInnes "iter = %d, Function value %g, Gradient norm %g \n", 55e7e93795SLois Curfman McInnes its,snes->fc,fgnorm); 56e7e93795SLois Curfman McInnes } 57e7e93795SLois Curfman McInnes else if (fgnorm > 1.e-11) { 58e7e93795SLois Curfman McInnes MPIU_printf(snes->comm, 59e7e93795SLois Curfman McInnes "iter = %d, Function value %g, Gradient norm %5.3e \n", 60e7e93795SLois Curfman McInnes its,snes->fc,fgnorm); 61e7e93795SLois Curfman McInnes } 62e7e93795SLois Curfman McInnes else { 63e7e93795SLois Curfman McInnes MPIU_printf(snes->comm, 64e7e93795SLois Curfman McInnes "iter = %d, Function value %g, Gradient norm < 1.e-11\n", 65e7e93795SLois Curfman McInnes its,snes->fc); 66e7e93795SLois Curfman McInnes } 67e7e93795SLois Curfman McInnes } else SETERRQ(1,"SNESDefaultSMonitor:Unknown method class"); 68e7e93795SLois Curfman McInnes return 0; 69e7e93795SLois Curfman McInnes } 70e7e93795SLois Curfman McInnes /* ---------------------------------------------------------------- */ 714b828684SBarry Smith /*@C 72e7e93795SLois Curfman McInnes SNESDefaultConverged - Default test for monitoring the convergence 73e7e93795SLois Curfman McInnes of the solvers for systems of nonlinear equations. 74e7e93795SLois Curfman McInnes 75e7e93795SLois Curfman McInnes Input Parameters: 76e7e93795SLois Curfman McInnes . snes - the SNES context 77e7e93795SLois Curfman McInnes . xnorm - 2-norm of current iterate 78e7e93795SLois Curfman McInnes . pnorm - 2-norm of current step 79e7e93795SLois Curfman McInnes . fnorm - 2-norm of function 80e7e93795SLois Curfman McInnes . dummy - unused context 81e7e93795SLois Curfman McInnes 82e7e93795SLois Curfman McInnes Returns: 83e7e93795SLois Curfman McInnes $ 2 if ( fnorm < atol ), 84e7e93795SLois Curfman McInnes $ 3 if ( pnorm < xtol*xnorm ), 85*5d2e0e51SBarry Smith $ 4 if ( fnorm < rtol*fnorm0 ), 86e7e93795SLois Curfman McInnes $ -2 if ( nfct > maxf ), 87e7e93795SLois Curfman McInnes $ 0 otherwise, 88e7e93795SLois Curfman McInnes 89e7e93795SLois Curfman McInnes where 90e7e93795SLois Curfman McInnes $ maxf - maximum number of function evaluations, 91e7e93795SLois Curfman McInnes $ set with SNESSetMaxFunctionEvaluations() 92e7e93795SLois Curfman McInnes $ nfct - number of function evaluations, 93e7e93795SLois Curfman McInnes $ atol - absolute function norm tolerance, 94e7e93795SLois Curfman McInnes $ set with SNESSetAbsoluteTolerance() 95e7e93795SLois Curfman McInnes $ xtol - relative function norm tolerance, 96e7e93795SLois Curfman McInnes $ set with SNESSetRelativeTolerance() 97e7e93795SLois Curfman McInnes 98e7e93795SLois Curfman McInnes .keywords: SNES, nonlinear, default, converged, convergence 99e7e93795SLois Curfman McInnes 100e7e93795SLois Curfman McInnes .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 101e7e93795SLois Curfman McInnes @*/ 102*5d2e0e51SBarry Smith int SNESDefaultConverged(SNES snes,double xnorm,double pnorm,double fnorm,void *dummy) 103e7e93795SLois Curfman McInnes { 104e7e93795SLois Curfman McInnes if (snes->method_class != SNES_NONLINEAR_EQUATIONS) SETERRQ(1, 10548d91487SBarry Smith "SNESDefaultConverged:For SNES_NONLINEAR_EQUATIONS only"); 106e7e93795SLois Curfman McInnes /* Note: Reserve return code 1, -1 for compatibility with 107e7e93795SLois Curfman McInnes SNESTrustRegionDefaultConverged */ 108*5d2e0e51SBarry Smith if (snes->iter == 1) { /* first iteration so set ttol */ 109*5d2e0e51SBarry Smith snes->ttol = fnorm*snes->rtol; 110*5d2e0e51SBarry Smith } 111*5d2e0e51SBarry Smith else { 112*5d2e0e51SBarry Smith if (fnorm <= snes->ttol) { 113*5d2e0e51SBarry Smith PLogInfo((PetscObject)snes, 114*5d2e0e51SBarry Smith "SNES:Converged due to function norm %g < %g (relative tolerance)\n",fnorm,snes->ttol); 115*5d2e0e51SBarry Smith return 4; 116*5d2e0e51SBarry Smith } 117*5d2e0e51SBarry Smith } 118*5d2e0e51SBarry Smith 119e7e93795SLois Curfman McInnes if (fnorm < snes->atol) { 120e7e93795SLois Curfman McInnes PLogInfo((PetscObject)snes, 1210de55854SLois Curfman McInnes "SNES: Converged due to function norm %g < %g\n",fnorm,snes->atol); 122e7e93795SLois Curfman McInnes return 2; 123e7e93795SLois Curfman McInnes } 124e7e93795SLois Curfman McInnes if (pnorm < snes->xtol*(xnorm)) { 125e7e93795SLois Curfman McInnes PLogInfo((PetscObject)snes, 126e7e93795SLois Curfman McInnes "SNES: Converged due to small update length: %g < %g * %g\n", 127e7e93795SLois Curfman McInnes pnorm,snes->xtol,xnorm); 128e7e93795SLois Curfman McInnes return 3; 129e7e93795SLois Curfman McInnes } 130e7e93795SLois Curfman McInnes if (snes->nfuncs > snes->max_funcs) { 131e7e93795SLois Curfman McInnes PLogInfo((PetscObject)snes, 132e7e93795SLois Curfman McInnes "SNES: Exceeded maximum number of function evaluations: %d > %d\n", 133e7e93795SLois Curfman McInnes snes->nfuncs, snes->max_funcs ); 134e7e93795SLois Curfman McInnes return -2; 135e7e93795SLois Curfman McInnes } 136e7e93795SLois Curfman McInnes return 0; 137e7e93795SLois Curfman McInnes } 138e7e93795SLois Curfman McInnes /* ------------------------------------------------------------ */ 139e7e93795SLois Curfman McInnes /*@ 140e7e93795SLois Curfman McInnes SNES_KSP_SetConvergenceTestEW - Sets alternative convergence test for 141e7e93795SLois Curfman McInnes for the linear solvers within an inexact Newton method. 142e7e93795SLois Curfman McInnes 143e7e93795SLois Curfman McInnes Input Parameter: 144e7e93795SLois Curfman McInnes . snes - SNES context 145e7e93795SLois Curfman McInnes 146e7e93795SLois Curfman McInnes Notes: 147e7e93795SLois Curfman McInnes Currently, the default is to use a constant relative tolerance for 148e7e93795SLois Curfman McInnes the inner linear solvers. Alternatively, one can use the 149e7e93795SLois Curfman McInnes Eisenstat-Walker method, where the relative convergence tolerance 150e7e93795SLois Curfman McInnes is reset at each Newton iteration according progress of the nonlinear 151e7e93795SLois Curfman McInnes solver. 152e7e93795SLois Curfman McInnes 153e7e93795SLois Curfman McInnes Reference: 154e7e93795SLois Curfman McInnes S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 155e7e93795SLois Curfman McInnes inexact Newton method", Utah State University Math. Stat. Dept. Res. 156e7e93795SLois Curfman McInnes Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 157e7e93795SLois Curfman McInnes 158e7e93795SLois Curfman McInnes .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 159e7e93795SLois Curfman McInnes @*/ 160e7e93795SLois Curfman McInnes int SNES_KSP_SetConvergenceTestEW(SNES snes) 161e7e93795SLois Curfman McInnes { 162e7e93795SLois Curfman McInnes snes->ksp_ewconv = 1; 163e7e93795SLois Curfman McInnes return 0; 164e7e93795SLois Curfman McInnes } 165e7e93795SLois Curfman McInnes 166e7e93795SLois Curfman McInnes /*@ 167e7e93795SLois Curfman McInnes SNES_KSP_SetParametersEW - Sets parameters for Eisenstat-Walker 168e7e93795SLois Curfman McInnes convergence criteria for the linear solvers within an inexact 169e7e93795SLois Curfman McInnes Newton method. 170e7e93795SLois Curfman McInnes 171e7e93795SLois Curfman McInnes Input Parameters: 172e7e93795SLois Curfman McInnes . snes - SNES context 173e7e93795SLois Curfman McInnes . version - version 1 or 2 (default is 2) 174e7e93795SLois Curfman McInnes . rtol_0 - initial relative tolerance 175e7e93795SLois Curfman McInnes $ (0 <= rtol_0 < 1) 176e7e93795SLois Curfman McInnes . rtol_max - maximum relative tolerance 177e7e93795SLois Curfman McInnes $ (0 <= rtol_max < 1) 178e7e93795SLois Curfman McInnes . alpha - power for version 2 rtol computation 179e7e93795SLois Curfman McInnes $ (1 < alpha <= 2) 180e7e93795SLois Curfman McInnes . alpha2 - power for safeguard 181e7e93795SLois Curfman McInnes . gamma2 - multiplicative factor for version 2 rtol computation 182e7e93795SLois Curfman McInnes $ (0 <= gamma2 <= 1) 183e7e93795SLois Curfman McInnes . threshold - threshold for imposing safeguard 184e7e93795SLois Curfman McInnes $ (0 < threshold < 1) 185e7e93795SLois Curfman McInnes 186e7e93795SLois Curfman McInnes Note: 187e7e93795SLois Curfman McInnes Use PETSC_DEFAULT to retain the default for any of the parameters. 188e7e93795SLois Curfman McInnes 189e7e93795SLois Curfman McInnes Reference: 190e7e93795SLois Curfman McInnes S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 191e7e93795SLois Curfman McInnes inexact Newton method", Utah State University Math. Stat. Dept. Res. 192e7e93795SLois Curfman McInnes Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 193e7e93795SLois Curfman McInnes 194e7e93795SLois Curfman McInnes .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 195e7e93795SLois Curfman McInnes 196e7e93795SLois Curfman McInnes .seealso: SNES_KSP_SetConvergenceTestEW() 197e7e93795SLois Curfman McInnes @*/ 198e7e93795SLois Curfman McInnes int SNES_KSP_SetParametersEW(SNES snes,int version,double rtol_0, 199e7e93795SLois Curfman McInnes double rtol_max,double gamma2,double alpha, 200e7e93795SLois Curfman McInnes double alpha2,double threshold) 201e7e93795SLois Curfman McInnes { 202e7e93795SLois Curfman McInnes SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx; 203e7e93795SLois Curfman McInnes if (!kctx) SETERRQ(1,"SNES_KSP_SetParametersEW:No context"); 204e7e93795SLois Curfman McInnes if (version != PETSC_DEFAULT) kctx->version = version; 205e7e93795SLois Curfman McInnes if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 206e7e93795SLois Curfman McInnes if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 207e7e93795SLois Curfman McInnes if (gamma2 != PETSC_DEFAULT) kctx->gamma = gamma2; 208e7e93795SLois Curfman McInnes if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 209e7e93795SLois Curfman McInnes if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 210e7e93795SLois Curfman McInnes if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 211e7e93795SLois Curfman McInnes if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) SETERRQ(1, 212e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: 0.0 <= rtol_0 < 1.0\n"); 213e7e93795SLois Curfman McInnes if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) SETERRQ(1, 214e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: 0.0 <= rtol_max < 1.0\n"); 215e7e93795SLois Curfman McInnes if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) SETERRQ(1, 216e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: 0.0 < threshold < 1.0\n"); 217e7e93795SLois Curfman McInnes if (kctx->gamma < 0.0 || kctx->gamma > 1.0) SETERRQ(1, 218e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: 0.0 <= alpha <= 1.0\n"); 219e7e93795SLois Curfman McInnes if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) SETERRQ(1, 220e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: 1.0 < alpha <= 2.0\n"); 221e7e93795SLois Curfman McInnes if (kctx->version != 1 && kctx->version !=2) SETERRQ(1, 222e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: Only versions 1 and 2 are supported"); 223e7e93795SLois Curfman McInnes return 0; 224e7e93795SLois Curfman McInnes } 225e7e93795SLois Curfman McInnes 226e7e93795SLois Curfman McInnes int SNES_KSP_EW_ComputeRelativeTolerance_Private(SNES snes,KSP ksp) 227e7e93795SLois Curfman McInnes { 228e7e93795SLois Curfman McInnes SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx; 229e7e93795SLois Curfman McInnes double rtol, stol; 230e7e93795SLois Curfman McInnes int ierr; 231e7e93795SLois Curfman McInnes if (!kctx) 232e7e93795SLois Curfman McInnes SETERRQ(1,"SNES_KSP_EW_ComputeRelativeTolerance_Private:No context"); 233e7e93795SLois Curfman McInnes if (snes->iter == 1) { 234e7e93795SLois Curfman McInnes rtol = kctx->rtol_0; 235e7e93795SLois Curfman McInnes } else { 236e7e93795SLois Curfman McInnes if (kctx->version == 1) { 237e7e93795SLois Curfman McInnes rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 238e7e93795SLois Curfman McInnes if (rtol < 0.0) rtol = -rtol; 239e7e93795SLois Curfman McInnes stol = pow(kctx->rtol_last,kctx->alpha2); 2400452661fSBarry Smith if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 241e7e93795SLois Curfman McInnes } else if (kctx->version == 2) { 242e7e93795SLois Curfman McInnes rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 243e7e93795SLois Curfman McInnes stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 2440452661fSBarry Smith if (stol > kctx->threshold) rtol = PetscMax(rtol,stol); 245e7e93795SLois Curfman McInnes } else SETERRQ(1, 24648d91487SBarry Smith "SNES_KSP_EW_Converged_Private:Only versions 1 or 2 are supported"); 247e7e93795SLois Curfman McInnes } 2480452661fSBarry Smith rtol = PetscMin(rtol,kctx->rtol_max); 249e7e93795SLois Curfman McInnes kctx->rtol_last = rtol; 250e7e93795SLois Curfman McInnes PLogInfo((PetscObject)snes, 251e7e93795SLois Curfman McInnes "SNES: iter %d, Eisenstat-Walker (version %d) KSP rtol = %g\n", 252e7e93795SLois Curfman McInnes snes->iter,kctx->version,rtol); 253e7e93795SLois Curfman McInnes ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT); 254e7e93795SLois Curfman McInnes CHKERRQ(ierr); 255e7e93795SLois Curfman McInnes kctx->norm_last = snes->norm; 256e7e93795SLois Curfman McInnes return 0; 257e7e93795SLois Curfman McInnes } 258e7e93795SLois Curfman McInnes 259e7e93795SLois Curfman McInnes int SNES_KSP_EW_Converged_Private(KSP ksp,int n,double rnorm,void *ctx) 260e7e93795SLois Curfman McInnes { 261e7e93795SLois Curfman McInnes SNES snes = (SNES)ctx; 262e7e93795SLois Curfman McInnes SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx; 263e7e93795SLois Curfman McInnes int convinfo; 264e7e93795SLois Curfman McInnes 26548d91487SBarry Smith if (!kctx) SETERRQ(1,"SNES_KSP_EW_Converged_Private:No convergence context"); 266e7e93795SLois Curfman McInnes if (n == 0) SNES_KSP_EW_ComputeRelativeTolerance_Private(snes,ksp); 267e7e93795SLois Curfman McInnes convinfo = KSPDefaultConverged(ksp,n,rnorm,ctx); 268e7e93795SLois Curfman McInnes kctx->lresid_last = rnorm; 269e7e93795SLois Curfman McInnes if (convinfo) 270e7e93795SLois Curfman McInnes PLogInfo((PetscObject)snes,"SNES: KSP iterations=%d, rnorm=%g\n",n,rnorm); 271e7e93795SLois Curfman McInnes return convinfo; 272e7e93795SLois Curfman McInnes } 273e7e93795SLois Curfman McInnes 274e7e93795SLois Curfman McInnes 275