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