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