1 2 #include "src/snes/impls/tr/tr.h" /*I "petscsnes.h" I*/ 3 4 /* 5 This convergence test determines if the two norm of the 6 solution lies outside the trust region, if so it halts. 7 */ 8 #undef __FUNCT__ 9 #define __FUNCT__ "SNES_TR_KSPConverged_Private" 10 int SNES_TR_KSPConverged_Private(KSP ksp,int n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx) 11 { 12 SNES snes = (SNES) ctx; 13 SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx; 14 SNES_TR *neP = (SNES_TR*)snes->data; 15 Vec x; 16 PetscReal nrm; 17 int ierr; 18 19 PetscFunctionBegin; 20 if (snes->ksp_ewconv) { 21 if (!kctx) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Eisenstat-Walker onvergence context not created"); 22 if (!n) {ierr = SNES_KSP_EW_ComputeRelativeTolerance_Private(snes,ksp);CHKERRQ(ierr);} 23 kctx->lresid_last = rnorm; 24 } 25 ierr = KSPDefaultConverged(ksp,n,rnorm,reason,ctx);CHKERRQ(ierr); 26 if (*reason) { 27 PetscLogInfo(snes,"SNES_TR_KSPConverged_Private: regular convergence test KSP iterations=%d, rnorm=%g\n",n,rnorm); 28 } 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 PetscLogInfo(snes,"SNES_TR_KSPConverged_Private: KSP iterations=%d, rnorm=%g\n",n,rnorm); 35 PetscLogInfo(snes,"SNES_TR_KSPConverged_Private: Ending linear iteration early, delta=%g, length=%g\n",neP->delta,nrm); 36 *reason = KSP_CONVERGED_STEP_LENGTH; 37 } 38 PetscFunctionReturn(0); 39 } 40 41 /* 42 SNESSolve_TR - Implements Newton's Method with a very simple trust 43 region approach for solving systems of nonlinear equations. 44 45 46 */ 47 #undef __FUNCT__ 48 #define __FUNCT__ "SNESSolve_TR" 49 static int SNESSolve_TR(SNES snes) 50 { 51 SNES_TR *neP = (SNES_TR*)snes->data; 52 Vec X,F,Y,G,TMP,Ytmp; 53 int maxits,i,ierr,lits,breakout = 0; 54 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 55 PetscReal rho,fnorm,gnorm,gpnorm,xnorm,delta,nrm,ynorm,norm1; 56 PetscScalar mone = -1.0,cnorm; 57 KSP ksp; 58 SNESConvergedReason reason; 59 PetscTruth conv; 60 61 PetscFunctionBegin; 62 maxits = snes->max_its; /* maximum number of iterations */ 63 X = snes->vec_sol; /* solution vector */ 64 F = snes->vec_func; /* residual vector */ 65 Y = snes->work[0]; /* work vectors */ 66 G = snes->work[1]; 67 Ytmp = snes->work[2]; 68 69 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 70 snes->iter = 0; 71 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 72 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm = || X || */ 73 74 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); /* F(X) */ 75 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- || F || */ 76 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 77 snes->norm = fnorm; 78 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 79 delta = neP->delta0*fnorm; 80 neP->delta = delta; 81 SNESLogConvHistory(snes,fnorm,0); 82 SNESMonitor(snes,0,fnorm); 83 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 84 85 if (fnorm < snes->atol) {snes->reason = SNES_CONVERGED_FNORM_ABS; PetscFunctionReturn(0);} 86 87 /* set parameter for default relative tolerance convergence test */ 88 snes->ttol = fnorm*snes->rtol; 89 90 /* Set the stopping criteria to use the More' trick. */ 91 ierr = PetscOptionsHasName(PETSC_NULL,"-snes_tr_ksp_regular_convergence_test",&conv);CHKERRQ(ierr); 92 if (!conv) { 93 ierr = KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,(void*)snes);CHKERRQ(ierr); 94 PetscLogInfo(snes,"SNESSolve_TR: Using Krylov convergence test SNES_TR_KSPConverged_Private\n"); 95 } 96 97 for (i=0; i<maxits; i++) { 98 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 99 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 100 101 /* Solve J Y = F, where J is Jacobian matrix */ 102 ierr = KSPSolve(snes->ksp,F,Ytmp);CHKERRQ(ierr); 103 ierr = KSPGetIterationNumber(ksp,&lits);CHKERRQ(ierr); 104 snes->linear_its += lits; 105 PetscLogInfo(snes,"SNESSolve_TR: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 106 ierr = VecNorm(Ytmp,NORM_2,&nrm);CHKERRQ(ierr); 107 norm1 = nrm; 108 while(1) { 109 ierr = VecCopy(Ytmp,Y);CHKERRQ(ierr); 110 nrm = norm1; 111 112 /* Scale Y if need be and predict new value of F norm */ 113 if (nrm >= delta) { 114 nrm = delta/nrm; 115 gpnorm = (1.0 - nrm)*fnorm; 116 cnorm = nrm; 117 PetscLogInfo(snes,"SNESSolve_TR: Scaling direction by %g\n",nrm); 118 ierr = VecScale(&cnorm,Y);CHKERRQ(ierr); 119 nrm = gpnorm; 120 ynorm = delta; 121 } else { 122 gpnorm = 0.0; 123 PetscLogInfo(snes,"SNESSolve_TR: Direction is in Trust Region\n"); 124 ynorm = nrm; 125 } 126 ierr = VecAYPX(&mone,X,Y);CHKERRQ(ierr); /* Y <- X - Y */ 127 ierr = VecCopy(X,snes->vec_sol_update_always);CHKERRQ(ierr); 128 ierr = SNESComputeFunction(snes,Y,G);CHKERRQ(ierr); /* F(X) */ 129 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); /* gnorm <- || g || */ 130 if (fnorm == gpnorm) rho = 0.0; 131 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 132 133 /* Update size of trust region */ 134 if (rho < neP->mu) delta *= neP->delta1; 135 else if (rho < neP->eta) delta *= neP->delta2; 136 else delta *= neP->delta3; 137 PetscLogInfo(snes,"SNESSolve_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 138 PetscLogInfo(snes,"SNESSolve_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta); 139 neP->delta = delta; 140 if (rho > neP->sigma) break; 141 PetscLogInfo(snes,"SNESSolve_TR: Trying again in smaller region\n"); 142 /* check to see if progress is hopeless */ 143 neP->itflag = 0; 144 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 145 if (reason) { 146 /* We're not progressing, so return with the current iterate */ 147 SNESMonitor(snes,i+1,fnorm); 148 breakout = 1; 149 break; 150 } 151 snes->numFailures++; 152 } 153 if (!breakout) { 154 fnorm = gnorm; 155 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 156 snes->iter = i+1; 157 snes->norm = fnorm; 158 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 159 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 160 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 161 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm = || X || */ 162 SNESLogConvHistory(snes,fnorm,lits); 163 SNESMonitor(snes,i+1,fnorm); 164 165 /* Test for convergence */ 166 neP->itflag = 1; 167 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 168 if (reason) { 169 break; 170 } 171 } else { 172 break; 173 } 174 } 175 /* Verify solution is in corect location */ 176 if (X != snes->vec_sol) { 177 ierr = VecCopy(X,snes->vec_sol);CHKERRQ(ierr); 178 } 179 if (F != snes->vec_func) { 180 ierr = VecCopy(F,snes->vec_func);CHKERRQ(ierr); 181 } 182 snes->vec_sol_always = snes->vec_sol; 183 snes->vec_func_always = snes->vec_func; 184 if (i == maxits) { 185 PetscLogInfo(snes,"SNESSolve_TR: Maximum number of iterations has been reached: %d\n",maxits); 186 reason = SNES_DIVERGED_MAX_IT; 187 } 188 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 189 snes->reason = reason; 190 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 191 PetscFunctionReturn(0); 192 } 193 /*------------------------------------------------------------*/ 194 #undef __FUNCT__ 195 #define __FUNCT__ "SNESSetUp_TR" 196 static int SNESSetUp_TR(SNES snes) 197 { 198 int ierr; 199 200 PetscFunctionBegin; 201 snes->nwork = 4; 202 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 203 PetscLogObjectParents(snes,snes->nwork,snes->work); 204 snes->vec_sol_update_always = snes->work[3]; 205 PetscFunctionReturn(0); 206 } 207 /*------------------------------------------------------------*/ 208 #undef __FUNCT__ 209 #define __FUNCT__ "SNESDestroy_TR" 210 static int SNESDestroy_TR(SNES snes) 211 { 212 int ierr; 213 214 PetscFunctionBegin; 215 if (snes->nwork) { 216 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 217 } 218 ierr = PetscFree(snes->data);CHKERRQ(ierr); 219 PetscFunctionReturn(0); 220 } 221 /*------------------------------------------------------------*/ 222 223 #undef __FUNCT__ 224 #define __FUNCT__ "SNESSetFromOptions_TR" 225 static int SNESSetFromOptions_TR(SNES snes) 226 { 227 SNES_TR *ctx = (SNES_TR *)snes->data; 228 int ierr; 229 230 PetscFunctionBegin; 231 ierr = PetscOptionsHead("SNES trust region options for nonlinear equations");CHKERRQ(ierr); 232 ierr = PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,0);CHKERRQ(ierr); 233 ierr = PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,0);CHKERRQ(ierr); 234 ierr = PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,0);CHKERRQ(ierr); 235 ierr = PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,0);CHKERRQ(ierr); 236 ierr = PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,0);CHKERRQ(ierr); 237 ierr = PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,0);CHKERRQ(ierr); 238 ierr = PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,0);CHKERRQ(ierr); 239 ierr = PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,0);CHKERRQ(ierr); 240 ierr = PetscOptionsTail();CHKERRQ(ierr); 241 PetscFunctionReturn(0); 242 } 243 244 #undef __FUNCT__ 245 #define __FUNCT__ "SNESView_TR" 246 static int SNESView_TR(SNES snes,PetscViewer viewer) 247 { 248 SNES_TR *tr = (SNES_TR *)snes->data; 249 int ierr; 250 PetscTruth iascii; 251 252 PetscFunctionBegin; 253 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 254 if (iascii) { 255 ierr = PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 256 ierr = PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 257 } else { 258 SETERRQ1(1,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name); 259 } 260 PetscFunctionReturn(0); 261 } 262 263 /* ---------------------------------------------------------------- */ 264 #undef __FUNCT__ 265 #define __FUNCT__ "SNESConverged_TR" 266 /*@C 267 SNESConverged_TR - Monitors the convergence of the trust region 268 method SNESTR for solving systems of nonlinear equations (default). 269 270 Collective on SNES 271 272 Input Parameters: 273 + snes - the SNES context 274 . xnorm - 2-norm of current iterate 275 . pnorm - 2-norm of current step 276 . fnorm - 2-norm of function 277 - dummy - unused context 278 279 Output Parameter: 280 . reason - one of 281 $ SNES_CONVERGED_FNORM_ABS - (fnorm < atol), 282 $ SNES_CONVERGED_PNORM_RELATIVE - (pnorm < xtol*xnorm), 283 $ SNES_CONVERGED_FNORM_RELATIVE - (fnorm < rtol*fnorm0), 284 $ SNES_DIVERGED_FUNCTION_COUNT - (nfct > maxf), 285 $ SNES_DIVERGED_FNORM_NAN - (fnorm == NaN), 286 $ SNES_CONVERGED_TR_DELTA - (delta < xnorm*deltatol), 287 $ SNES_CONVERGED_ITERATING - (otherwise) 288 289 where 290 + delta - trust region paramenter 291 . deltatol - trust region size tolerance, 292 set with SNESSetTrustRegionTolerance() 293 . maxf - maximum number of function evaluations, 294 set with SNESSetTolerances() 295 . nfct - number of function evaluations, 296 . atol - absolute function norm tolerance, 297 set with SNESSetTolerances() 298 - xtol - relative function norm tolerance, 299 set with SNESSetTolerances() 300 301 Level: intermediate 302 303 .keywords: SNES, nonlinear, default, converged, convergence 304 305 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 306 @*/ 307 int SNESConverged_TR(SNES snes,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) 308 { 309 SNES_TR *neP = (SNES_TR *)snes->data; 310 int ierr; 311 312 PetscFunctionBegin; 313 if (fnorm != fnorm) { 314 PetscLogInfo(snes,"SNESConverged_TR:Failed to converged, function norm is NaN\n"); 315 *reason = SNES_DIVERGED_FNORM_NAN; 316 } else if (neP->delta < xnorm * snes->deltatol) { 317 PetscLogInfo(snes,"SNESConverged_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 318 *reason = SNES_CONVERGED_TR_DELTA; 319 } else if (neP->itflag) { 320 ierr = SNESConverged_LS(snes,xnorm,pnorm,fnorm,reason,dummy);CHKERRQ(ierr); 321 } else if (snes->nfuncs > snes->max_funcs) { 322 PetscLogInfo(snes,"SNESConverged_TR: Exceeded maximum number of function evaluations: %d > %d\n",snes->nfuncs,snes->max_funcs); 323 *reason = SNES_DIVERGED_FUNCTION_COUNT; 324 } else { 325 *reason = SNES_CONVERGED_ITERATING; 326 } 327 PetscFunctionReturn(0); 328 } 329 /* ------------------------------------------------------------ */ 330 /*MC 331 SNESTR - Newton based nonlinear solver that uses a trust region 332 333 Options Database: 334 + -snes_trtol <tol> Trust region tolerance 335 . -snes_tr_mu <mu> 336 . -snes_tr_eta <eta> 337 . -snes_tr_sigma <sigma> 338 . -snes_tr_delta0 <delta0> 339 . -snes_tr_delta1 <delta1> 340 . -snes_tr_delta2 <delta2> 341 - -snes_tr_delta3 <delta3> 342 343 The basic algorithm is taken from "The Minpack Project", by More', 344 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 345 of Mathematical Software", Wayne Cowell, editor. 346 347 This is intended as a model implementation, since it does not 348 necessarily have many of the bells and whistles of other 349 implementations. 350 351 Level: intermediate 352 353 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESSetTrustRegionTolerance() 354 355 M*/ 356 EXTERN_C_BEGIN 357 #undef __FUNCT__ 358 #define __FUNCT__ "SNESCreate_TR" 359 int SNESCreate_TR(SNES snes) 360 { 361 SNES_TR *neP; 362 int ierr; 363 364 PetscFunctionBegin; 365 snes->setup = SNESSetUp_TR; 366 snes->solve = SNESSolve_TR; 367 snes->destroy = SNESDestroy_TR; 368 snes->converged = SNESConverged_TR; 369 snes->setfromoptions = SNESSetFromOptions_TR; 370 snes->view = SNESView_TR; 371 snes->nwork = 0; 372 373 ierr = PetscNew(SNES_TR,&neP);CHKERRQ(ierr); 374 PetscLogObjectMemory(snes,sizeof(SNES_TR)); 375 snes->data = (void*)neP; 376 neP->mu = 0.25; 377 neP->eta = 0.75; 378 neP->delta = 0.0; 379 neP->delta0 = 0.2; 380 neP->delta1 = 0.3; 381 neP->delta2 = 0.75; 382 neP->delta3 = 2.0; 383 neP->sigma = 0.0001; 384 neP->itflag = 0; 385 neP->rnorm0 = 0; 386 neP->ttol = 0; 387 PetscFunctionReturn(0); 388 } 389 EXTERN_C_END 390 391