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