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_EQ_TR_KSPConverged_Private" 11 int SNES_EQ_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_EQ_TR *neP = (SNES_EQ_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_EQ_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_EQ_TR_KSPConverged_Private: KSP iterations=%d, rnorm=%g\n",n,rnorm); 36 PetscLogInfo(snes,"SNES_EQ_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_EQ_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_EQ_TR" 56 static int SNESSolve_EQ_TR(SNES snes,int *its) 57 { 58 SNES_EQ_TR *neP = (SNES_EQ_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_EQ_TR_KSPConverged_Private,(void *)snes);CHKERRQ(ierr); 103 PetscLogInfo(snes,"SNESSolve_EQ_TR: Using Krylov convergence test SNES_EQ_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_EQ_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_EQ_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_EQ_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_EQ_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 146 PetscLogInfo(snes,"SNESSolve_EQ_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_EQ_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_EQ_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_EQ_TR" 206 static int SNESSetUp_EQ_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_EQ_TR" 220 static int SNESDestroy_EQ_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_EQ_TR" 235 static int SNESSetFromOptions_EQ_TR(SNES snes) 236 { 237 SNES_EQ_TR *ctx = (SNES_EQ_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_eq_tr_mu","mu","None",ctx->mu,&ctx->mu,0);CHKERRQ(ierr); 244 ierr = PetscOptionsReal("-snes_eq_tr_eta","eta","None",ctx->eta,&ctx->eta,0);CHKERRQ(ierr); 245 ierr = PetscOptionsReal("-snes_eq_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,0);CHKERRQ(ierr); 246 ierr = PetscOptionsReal("-snes_eq_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,0);CHKERRQ(ierr); 247 ierr = PetscOptionsReal("-snes_eq_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,0);CHKERRQ(ierr); 248 ierr = PetscOptionsReal("-snes_eq_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,0);CHKERRQ(ierr); 249 ierr = PetscOptionsReal("-snes_eq_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_EQ_TR" 256 static int SNESView_EQ_TR(SNES snes,PetscViewer viewer) 257 { 258 SNES_EQ_TR *tr = (SNES_EQ_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_EQ_TR" 276 /*@C 277 SNESConverged_EQ_TR - Monitors the convergence of the trust region 278 method SNESEQTR 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_EQ_TR(SNES snes,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) 318 { 319 SNES_EQ_TR *neP = (SNES_EQ_TR *)snes->data; 320 int ierr; 321 322 PetscFunctionBegin; 323 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 324 SETERRQ(PETSC_ERR_ARG_WRONG,"For SNES_NONLINEAR_EQUATIONS only"); 325 } 326 327 if (fnorm != fnorm) { 328 PetscLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 329 *reason = SNES_DIVERGED_FNORM_NAN; 330 } else if (neP->delta < xnorm * snes->deltatol) { 331 PetscLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 332 *reason = SNES_CONVERGED_TR_DELTA; 333 } else if (neP->itflag) { 334 ierr = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,reason,dummy);CHKERRQ(ierr); 335 } else if (snes->nfuncs > snes->max_funcs) { 336 PetscLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n",snes->nfuncs,snes->max_funcs); 337 *reason = SNES_DIVERGED_FUNCTION_COUNT; 338 } else { 339 *reason = SNES_CONVERGED_ITERATING; 340 } 341 PetscFunctionReturn(0); 342 } 343 /* ------------------------------------------------------------ */ 344 EXTERN_C_BEGIN 345 #undef __FUNCT__ 346 #define __FUNCT__ "SNESCreate_EQ_TR" 347 int SNESCreate_EQ_TR(SNES snes) 348 { 349 SNES_EQ_TR *neP; 350 int ierr; 351 352 PetscFunctionBegin; 353 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 354 SETERRQ(PETSC_ERR_ARG_WRONG,"For SNES_NONLINEAR_EQUATIONS only"); 355 } 356 snes->setup = SNESSetUp_EQ_TR; 357 snes->solve = SNESSolve_EQ_TR; 358 snes->destroy = SNESDestroy_EQ_TR; 359 snes->converged = SNESConverged_EQ_TR; 360 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 361 snes->view = SNESView_EQ_TR; 362 snes->nwork = 0; 363 364 ierr = PetscNew(SNES_EQ_TR,&neP);CHKERRQ(ierr); 365 PetscLogObjectMemory(snes,sizeof(SNES_EQ_TR)); 366 snes->data = (void*)neP; 367 neP->mu = 0.25; 368 neP->eta = 0.75; 369 neP->delta = 0.0; 370 neP->delta0 = 0.2; 371 neP->delta1 = 0.3; 372 neP->delta2 = 0.75; 373 neP->delta3 = 2.0; 374 neP->sigma = 0.0001; 375 neP->itflag = 0; 376 neP->rnorm0 = 0; 377 neP->ttol = 0; 378 PetscFunctionReturn(0); 379 } 380 EXTERN_C_END 381 382