1 2 #include <../src/snes/impls/tr/trimpl.h> /*I "petscsnes.h" I*/ 3 4 typedef struct { 5 SNES snes; 6 /* Information on the regular SNES convergence test; which may have been user provided */ 7 PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*); 8 PetscErrorCode (*convdestroy)(void*); 9 void *convctx; 10 } SNES_TR_KSPConverged_Ctx; 11 12 static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx) 13 { 14 SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx; 15 SNES snes = ctx->snes; 16 SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data; 17 Vec x; 18 PetscReal nrm; 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 ierr = (*ctx->convtest)(ksp,n,rnorm,reason,ctx->convctx);CHKERRQ(ierr); 23 if (*reason) { 24 ierr = PetscInfo2(snes,"Default or user provided convergence test KSP iterations=%D, rnorm=%g\n",n,(double)rnorm);CHKERRQ(ierr); 25 } 26 /* Determine norm of solution */ 27 ierr = KSPBuildSolution(ksp,0,&x);CHKERRQ(ierr); 28 ierr = VecNorm(x,NORM_2,&nrm);CHKERRQ(ierr); 29 if (nrm >= neP->delta) { 30 ierr = PetscInfo2(snes,"Ending linear iteration early, delta=%g, length=%g\n",(double)neP->delta,(double)nrm);CHKERRQ(ierr); 31 *reason = KSP_CONVERGED_STEP_LENGTH; 32 } 33 PetscFunctionReturn(0); 34 } 35 36 static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx) 37 { 38 SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx; 39 PetscErrorCode ierr; 40 41 PetscFunctionBegin; 42 ierr = (*ctx->convdestroy)(ctx->convctx);CHKERRQ(ierr); 43 ierr = PetscFree(ctx);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 /* ---------------------------------------------------------------- */ 48 /* 49 SNESTR_Converged_Private -test convergence JUST for 50 the trust region tolerance. 51 52 */ 53 static PetscErrorCode SNESTR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) 54 { 55 SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 *reason = SNES_CONVERGED_ITERATING; 60 if (neP->delta < xnorm * snes->deltatol) { 61 ierr = PetscInfo3(snes,"Converged due to trust region param %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol);CHKERRQ(ierr); 62 *reason = SNES_CONVERGED_TR_DELTA; 63 } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) { 64 ierr = PetscInfo1(snes,"Exceeded maximum number of function evaluations: %D\n",snes->max_funcs);CHKERRQ(ierr); 65 *reason = SNES_DIVERGED_FUNCTION_COUNT; 66 } 67 PetscFunctionReturn(0); 68 } 69 70 71 /* 72 SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust 73 region approach for solving systems of nonlinear equations. 74 75 76 */ 77 static PetscErrorCode SNESSolve_NEWTONTR(SNES snes) 78 { 79 SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data; 80 Vec X,F,Y,G,Ytmp; 81 PetscErrorCode ierr; 82 PetscInt maxits,i,lits; 83 PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1; 84 PetscScalar cnorm; 85 KSP ksp; 86 SNESConvergedReason reason = SNES_CONVERGED_ITERATING; 87 PetscBool breakout = PETSC_FALSE; 88 SNES_TR_KSPConverged_Ctx *ctx; 89 PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*); 90 91 PetscFunctionBegin; 92 if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name); 93 94 maxits = snes->max_its; /* maximum number of iterations */ 95 X = snes->vec_sol; /* solution vector */ 96 F = snes->vec_func; /* residual vector */ 97 Y = snes->work[0]; /* work vectors */ 98 G = snes->work[1]; 99 Ytmp = snes->work[2]; 100 101 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 102 snes->iter = 0; 103 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 104 105 /* Set the linear stopping criteria to use the More' trick. */ 106 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 107 ierr = KSPGetConvergenceTest(ksp,&convtest,NULL,NULL);CHKERRQ(ierr); 108 if (convtest != SNESTR_KSPConverged_Private) { 109 ierr = PetscNew(&ctx);CHKERRQ(ierr); 110 ctx->snes = snes; 111 ierr = KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);CHKERRQ(ierr); 112 ierr = KSPSetConvergenceTest(ksp,SNESTR_KSPConverged_Private,ctx,SNESTR_KSPConverged_Destroy);CHKERRQ(ierr); 113 ierr = PetscInfo(snes,"Using Krylov convergence test SNESTR_KSPConverged_Private\n");CHKERRQ(ierr); 114 } 115 116 if (!snes->vec_func_init_set) { 117 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); /* F(X) */ 118 } else snes->vec_func_init_set = PETSC_FALSE; 119 120 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- || F || */ 121 SNESCheckFunctionNorm(snes,fnorm); 122 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* fnorm <- || F || */ 123 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 124 snes->norm = fnorm; 125 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 126 delta = xnorm ? neP->delta0*xnorm : neP->delta0; 127 neP->delta = delta; 128 ierr = SNESLogConvergenceHistory(snes,fnorm,0);CHKERRQ(ierr); 129 ierr = SNESMonitor(snes,0,fnorm);CHKERRQ(ierr); 130 131 /* test convergence */ 132 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 133 if (snes->reason) PetscFunctionReturn(0); 134 135 136 for (i=0; i<maxits; i++) { 137 138 /* Call general purpose update function */ 139 if (snes->ops->update) { 140 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 141 } 142 143 /* Solve J Y = F, where J is Jacobian matrix */ 144 ierr = SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 145 SNESCheckJacobianDomainerror(snes); 146 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);CHKERRQ(ierr); 147 ierr = KSPSolve(snes->ksp,F,Ytmp);CHKERRQ(ierr); 148 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 149 150 snes->linear_its += lits; 151 152 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 153 ierr = VecNorm(Ytmp,NORM_2,&nrm);CHKERRQ(ierr); 154 norm1 = nrm; 155 while (1) { 156 ierr = VecCopy(Ytmp,Y);CHKERRQ(ierr); 157 nrm = norm1; 158 159 /* Scale Y if need be and predict new value of F norm */ 160 if (nrm >= delta) { 161 nrm = delta/nrm; 162 gpnorm = (1.0 - nrm)*fnorm; 163 cnorm = nrm; 164 ierr = PetscInfo1(snes,"Scaling direction by %g\n",(double)nrm);CHKERRQ(ierr); 165 ierr = VecScale(Y,cnorm);CHKERRQ(ierr); 166 nrm = gpnorm; 167 ynorm = delta; 168 } else { 169 gpnorm = 0.0; 170 ierr = PetscInfo(snes,"Direction is in Trust Region\n");CHKERRQ(ierr); 171 ynorm = nrm; 172 } 173 ierr = VecCopy(Y,snes->vec_sol_update);CHKERRQ(ierr); 174 ierr = VecAYPX(Y,-1.0,X);CHKERRQ(ierr); /* Y <- X - Y */ 175 ierr = SNESComputeFunction(snes,Y,G);CHKERRQ(ierr); /* F(X) */ 176 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); /* gnorm <- || g || */ 177 SNESCheckFunctionNorm(snes,gnorm); 178 if (fnorm == gpnorm) rho = 0.0; 179 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 180 181 /* Update size of trust region */ 182 if (rho < neP->mu) delta *= neP->delta1; 183 else if (rho < neP->eta) delta *= neP->delta2; 184 else delta *= neP->delta3; 185 ierr = PetscInfo3(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm);CHKERRQ(ierr); 186 ierr = PetscInfo3(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta);CHKERRQ(ierr); 187 188 neP->delta = delta; 189 if (rho > neP->sigma) break; 190 ierr = PetscInfo(snes,"Trying again in smaller region\n");CHKERRQ(ierr); 191 /* check to see if progress is hopeless */ 192 neP->itflag = PETSC_FALSE; 193 ierr = SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 194 if (!reason) { ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); } 195 if (reason) { 196 /* We're not progressing, so return with the current iterate */ 197 ierr = SNESMonitor(snes,i+1,fnorm);CHKERRQ(ierr); 198 breakout = PETSC_TRUE; 199 break; 200 } 201 snes->numFailures++; 202 } 203 if (!breakout) { 204 /* Update function and solution vectors */ 205 fnorm = gnorm; 206 ierr = VecCopy(G,F);CHKERRQ(ierr); 207 ierr = VecCopy(Y,X);CHKERRQ(ierr); 208 /* Monitor convergence */ 209 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 210 snes->iter = i+1; 211 snes->norm = fnorm; 212 snes->xnorm = xnorm; 213 snes->ynorm = ynorm; 214 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 215 ierr = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr); 216 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 217 /* Test for convergence, xnorm = || X || */ 218 neP->itflag = PETSC_TRUE; 219 if (snes->ops->converged != SNESConvergedSkip) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 220 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 221 if (reason) break; 222 } else break; 223 } 224 if (i == maxits) { 225 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 226 if (!reason) reason = SNES_DIVERGED_MAX_IT; 227 } 228 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 229 snes->reason = reason; 230 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 231 PetscFunctionReturn(0); 232 } 233 /*------------------------------------------------------------*/ 234 static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes) 235 { 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr); 240 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 241 PetscFunctionReturn(0); 242 } 243 244 PetscErrorCode SNESReset_NEWTONTR(SNES snes) 245 { 246 247 PetscFunctionBegin; 248 PetscFunctionReturn(0); 249 } 250 251 static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes) 252 { 253 PetscErrorCode ierr; 254 255 PetscFunctionBegin; 256 ierr = SNESReset_NEWTONTR(snes);CHKERRQ(ierr); 257 ierr = PetscFree(snes->data);CHKERRQ(ierr); 258 PetscFunctionReturn(0); 259 } 260 /*------------------------------------------------------------*/ 261 262 static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes) 263 { 264 SNES_NEWTONTR *ctx = (SNES_NEWTONTR*)snes->data; 265 PetscErrorCode ierr; 266 267 PetscFunctionBegin; 268 ierr = PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");CHKERRQ(ierr); 269 ierr = PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);CHKERRQ(ierr); 270 ierr = PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL);CHKERRQ(ierr); 271 ierr = PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL);CHKERRQ(ierr); 272 ierr = PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL);CHKERRQ(ierr); 273 ierr = PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);CHKERRQ(ierr); 274 ierr = PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL);CHKERRQ(ierr); 275 ierr = PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL);CHKERRQ(ierr); 276 ierr = PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL);CHKERRQ(ierr); 277 ierr = PetscOptionsTail();CHKERRQ(ierr); 278 PetscFunctionReturn(0); 279 } 280 281 static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer) 282 { 283 SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 284 PetscErrorCode ierr; 285 PetscBool iascii; 286 287 PetscFunctionBegin; 288 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 289 if (iascii) { 290 ierr = PetscViewerASCIIPrintf(viewer," Trust region tolerance (-snes_trtol)\n",(double)snes->deltatol);CHKERRQ(ierr); 291 ierr = PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma);CHKERRQ(ierr); 292 ierr = PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3);CHKERRQ(ierr); 293 } 294 PetscFunctionReturn(0); 295 } 296 /* ------------------------------------------------------------ */ 297 /*MC 298 SNESNEWTONTR - Newton based nonlinear solver that uses a trust region 299 300 Options Database: 301 + -snes_trtol <tol> - trust region tolerance 302 . -snes_tr_mu <mu> - trust region parameter 303 . -snes_tr_eta <eta> - trust region parameter 304 . -snes_tr_sigma <sigma> - trust region parameter 305 . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x) 306 . -snes_tr_delta1 <delta1> - trust region parameter 307 . -snes_tr_delta2 <delta2> - trust region parameter 308 - -snes_tr_delta3 <delta3> - trust region parameter 309 310 The basic algorithm is taken from "The Minpack Project", by More', 311 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 312 of Mathematical Software", Wayne Cowell, editor. 313 314 Level: intermediate 315 316 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance() 317 318 M*/ 319 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes) 320 { 321 SNES_NEWTONTR *neP; 322 PetscErrorCode ierr; 323 324 PetscFunctionBegin; 325 snes->ops->setup = SNESSetUp_NEWTONTR; 326 snes->ops->solve = SNESSolve_NEWTONTR; 327 snes->ops->destroy = SNESDestroy_NEWTONTR; 328 snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR; 329 snes->ops->view = SNESView_NEWTONTR; 330 snes->ops->reset = SNESReset_NEWTONTR; 331 332 snes->usesksp = PETSC_TRUE; 333 snes->usesnpc = PETSC_FALSE; 334 335 snes->alwayscomputesfinalresidual = PETSC_TRUE; 336 337 ierr = PetscNewLog(snes,&neP);CHKERRQ(ierr); 338 snes->data = (void*)neP; 339 neP->mu = 0.25; 340 neP->eta = 0.75; 341 neP->delta = 0.0; 342 neP->delta0 = 0.2; 343 neP->delta1 = 0.3; 344 neP->delta2 = 0.75; 345 neP->delta3 = 2.0; 346 neP->sigma = 0.0001; 347 neP->itflag = PETSC_FALSE; 348 neP->rnorm0 = 0.0; 349 neP->ttol = 0.0; 350 PetscFunctionReturn(0); 351 } 352 353