1e7e93795SLois Curfman McInnes #ifndef lint 2*4b828684SBarry Smith static char vcid[] = "$Id: snesut.c,v 1.1 1995/08/22 19:40:03 curfman 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 #if defined(HAVE_STRING_H) 8e7e93795SLois Curfman McInnes #include <string.h> 9e7e93795SLois Curfman McInnes #endif 10e7e93795SLois Curfman McInnes 11*4b828684SBarry 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 /* ---------------------------------------------------------------- */ 74*4b828684SBarry 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, 113e7e93795SLois Curfman McInnes "SNES: Converged due to function norm %g < %g\n", 114e7e93795SLois Curfman McInnes fnorm,snes->atol); 115e7e93795SLois Curfman McInnes return 2; 116e7e93795SLois Curfman McInnes } 117e7e93795SLois Curfman McInnes if (pnorm < snes->xtol*(xnorm)) { 118e7e93795SLois Curfman McInnes PLogInfo((PetscObject)snes, 119e7e93795SLois Curfman McInnes "SNES: Converged due to small update length: %g < %g * %g\n", 120e7e93795SLois Curfman McInnes pnorm,snes->xtol,xnorm); 121e7e93795SLois Curfman McInnes return 3; 122e7e93795SLois Curfman McInnes } 123e7e93795SLois Curfman McInnes if (snes->nfuncs > snes->max_funcs) { 124e7e93795SLois Curfman McInnes PLogInfo((PetscObject)snes, 125e7e93795SLois Curfman McInnes "SNES: Exceeded maximum number of function evaluations: %d > %d\n", 126e7e93795SLois Curfman McInnes snes->nfuncs, snes->max_funcs ); 127e7e93795SLois Curfman McInnes return -2; 128e7e93795SLois Curfman McInnes } 129e7e93795SLois Curfman McInnes return 0; 130e7e93795SLois Curfman McInnes } 131e7e93795SLois Curfman McInnes /* ------------------------------------------------------------ */ 132e7e93795SLois Curfman McInnes /*@ 133e7e93795SLois Curfman McInnes SNES_KSP_SetConvergenceTestEW - Sets alternative convergence test for 134e7e93795SLois Curfman McInnes for the linear solvers within an inexact Newton method. 135e7e93795SLois Curfman McInnes 136e7e93795SLois Curfman McInnes Input Parameter: 137e7e93795SLois Curfman McInnes . snes - SNES context 138e7e93795SLois Curfman McInnes 139e7e93795SLois Curfman McInnes Notes: 140e7e93795SLois Curfman McInnes Currently, the default is to use a constant relative tolerance for 141e7e93795SLois Curfman McInnes the inner linear solvers. Alternatively, one can use the 142e7e93795SLois Curfman McInnes Eisenstat-Walker method, where the relative convergence tolerance 143e7e93795SLois Curfman McInnes is reset at each Newton iteration according progress of the nonlinear 144e7e93795SLois Curfman McInnes solver. 145e7e93795SLois Curfman McInnes 146e7e93795SLois Curfman McInnes Reference: 147e7e93795SLois Curfman McInnes S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 148e7e93795SLois Curfman McInnes inexact Newton method", Utah State University Math. Stat. Dept. Res. 149e7e93795SLois Curfman McInnes Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 150e7e93795SLois Curfman McInnes 151e7e93795SLois Curfman McInnes .keywords: SNES, KSP, Eisenstat, Walker, convergence, test, inexact, Newton 152e7e93795SLois Curfman McInnes @*/ 153e7e93795SLois Curfman McInnes int SNES_KSP_SetConvergenceTestEW(SNES snes) 154e7e93795SLois Curfman McInnes { 155e7e93795SLois Curfman McInnes snes->ksp_ewconv = 1; 156e7e93795SLois Curfman McInnes return 0; 157e7e93795SLois Curfman McInnes } 158e7e93795SLois Curfman McInnes 159e7e93795SLois Curfman McInnes /*@ 160e7e93795SLois Curfman McInnes SNES_KSP_SetParametersEW - Sets parameters for Eisenstat-Walker 161e7e93795SLois Curfman McInnes convergence criteria for the linear solvers within an inexact 162e7e93795SLois Curfman McInnes Newton method. 163e7e93795SLois Curfman McInnes 164e7e93795SLois Curfman McInnes Input Parameters: 165e7e93795SLois Curfman McInnes . snes - SNES context 166e7e93795SLois Curfman McInnes . version - version 1 or 2 (default is 2) 167e7e93795SLois Curfman McInnes . rtol_0 - initial relative tolerance 168e7e93795SLois Curfman McInnes $ (0 <= rtol_0 < 1) 169e7e93795SLois Curfman McInnes . rtol_max - maximum relative tolerance 170e7e93795SLois Curfman McInnes $ (0 <= rtol_max < 1) 171e7e93795SLois Curfman McInnes . alpha - power for version 2 rtol computation 172e7e93795SLois Curfman McInnes $ (1 < alpha <= 2) 173e7e93795SLois Curfman McInnes . alpha2 - power for safeguard 174e7e93795SLois Curfman McInnes . gamma2 - multiplicative factor for version 2 rtol computation 175e7e93795SLois Curfman McInnes $ (0 <= gamma2 <= 1) 176e7e93795SLois Curfman McInnes . threshold - threshold for imposing safeguard 177e7e93795SLois Curfman McInnes $ (0 < threshold < 1) 178e7e93795SLois Curfman McInnes 179e7e93795SLois Curfman McInnes Note: 180e7e93795SLois Curfman McInnes Use PETSC_DEFAULT to retain the default for any of the parameters. 181e7e93795SLois Curfman McInnes 182e7e93795SLois Curfman McInnes Reference: 183e7e93795SLois Curfman McInnes S. C. Eisenstat and H. F. Walker, "Choosing the forcing terms in an 184e7e93795SLois Curfman McInnes inexact Newton method", Utah State University Math. Stat. Dept. Res. 185e7e93795SLois Curfman McInnes Report 6/94/75, June, 1994, to appear in SIAM J. Sci. Comput. 186e7e93795SLois Curfman McInnes 187e7e93795SLois Curfman McInnes .keywords: SNES, KSP, Eisenstat, Walker, set, parameters 188e7e93795SLois Curfman McInnes 189e7e93795SLois Curfman McInnes .seealso: SNES_KSP_SetConvergenceTestEW() 190e7e93795SLois Curfman McInnes @*/ 191e7e93795SLois Curfman McInnes int SNES_KSP_SetParametersEW(SNES snes,int version,double rtol_0, 192e7e93795SLois Curfman McInnes double rtol_max,double gamma2,double alpha, 193e7e93795SLois Curfman McInnes double alpha2,double threshold) 194e7e93795SLois Curfman McInnes { 195e7e93795SLois Curfman McInnes SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx; 196e7e93795SLois Curfman McInnes if (!kctx) SETERRQ(1,"SNES_KSP_SetParametersEW: No context"); 197e7e93795SLois Curfman McInnes if (version != PETSC_DEFAULT) kctx->version = version; 198e7e93795SLois Curfman McInnes if (rtol_0 != PETSC_DEFAULT) kctx->rtol_0 = rtol_0; 199e7e93795SLois Curfman McInnes if (rtol_max != PETSC_DEFAULT) kctx->rtol_max = rtol_max; 200e7e93795SLois Curfman McInnes if (gamma2 != PETSC_DEFAULT) kctx->gamma = gamma2; 201e7e93795SLois Curfman McInnes if (alpha != PETSC_DEFAULT) kctx->alpha = alpha; 202e7e93795SLois Curfman McInnes if (alpha2 != PETSC_DEFAULT) kctx->alpha2 = alpha2; 203e7e93795SLois Curfman McInnes if (threshold != PETSC_DEFAULT) kctx->threshold = threshold; 204e7e93795SLois Curfman McInnes if (kctx->rtol_0 < 0.0 || kctx->rtol_0 >= 1.0) SETERRQ(1, 205e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: 0.0 <= rtol_0 < 1.0\n"); 206e7e93795SLois Curfman McInnes if (kctx->rtol_max < 0.0 || kctx->rtol_max >= 1.0) SETERRQ(1, 207e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: 0.0 <= rtol_max < 1.0\n"); 208e7e93795SLois Curfman McInnes if (kctx->threshold <= 0.0 || kctx->threshold >= 1.0) SETERRQ(1, 209e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: 0.0 < threshold < 1.0\n"); 210e7e93795SLois Curfman McInnes if (kctx->gamma < 0.0 || kctx->gamma > 1.0) SETERRQ(1, 211e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: 0.0 <= alpha <= 1.0\n"); 212e7e93795SLois Curfman McInnes if (kctx->alpha <= 1.0 || kctx->alpha > 2.0) SETERRQ(1, 213e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: 1.0 < alpha <= 2.0\n"); 214e7e93795SLois Curfman McInnes if (kctx->version != 1 && kctx->version !=2) SETERRQ(1, 215e7e93795SLois Curfman McInnes "SNES_KSP_SetParametersEW: Only versions 1 and 2 are supported"); 216e7e93795SLois Curfman McInnes return 0; 217e7e93795SLois Curfman McInnes } 218e7e93795SLois Curfman McInnes 219e7e93795SLois Curfman McInnes int SNES_KSP_EW_ComputeRelativeTolerance_Private(SNES snes,KSP ksp) 220e7e93795SLois Curfman McInnes { 221e7e93795SLois Curfman McInnes SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx; 222e7e93795SLois Curfman McInnes double rtol, stol; 223e7e93795SLois Curfman McInnes int ierr; 224e7e93795SLois Curfman McInnes if (!kctx) 225e7e93795SLois Curfman McInnes SETERRQ(1,"SNES_KSP_EW_ComputeRelativeTolerance_Private: No context"); 226e7e93795SLois Curfman McInnes if (snes->iter == 1) { 227e7e93795SLois Curfman McInnes rtol = kctx->rtol_0; 228e7e93795SLois Curfman McInnes } else { 229e7e93795SLois Curfman McInnes if (kctx->version == 1) { 230e7e93795SLois Curfman McInnes rtol = (snes->norm - kctx->lresid_last)/kctx->norm_last; 231e7e93795SLois Curfman McInnes if (rtol < 0.0) rtol = -rtol; 232e7e93795SLois Curfman McInnes stol = pow(kctx->rtol_last,kctx->alpha2); 233e7e93795SLois Curfman McInnes if (stol > kctx->threshold) rtol = PETSCMAX(rtol,stol); 234e7e93795SLois Curfman McInnes } else if (kctx->version == 2) { 235e7e93795SLois Curfman McInnes rtol = kctx->gamma * pow(snes->norm/kctx->norm_last,kctx->alpha); 236e7e93795SLois Curfman McInnes stol = kctx->gamma * pow(kctx->rtol_last,kctx->alpha); 237e7e93795SLois Curfman McInnes if (stol > kctx->threshold) rtol = PETSCMAX(rtol,stol); 238e7e93795SLois Curfman McInnes } else SETERRQ(1, 239e7e93795SLois Curfman McInnes "SNES_KSP_EW_Converged_Private: Only versions 1 and 2 are supported"); 240e7e93795SLois Curfman McInnes } 241e7e93795SLois Curfman McInnes rtol = PETSCMIN(rtol,kctx->rtol_max); 242e7e93795SLois Curfman McInnes kctx->rtol_last = rtol; 243e7e93795SLois Curfman McInnes PLogInfo((PetscObject)snes, 244e7e93795SLois Curfman McInnes "SNES: iter %d, Eisenstat-Walker (version %d) KSP rtol = %g\n", 245e7e93795SLois Curfman McInnes snes->iter,kctx->version,rtol); 246e7e93795SLois Curfman McInnes ierr = KSPSetTolerances(ksp,rtol,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT); 247e7e93795SLois Curfman McInnes CHKERRQ(ierr); 248e7e93795SLois Curfman McInnes kctx->norm_last = snes->norm; 249e7e93795SLois Curfman McInnes return 0; 250e7e93795SLois Curfman McInnes } 251e7e93795SLois Curfman McInnes 252e7e93795SLois Curfman McInnes int SNES_KSP_EW_Converged_Private(KSP ksp,int n,double rnorm,void *ctx) 253e7e93795SLois Curfman McInnes { 254e7e93795SLois Curfman McInnes SNES snes = (SNES)ctx; 255e7e93795SLois Curfman McInnes SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx; 256e7e93795SLois Curfman McInnes int convinfo; 257e7e93795SLois Curfman McInnes 258e7e93795SLois Curfman McInnes if (!kctx) SETERRQ(1, 259e7e93795SLois Curfman McInnes "SNES_KSP_EW_Converged_Private: Convergence context does not exist"); 260e7e93795SLois Curfman McInnes if (n == 0) SNES_KSP_EW_ComputeRelativeTolerance_Private(snes,ksp); 261e7e93795SLois Curfman McInnes convinfo = KSPDefaultConverged(ksp,n,rnorm,ctx); 262e7e93795SLois Curfman McInnes kctx->lresid_last = rnorm; 263e7e93795SLois Curfman McInnes if (convinfo) 264e7e93795SLois Curfman McInnes PLogInfo((PetscObject)snes,"SNES: KSP iterations=%d, rnorm=%g\n",n,rnorm); 265e7e93795SLois Curfman McInnes return convinfo; 266e7e93795SLois Curfman McInnes } 267e7e93795SLois Curfman McInnes 268e7e93795SLois Curfman McInnes 269