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