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 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) {*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 = 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 i--; 196 reason = SNES_DIVERGED_MAX_IT; 197 } 198 *its = i+1; 199 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 200 snes->reason = reason; 201 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 202 PetscFunctionReturn(0); 203 } 204 /*------------------------------------------------------------*/ 205 #undef __FUNCT__ 206 #define __FUNCT__ "SNESSetUp_TR" 207 static int SNESSetUp_TR(SNES snes) 208 { 209 int ierr; 210 211 PetscFunctionBegin; 212 snes->nwork = 4; 213 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 214 PetscLogObjectParents(snes,snes->nwork,snes->work); 215 snes->vec_sol_update_always = snes->work[3]; 216 PetscFunctionReturn(0); 217 } 218 /*------------------------------------------------------------*/ 219 #undef __FUNCT__ 220 #define __FUNCT__ "SNESDestroy_TR" 221 static int SNESDestroy_TR(SNES snes) 222 { 223 int ierr; 224 225 PetscFunctionBegin; 226 if (snes->nwork) { 227 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 228 } 229 ierr = PetscFree(snes->data);CHKERRQ(ierr); 230 PetscFunctionReturn(0); 231 } 232 /*------------------------------------------------------------*/ 233 234 #undef __FUNCT__ 235 #define __FUNCT__ "SNESSetFromOptions_TR" 236 static int SNESSetFromOptions_TR(SNES snes) 237 { 238 SNES_TR *ctx = (SNES_TR *)snes->data; 239 int ierr; 240 241 PetscFunctionBegin; 242 ierr = PetscOptionsHead("SNES trust region options for nonlinear equations");CHKERRQ(ierr); 243 ierr = PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,0);CHKERRQ(ierr); 244 ierr = PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,0);CHKERRQ(ierr); 245 ierr = PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,0);CHKERRQ(ierr); 246 ierr = PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,0);CHKERRQ(ierr); 247 ierr = PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,0);CHKERRQ(ierr); 248 ierr = PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,0);CHKERRQ(ierr); 249 ierr = PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,0);CHKERRQ(ierr); 250 ierr = PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,0);CHKERRQ(ierr); 251 ierr = PetscOptionsTail();CHKERRQ(ierr); 252 PetscFunctionReturn(0); 253 } 254 255 #undef __FUNCT__ 256 #define __FUNCT__ "SNESView_TR" 257 static int SNESView_TR(SNES snes,PetscViewer viewer) 258 { 259 SNES_TR *tr = (SNES_TR *)snes->data; 260 int ierr; 261 PetscTruth isascii; 262 263 PetscFunctionBegin; 264 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 265 if (isascii) { 266 ierr = PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 267 ierr = PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 268 } else { 269 SETERRQ1(1,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name); 270 } 271 PetscFunctionReturn(0); 272 } 273 274 /* ---------------------------------------------------------------- */ 275 #undef __FUNCT__ 276 #define __FUNCT__ "SNESConverged_TR" 277 /*@C 278 SNESConverged_TR - Monitors the convergence of the trust region 279 method SNESTR for solving systems of nonlinear equations (default). 280 281 Collective on SNES 282 283 Input Parameters: 284 + snes - the SNES context 285 . xnorm - 2-norm of current iterate 286 . pnorm - 2-norm of current step 287 . fnorm - 2-norm of function 288 - dummy - unused context 289 290 Output Parameter: 291 . reason - one of 292 $ SNES_CONVERGED_FNORM_ABS - (fnorm < atol), 293 $ SNES_CONVERGED_PNORM_RELATIVE - (pnorm < xtol*xnorm), 294 $ SNES_CONVERGED_FNORM_RELATIVE - (fnorm < rtol*fnorm0), 295 $ SNES_DIVERGED_FUNCTION_COUNT - (nfct > maxf), 296 $ SNES_DIVERGED_FNORM_NAN - (fnorm == NaN), 297 $ SNES_CONVERGED_TR_DELTA - (delta < xnorm*deltatol), 298 $ SNES_CONVERGED_ITERATING - (otherwise) 299 300 where 301 + delta - trust region paramenter 302 . deltatol - trust region size tolerance, 303 set with SNESSetTrustRegionTolerance() 304 . maxf - maximum number of function evaluations, 305 set with SNESSetTolerances() 306 . nfct - number of function evaluations, 307 . atol - absolute function norm tolerance, 308 set with SNESSetTolerances() 309 - xtol - relative function norm tolerance, 310 set with SNESSetTolerances() 311 312 Level: intermediate 313 314 .keywords: SNES, nonlinear, default, converged, convergence 315 316 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 317 @*/ 318 int SNESConverged_TR(SNES snes,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) 319 { 320 SNES_TR *neP = (SNES_TR *)snes->data; 321 int ierr; 322 323 PetscFunctionBegin; 324 if (fnorm != fnorm) { 325 PetscLogInfo(snes,"SNESConverged_TR:Failed to converged, function norm is NaN\n"); 326 *reason = SNES_DIVERGED_FNORM_NAN; 327 } else if (neP->delta < xnorm * snes->deltatol) { 328 PetscLogInfo(snes,"SNESConverged_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 329 *reason = SNES_CONVERGED_TR_DELTA; 330 } else if (neP->itflag) { 331 ierr = SNESConverged_LS(snes,xnorm,pnorm,fnorm,reason,dummy);CHKERRQ(ierr); 332 } else if (snes->nfuncs > snes->max_funcs) { 333 PetscLogInfo(snes,"SNESConverged_TR: Exceeded maximum number of function evaluations: %d > %d\n",snes->nfuncs,snes->max_funcs); 334 *reason = SNES_DIVERGED_FUNCTION_COUNT; 335 } else { 336 *reason = SNES_CONVERGED_ITERATING; 337 } 338 PetscFunctionReturn(0); 339 } 340 /* ------------------------------------------------------------ */ 341 EXTERN_C_BEGIN 342 #undef __FUNCT__ 343 #define __FUNCT__ "SNESCreate_TR" 344 int SNESCreate_TR(SNES snes) 345 { 346 SNES_TR *neP; 347 int ierr; 348 349 PetscFunctionBegin; 350 snes->setup = SNESSetUp_TR; 351 snes->solve = SNESSolve_TR; 352 snes->destroy = SNESDestroy_TR; 353 snes->converged = SNESConverged_TR; 354 snes->setfromoptions = SNESSetFromOptions_TR; 355 snes->view = SNESView_TR; 356 snes->nwork = 0; 357 358 ierr = PetscNew(SNES_TR,&neP);CHKERRQ(ierr); 359 PetscLogObjectMemory(snes,sizeof(SNES_TR)); 360 snes->data = (void*)neP; 361 neP->mu = 0.25; 362 neP->eta = 0.75; 363 neP->delta = 0.0; 364 neP->delta0 = 0.2; 365 neP->delta1 = 0.3; 366 neP->delta2 = 0.75; 367 neP->delta3 = 2.0; 368 neP->sigma = 0.0001; 369 neP->itflag = 0; 370 neP->rnorm0 = 0; 371 neP->ttol = 0; 372 PetscFunctionReturn(0); 373 } 374 EXTERN_C_END 375 376