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