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