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