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 if (fnorm == gpnorm) rho = 0.0; 178 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 179 180 /* Update size of trust region */ 181 if (rho < neP->mu) delta *= neP->delta1; 182 else if (rho < neP->eta) delta *= neP->delta2; 183 else delta *= neP->delta3; 184 ierr = PetscInfo3(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm);CHKERRQ(ierr); 185 ierr = PetscInfo3(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta);CHKERRQ(ierr); 186 187 neP->delta = delta; 188 if (rho > neP->sigma) break; 189 ierr = PetscInfo(snes,"Trying again in smaller region\n");CHKERRQ(ierr); 190 /* check to see if progress is hopeless */ 191 neP->itflag = PETSC_FALSE; 192 ierr = SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 193 if (!reason) { ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); } 194 if (reason) { 195 /* We're not progressing, so return with the current iterate */ 196 ierr = SNESMonitor(snes,i+1,fnorm);CHKERRQ(ierr); 197 breakout = PETSC_TRUE; 198 break; 199 } 200 snes->numFailures++; 201 } 202 if (!breakout) { 203 /* Update function and solution vectors */ 204 fnorm = gnorm; 205 ierr = VecCopy(G,F);CHKERRQ(ierr); 206 ierr = VecCopy(Y,X);CHKERRQ(ierr); 207 /* Monitor convergence */ 208 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 209 snes->iter = i+1; 210 snes->norm = fnorm; 211 snes->xnorm = xnorm; 212 snes->ynorm = ynorm; 213 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 214 ierr = SNESLogConvergenceHistory(snes,snes->norm,lits);CHKERRQ(ierr); 215 ierr = SNESMonitor(snes,snes->iter,snes->norm);CHKERRQ(ierr); 216 /* Test for convergence, xnorm = || X || */ 217 neP->itflag = PETSC_TRUE; 218 if (snes->ops->converged != SNESConvergedSkip) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 219 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 220 if (reason) break; 221 } else break; 222 } 223 if (i == maxits) { 224 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 225 if (!reason) reason = SNES_DIVERGED_MAX_IT; 226 } 227 ierr = PetscObjectSAWsTakeAccess((PetscObject)snes);CHKERRQ(ierr); 228 snes->reason = reason; 229 ierr = PetscObjectSAWsGrantAccess((PetscObject)snes);CHKERRQ(ierr); 230 PetscFunctionReturn(0); 231 } 232 /*------------------------------------------------------------*/ 233 static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes) 234 { 235 PetscErrorCode ierr; 236 237 PetscFunctionBegin; 238 ierr = SNESSetWorkVecs(snes,3);CHKERRQ(ierr); 239 ierr = SNESSetUpMatrices(snes);CHKERRQ(ierr); 240 PetscFunctionReturn(0); 241 } 242 243 PetscErrorCode SNESReset_NEWTONTR(SNES snes) 244 { 245 246 PetscFunctionBegin; 247 PetscFunctionReturn(0); 248 } 249 250 static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes) 251 { 252 PetscErrorCode ierr; 253 254 PetscFunctionBegin; 255 ierr = SNESReset_NEWTONTR(snes);CHKERRQ(ierr); 256 ierr = PetscFree(snes->data);CHKERRQ(ierr); 257 PetscFunctionReturn(0); 258 } 259 /*------------------------------------------------------------*/ 260 261 static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes) 262 { 263 SNES_NEWTONTR *ctx = (SNES_NEWTONTR*)snes->data; 264 PetscErrorCode ierr; 265 266 PetscFunctionBegin; 267 ierr = PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");CHKERRQ(ierr); 268 ierr = PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);CHKERRQ(ierr); 269 ierr = PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL);CHKERRQ(ierr); 270 ierr = PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL);CHKERRQ(ierr); 271 ierr = PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL);CHKERRQ(ierr); 272 ierr = PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);CHKERRQ(ierr); 273 ierr = PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL);CHKERRQ(ierr); 274 ierr = PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL);CHKERRQ(ierr); 275 ierr = PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL);CHKERRQ(ierr); 276 ierr = PetscOptionsTail();CHKERRQ(ierr); 277 PetscFunctionReturn(0); 278 } 279 280 static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer) 281 { 282 SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data; 283 PetscErrorCode ierr; 284 PetscBool iascii; 285 286 PetscFunctionBegin; 287 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 288 if (iascii) { 289 ierr = PetscViewerASCIIPrintf(viewer," Trust region tolerance (-snes_trtol)\n",(double)snes->deltatol);CHKERRQ(ierr); 290 ierr = PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma);CHKERRQ(ierr); 291 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); 292 } 293 PetscFunctionReturn(0); 294 } 295 /* ------------------------------------------------------------ */ 296 /*MC 297 SNESNEWTONTR - Newton based nonlinear solver that uses a trust region 298 299 Options Database: 300 + -snes_trtol <tol> - trust region tolerance 301 . -snes_tr_mu <mu> - trust region parameter 302 . -snes_tr_eta <eta> - trust region parameter 303 . -snes_tr_sigma <sigma> - trust region parameter 304 . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x) 305 . -snes_tr_delta1 <delta1> - trust region parameter 306 . -snes_tr_delta2 <delta2> - trust region parameter 307 - -snes_tr_delta3 <delta3> - trust region parameter 308 309 The basic algorithm is taken from "The Minpack Project", by More', 310 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 311 of Mathematical Software", Wayne Cowell, editor. 312 313 Level: intermediate 314 315 .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance() 316 317 M*/ 318 PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes) 319 { 320 SNES_NEWTONTR *neP; 321 PetscErrorCode ierr; 322 323 PetscFunctionBegin; 324 snes->ops->setup = SNESSetUp_NEWTONTR; 325 snes->ops->solve = SNESSolve_NEWTONTR; 326 snes->ops->destroy = SNESDestroy_NEWTONTR; 327 snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR; 328 snes->ops->view = SNESView_NEWTONTR; 329 snes->ops->reset = SNESReset_NEWTONTR; 330 331 snes->usesksp = PETSC_TRUE; 332 snes->usesnpc = PETSC_FALSE; 333 334 snes->alwayscomputesfinalresidual = PETSC_TRUE; 335 336 ierr = PetscNewLog(snes,&neP);CHKERRQ(ierr); 337 snes->data = (void*)neP; 338 neP->mu = 0.25; 339 neP->eta = 0.75; 340 neP->delta = 0.0; 341 neP->delta0 = 0.2; 342 neP->delta1 = 0.3; 343 neP->delta2 = 0.75; 344 neP->delta3 = 2.0; 345 neP->sigma = 0.0001; 346 neP->itflag = PETSC_FALSE; 347 neP->rnorm0 = 0.0; 348 neP->ttol = 0.0; 349 PetscFunctionReturn(0); 350 } 351 352