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_TR *neP = (SNES_TR*)snes->data; 20 Vec x; 21 PetscReal nrm; 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 ierr = KSPDefaultConverged(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,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",neP->delta,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 = KSPDefaultConvergedDestroy(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_TR *neP = (SNES_TR *)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",neP->delta,xnorm,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_TR - 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_TR" 86 static PetscErrorCode SNESSolve_TR(SNES snes) 87 { 88 SNES_TR *neP = (SNES_TR*)snes->data; 89 Vec X,F,Y,G,Ytmp; 90 PetscErrorCode ierr; 91 PetscInt maxits,i,lits; 92 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 93 PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1; 94 PetscScalar cnorm; 95 KSP ksp; 96 SNESConvergedReason reason = SNES_CONVERGED_ITERATING; 97 PetscBool conv = PETSC_FALSE,breakout = PETSC_FALSE; 98 99 PetscFunctionBegin; 100 maxits = snes->max_its; /* maximum number of iterations */ 101 X = snes->vec_sol; /* solution vector */ 102 F = snes->vec_func; /* residual vector */ 103 Y = snes->work[0]; /* work vectors */ 104 G = snes->work[1]; 105 Ytmp = snes->work[2]; 106 107 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 108 snes->iter = 0; 109 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 110 111 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); /* F(X) */ 112 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- || F || */ 113 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 114 snes->norm = fnorm; 115 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 116 delta = neP->delta0*fnorm; 117 neP->delta = delta; 118 SNESLogConvHistory(snes,fnorm,0); 119 SNESMonitor(snes,0,fnorm); 120 121 /* set parameter for default relative tolerance convergence test */ 122 snes->ttol = fnorm*snes->rtol; 123 /* test convergence */ 124 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 125 if (snes->reason) PetscFunctionReturn(0); 126 127 /* Set the stopping criteria to use the More' trick. */ 128 ierr = PetscOptionsGetBool(PETSC_NULL,"-snes_tr_ksp_regular_convergence_test",&conv,PETSC_NULL);CHKERRQ(ierr); 129 if (!conv) { 130 SNES_TR_KSPConverged_Ctx *ctx; 131 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 132 ierr = PetscNew(SNES_TR_KSPConverged_Ctx,&ctx);CHKERRQ(ierr); 133 ctx->snes = snes; 134 ierr = KSPDefaultConvergedCreate(&ctx->ctx);CHKERRQ(ierr); 135 ierr = KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,ctx,SNES_TR_KSPConverged_Destroy);CHKERRQ(ierr); 136 ierr = PetscInfo(snes,"Using Krylov convergence test SNES_TR_KSPConverged_Private\n");CHKERRQ(ierr); 137 } 138 139 for (i=0; i<maxits; i++) { 140 141 /* Call general purpose update function */ 142 if (snes->ops->update) { 143 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 144 } 145 146 /* Solve J Y = F, where J is Jacobian matrix */ 147 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 148 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 149 ierr = SNES_KSPSolve(snes,snes->ksp,F,Ytmp);CHKERRQ(ierr); 150 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 151 snes->linear_its += lits; 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",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 = VecAYPX(Y,-1.0,X);CHKERRQ(ierr); /* Y <- X - Y */ 174 ierr = VecCopy(X,snes->vec_sol_update);CHKERRQ(ierr); 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",fnorm,gnorm,ynorm);CHKERRQ(ierr); 185 ierr = PetscInfo3(snes,"gpred=%G, rho=%G, delta=%G\n",gpnorm,rho,delta);CHKERRQ(ierr); 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 = SNES_TR_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 SNESMonitor(snes,i+1,fnorm); 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 = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 208 snes->iter = i+1; 209 snes->norm = fnorm; 210 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 211 SNESLogConvHistory(snes,snes->norm,lits); 212 SNESMonitor(snes,snes->iter,snes->norm); 213 /* Test for convergence, xnorm = || X || */ 214 neP->itflag = PETSC_TRUE; 215 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 216 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 217 if (reason) break; 218 } else { 219 break; 220 } 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 = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 227 snes->reason = reason; 228 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 /*------------------------------------------------------------*/ 232 #undef __FUNCT__ 233 #define __FUNCT__ "SNESSetUp_TR" 234 static PetscErrorCode SNESSetUp_TR(SNES snes) 235 { 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 if (!snes->vec_sol_update) { 240 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 241 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 242 } 243 if (!snes->work) { 244 snes->nwork = 3; 245 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 246 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 247 } 248 PetscFunctionReturn(0); 249 } 250 /*------------------------------------------------------------*/ 251 #undef __FUNCT__ 252 #define __FUNCT__ "SNESDestroy_TR" 253 static PetscErrorCode SNESDestroy_TR(SNES snes) 254 { 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 if (snes->vec_sol_update) { 259 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 260 } 261 if (snes->work) { 262 ierr = VecDestroyVecs(snes->nwork,&snes->work);CHKERRQ(ierr); 263 } 264 ierr = PetscFree(snes->data);CHKERRQ(ierr); 265 PetscFunctionReturn(0); 266 } 267 /*------------------------------------------------------------*/ 268 269 #undef __FUNCT__ 270 #define __FUNCT__ "SNESSetFromOptions_TR" 271 static PetscErrorCode SNESSetFromOptions_TR(SNES snes) 272 { 273 SNES_TR *ctx = (SNES_TR *)snes->data; 274 PetscErrorCode ierr; 275 276 PetscFunctionBegin; 277 ierr = PetscOptionsHead("SNES trust region options for nonlinear equations");CHKERRQ(ierr); 278 ierr = PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,0);CHKERRQ(ierr); 279 ierr = PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,0);CHKERRQ(ierr); 280 ierr = PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,0);CHKERRQ(ierr); 281 ierr = PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,0);CHKERRQ(ierr); 282 ierr = PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,0);CHKERRQ(ierr); 283 ierr = PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,0);CHKERRQ(ierr); 284 ierr = PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,0);CHKERRQ(ierr); 285 ierr = PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,0);CHKERRQ(ierr); 286 ierr = PetscOptionsTail();CHKERRQ(ierr); 287 PetscFunctionReturn(0); 288 } 289 290 #undef __FUNCT__ 291 #define __FUNCT__ "SNESView_TR" 292 static PetscErrorCode SNESView_TR(SNES snes,PetscViewer viewer) 293 { 294 SNES_TR *tr = (SNES_TR *)snes->data; 295 PetscErrorCode ierr; 296 PetscBool iascii; 297 298 PetscFunctionBegin; 299 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 300 if (iascii) { 301 ierr = PetscViewerASCIIPrintf(viewer," mu=%G, eta=%G, sigma=%G\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 302 ierr = PetscViewerASCIIPrintf(viewer," delta0=%G, delta1=%G, delta2=%G, delta3=%G\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 303 } else { 304 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name); 305 } 306 PetscFunctionReturn(0); 307 } 308 /* ------------------------------------------------------------ */ 309 /*MC 310 SNESTR - Newton based nonlinear solver that uses a trust region 311 312 Options Database: 313 + -snes_trtol <tol> Trust region tolerance 314 . -snes_tr_mu <mu> 315 . -snes_tr_eta <eta> 316 . -snes_tr_sigma <sigma> 317 . -snes_tr_delta0 <delta0> 318 . -snes_tr_delta1 <delta1> 319 . -snes_tr_delta2 <delta2> 320 - -snes_tr_delta3 <delta3> 321 322 The basic algorithm is taken from "The Minpack Project", by More', 323 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 324 of Mathematical Software", Wayne Cowell, editor. 325 326 This is intended as a model implementation, since it does not 327 necessarily have many of the bells and whistles of other 328 implementations. 329 330 Level: intermediate 331 332 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESSetTrustRegionTolerance() 333 334 M*/ 335 EXTERN_C_BEGIN 336 #undef __FUNCT__ 337 #define __FUNCT__ "SNESCreate_TR" 338 PetscErrorCode SNESCreate_TR(SNES snes) 339 { 340 SNES_TR *neP; 341 PetscErrorCode ierr; 342 343 PetscFunctionBegin; 344 snes->ops->setup = SNESSetUp_TR; 345 snes->ops->solve = SNESSolve_TR; 346 snes->ops->destroy = SNESDestroy_TR; 347 snes->ops->setfromoptions = SNESSetFromOptions_TR; 348 snes->ops->view = SNESView_TR; 349 350 ierr = PetscNewLog(snes,SNES_TR,&neP);CHKERRQ(ierr); 351 snes->data = (void*)neP; 352 neP->mu = 0.25; 353 neP->eta = 0.75; 354 neP->delta = 0.0; 355 neP->delta0 = 0.2; 356 neP->delta1 = 0.3; 357 neP->delta2 = 0.75; 358 neP->delta3 = 2.0; 359 neP->sigma = 0.0001; 360 neP->itflag = PETSC_FALSE; 361 neP->rnorm0 = 0.0; 362 neP->ttol = 0.0; 363 PetscFunctionReturn(0); 364 } 365 EXTERN_C_END 366 367