1 #define PETSCSNES_DLL 2 3 #include "src/snes/impls/tr/tr.h" /*I "petscsnes.h" I*/ 4 5 typedef struct { 6 void *ctx; 7 SNES snes; 8 } SNES_TR_KSPConverged_Ctx; 9 10 /* 11 This convergence test determines if the two norm of the 12 solution lies outside the trust region, if so it halts. 13 */ 14 #undef __FUNCT__ 15 #define __FUNCT__ "SNES_TR_KSPConverged_Private" 16 PetscErrorCode SNES_TR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx) 17 { 18 SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx; 19 SNES snes = ctx->snes; 20 SNES_TR *neP = (SNES_TR*)snes->data; 21 Vec x; 22 PetscReal nrm; 23 PetscErrorCode ierr; 24 25 PetscFunctionBegin; 26 ierr = KSPDefaultConverged(ksp,n,rnorm,reason,ctx->ctx);CHKERRQ(ierr); 27 if (*reason) { 28 ierr = PetscInfo2(snes,"default convergence test KSP iterations=%D, rnorm=%G\n",n,rnorm);CHKERRQ(ierr); 29 } 30 /* Determine norm of solution */ 31 ierr = KSPBuildSolution(ksp,0,&x);CHKERRQ(ierr); 32 ierr = VecNorm(x,NORM_2,&nrm);CHKERRQ(ierr); 33 if (nrm >= neP->delta) { 34 ierr = PetscInfo2(snes,"Ending linear iteration early, delta=%G, length=%G\n",neP->delta,nrm);CHKERRQ(ierr); 35 *reason = KSP_CONVERGED_STEP_LENGTH; 36 } 37 PetscFunctionReturn(0); 38 } 39 40 #undef __FUNCT__ 41 #define __FUNCT__ "SNES_TR_KSPConverged_Destroy" 42 PetscErrorCode SNES_TR_KSPConverged_Destroy(void *cctx) 43 { 44 SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx *)cctx; 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 ierr = KSPDefaultConvergedDestroy(ctx->ctx);CHKERRQ(ierr); 49 ierr = PetscFree(ctx);CHKERRQ(ierr); 50 PetscFunctionReturn(0); 51 } 52 53 /* ---------------------------------------------------------------- */ 54 #undef __FUNCT__ 55 #define __FUNCT__ "SNES_TR_Converged_Private" 56 /* 57 SNES_TR_Converged_Private -test convergence JUST for 58 the trust region tolerance. 59 60 */ 61 static PetscErrorCode SNES_TR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) 62 { 63 SNES_TR *neP = (SNES_TR *)snes->data; 64 PetscErrorCode ierr; 65 66 PetscFunctionBegin; 67 *reason = SNES_CONVERGED_ITERATING; 68 if (neP->delta < xnorm * snes->deltatol) { 69 ierr = PetscInfo3(snes,"Converged due to trust region param %G<%G*%G\n",neP->delta,xnorm,snes->deltatol);CHKERRQ(ierr); 70 *reason = SNES_CONVERGED_TR_DELTA; 71 } else if (snes->nfuncs >= snes->max_funcs) { 72 ierr = PetscInfo1(snes,"Exceeded maximum number of function evaluations: %D\n",snes->max_funcs);CHKERRQ(ierr); 73 *reason = SNES_DIVERGED_FUNCTION_COUNT; 74 } 75 PetscFunctionReturn(0); 76 } 77 78 79 /* 80 SNESSolve_TR - Implements Newton's Method with a very simple trust 81 region approach for solving systems of nonlinear equations. 82 83 84 */ 85 #undef __FUNCT__ 86 #define __FUNCT__ "SNESSolve_TR" 87 static PetscErrorCode SNESSolve_TR(SNES snes) 88 { 89 SNES_TR *neP = (SNES_TR*)snes->data; 90 Vec X,F,Y,G,Ytmp; 91 PetscErrorCode ierr; 92 PetscInt maxits,i,lits; 93 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 94 PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1; 95 PetscScalar cnorm; 96 KSP ksp; 97 SNESConvergedReason reason = SNES_CONVERGED_ITERATING; 98 PetscTruth conv,breakout = PETSC_FALSE; 99 100 PetscFunctionBegin; 101 maxits = snes->max_its; /* maximum number of iterations */ 102 X = snes->vec_sol; /* solution vector */ 103 F = snes->vec_func; /* residual vector */ 104 Y = snes->work[0]; /* work vectors */ 105 G = snes->work[1]; 106 Ytmp = snes->work[2]; 107 108 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 109 snes->iter = 0; 110 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 111 112 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); /* F(X) */ 113 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- || F || */ 114 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 115 snes->norm = fnorm; 116 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 117 delta = neP->delta0*fnorm; 118 neP->delta = delta; 119 SNESLogConvHistory(snes,fnorm,0); 120 SNESMonitor(snes,0,fnorm); 121 122 /* set parameter for default relative tolerance convergence test */ 123 snes->ttol = fnorm*snes->rtol; 124 /* test convergence */ 125 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 126 if (snes->reason) PetscFunctionReturn(0); 127 128 /* Set the stopping criteria to use the More' trick. */ 129 ierr = PetscOptionsHasName(PETSC_NULL,"-snes_tr_ksp_regular_convergence_test",&conv);CHKERRQ(ierr); 130 if (!conv) { 131 SNES_TR_KSPConverged_Ctx *ctx; 132 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 133 ierr = PetscNew(SNES_TR_KSPConverged_Ctx,&ctx);CHKERRQ(ierr); 134 ctx->snes = snes; 135 ierr = KSPDefaultConvergedCreate(&ctx->ctx);CHKERRQ(ierr); 136 ierr = KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,ctx,SNES_TR_KSPConverged_Destroy);CHKERRQ(ierr); 137 ierr = PetscInfo(snes,"Using Krylov convergence test SNES_TR_KSPConverged_Private\n");CHKERRQ(ierr); 138 } 139 140 for (i=0; i<maxits; i++) { 141 142 /* Call general purpose update function */ 143 if (snes->ops->update) { 144 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 145 } 146 147 /* Solve J Y = F, where J is Jacobian matrix */ 148 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 149 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 150 ierr = SNES_KSPSolve(snes,snes->ksp,F,Ytmp);CHKERRQ(ierr); 151 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 152 snes->linear_its += lits; 153 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 154 ierr = VecNorm(Ytmp,NORM_2,&nrm);CHKERRQ(ierr); 155 norm1 = nrm; 156 while(1) { 157 ierr = VecCopy(Ytmp,Y);CHKERRQ(ierr); 158 nrm = norm1; 159 160 /* Scale Y if need be and predict new value of F norm */ 161 if (nrm >= delta) { 162 nrm = delta/nrm; 163 gpnorm = (1.0 - nrm)*fnorm; 164 cnorm = nrm; 165 ierr = PetscInfo1(snes,"Scaling direction by %G\n",nrm);CHKERRQ(ierr); 166 ierr = VecScale(Y,cnorm);CHKERRQ(ierr); 167 nrm = gpnorm; 168 ynorm = delta; 169 } else { 170 gpnorm = 0.0; 171 ierr = PetscInfo(snes,"Direction is in Trust Region\n");CHKERRQ(ierr); 172 ynorm = nrm; 173 } 174 ierr = VecAYPX(Y,-1.0,X);CHKERRQ(ierr); /* Y <- X - Y */ 175 ierr = VecCopy(X,snes->vec_sol_update);CHKERRQ(ierr); 176 ierr = SNESComputeFunction(snes,Y,G);CHKERRQ(ierr); /* F(X) */ 177 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); /* gnorm <- || g || */ 178 if (fnorm == gpnorm) rho = 0.0; 179 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 180 181 /* Update size of trust region */ 182 if (rho < neP->mu) delta *= neP->delta1; 183 else if (rho < neP->eta) delta *= neP->delta2; 184 else delta *= neP->delta3; 185 ierr = PetscInfo3(snes,"fnorm=%G, gnorm=%G, ynorm=%G\n",fnorm,gnorm,ynorm);CHKERRQ(ierr); 186 ierr = PetscInfo3(snes,"gpred=%G, rho=%G, delta=%G\n",gpnorm,rho,delta);CHKERRQ(ierr); 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 = SNES_TR_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 SNESMonitor(snes,i+1,fnorm); 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 = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 209 snes->iter = i+1; 210 snes->norm = fnorm; 211 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 212 SNESLogConvHistory(snes,snes->norm,lits); 213 SNESMonitor(snes,snes->iter,snes->norm); 214 /* Test for convergence, xnorm = || X || */ 215 neP->itflag = PETSC_TRUE; 216 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 217 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 218 if (reason) break; 219 } else { 220 break; 221 } 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 = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 228 snes->reason = reason; 229 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 230 PetscFunctionReturn(0); 231 } 232 /*------------------------------------------------------------*/ 233 #undef __FUNCT__ 234 #define __FUNCT__ "SNESSetUp_TR" 235 static PetscErrorCode SNESSetUp_TR(SNES snes) 236 { 237 PetscErrorCode ierr; 238 239 PetscFunctionBegin; 240 if (!snes->vec_sol_update) { 241 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 242 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 243 } 244 if (!snes->work) { 245 snes->nwork = 3; 246 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 247 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 248 } 249 PetscFunctionReturn(0); 250 } 251 /*------------------------------------------------------------*/ 252 #undef __FUNCT__ 253 #define __FUNCT__ "SNESDestroy_TR" 254 static PetscErrorCode SNESDestroy_TR(SNES snes) 255 { 256 PetscErrorCode ierr; 257 258 PetscFunctionBegin; 259 if (snes->vec_sol_update) { 260 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 261 snes->vec_sol_update = PETSC_NULL; 262 } 263 if (snes->nwork) { 264 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 265 snes->nwork = 0; 266 snes->work = PETSC_NULL; 267 } 268 ierr = PetscFree(snes->data);CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 /*------------------------------------------------------------*/ 272 273 #undef __FUNCT__ 274 #define __FUNCT__ "SNESSetFromOptions_TR" 275 static PetscErrorCode SNESSetFromOptions_TR(SNES snes) 276 { 277 SNES_TR *ctx = (SNES_TR *)snes->data; 278 PetscErrorCode ierr; 279 280 PetscFunctionBegin; 281 ierr = PetscOptionsHead("SNES trust region options for nonlinear equations");CHKERRQ(ierr); 282 ierr = PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,0);CHKERRQ(ierr); 283 ierr = PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,0);CHKERRQ(ierr); 284 ierr = PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,0);CHKERRQ(ierr); 285 ierr = PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,0);CHKERRQ(ierr); 286 ierr = PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,0);CHKERRQ(ierr); 287 ierr = PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,0);CHKERRQ(ierr); 288 ierr = PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,0);CHKERRQ(ierr); 289 ierr = PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,0);CHKERRQ(ierr); 290 ierr = PetscOptionsTail();CHKERRQ(ierr); 291 PetscFunctionReturn(0); 292 } 293 294 #undef __FUNCT__ 295 #define __FUNCT__ "SNESView_TR" 296 static PetscErrorCode SNESView_TR(SNES snes,PetscViewer viewer) 297 { 298 SNES_TR *tr = (SNES_TR *)snes->data; 299 PetscErrorCode ierr; 300 PetscTruth iascii; 301 302 PetscFunctionBegin; 303 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 304 if (iascii) { 305 ierr = PetscViewerASCIIPrintf(viewer," mu=%G, eta=%G, sigma=%G\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 306 ierr = PetscViewerASCIIPrintf(viewer," delta0=%G, delta1=%G, delta2=%G, delta3=%G\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 307 } else { 308 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name); 309 } 310 PetscFunctionReturn(0); 311 } 312 /* ------------------------------------------------------------ */ 313 /*MC 314 SNESTR - Newton based nonlinear solver that uses a trust region 315 316 Options Database: 317 + -snes_trtol <tol> Trust region tolerance 318 . -snes_tr_mu <mu> 319 . -snes_tr_eta <eta> 320 . -snes_tr_sigma <sigma> 321 . -snes_tr_delta0 <delta0> 322 . -snes_tr_delta1 <delta1> 323 . -snes_tr_delta2 <delta2> 324 - -snes_tr_delta3 <delta3> 325 326 The basic algorithm is taken from "The Minpack Project", by More', 327 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 328 of Mathematical Software", Wayne Cowell, editor. 329 330 This is intended as a model implementation, since it does not 331 necessarily have many of the bells and whistles of other 332 implementations. 333 334 Level: intermediate 335 336 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESSetTrustRegionTolerance() 337 338 M*/ 339 EXTERN_C_BEGIN 340 #undef __FUNCT__ 341 #define __FUNCT__ "SNESCreate_TR" 342 PetscErrorCode PETSCSNES_DLLEXPORT SNESCreate_TR(SNES snes) 343 { 344 SNES_TR *neP; 345 PetscErrorCode ierr; 346 347 PetscFunctionBegin; 348 snes->ops->setup = SNESSetUp_TR; 349 snes->ops->solve = SNESSolve_TR; 350 snes->ops->destroy = SNESDestroy_TR; 351 snes->ops->setfromoptions = SNESSetFromOptions_TR; 352 snes->ops->view = SNESView_TR; 353 354 ierr = PetscNewLog(snes,SNES_TR,&neP);CHKERRQ(ierr); 355 snes->data = (void*)neP; 356 neP->mu = 0.25; 357 neP->eta = 0.75; 358 neP->delta = 0.0; 359 neP->delta0 = 0.2; 360 neP->delta1 = 0.3; 361 neP->delta2 = 0.75; 362 neP->delta3 = 2.0; 363 neP->sigma = 0.0001; 364 neP->itflag = PETSC_FALSE; 365 neP->rnorm0 = 0; 366 neP->ttol = 0; 367 PetscFunctionReturn(0); 368 } 369 EXTERN_C_END 370 371