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