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