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