1 /*$Id: tr.c,v 1.114 2000/04/09 04:38:40 bsmith Exp bsmith $*/ 2 3 #include "src/snes/impls/tr/tr.h" /*I "snes.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 __FUNC__ 10 #define __FUNC__ /*<a name=""></a>*/"SNES_EQ_TR_KSPConverged_Private" 11 int SNES_EQ_TR_KSPConverged_Private(KSP ksp,int n,double 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 double norm; 18 int ierr; 19 20 PetscFunctionBegin; 21 if (snes->ksp_ewconv) { 22 if (!kctx) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"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 PLogInfo(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,&norm);CHKERRQ(ierr); 34 if (norm >= neP->delta) { 35 PLogInfo(snes,"SNES_EQ_TR_KSPConverged_Private: KSP iterations=%d, rnorm=%g\n",n,rnorm); 36 PLogInfo(snes,"SNES_EQ_TR_KSPConverged_Private: Ending linear iteration early, delta=%g, length=%g\n",neP->delta,norm); 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 __FUNC__ 55 #define __FUNC__ /*<a name=""></a>*/"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 double rho,fnorm,gnorm,gpnorm,xnorm,delta,norm,ynorm,norm1; 63 Scalar mone = -1.0,cnorm; 64 KSP ksp; 65 SLES sles; 66 SNESConvergedReason reason; 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 91 if (fnorm < snes->atol) {*its = 0; snes->reason = SNES_CONVERGED_FNORM_ABS; PetscFunctionReturn(0);} 92 93 /* set parameter for default relative tolerance convergence test */ 94 snes->ttol = fnorm*snes->rtol; 95 96 /* Set the stopping criteria to use the More' trick. */ 97 ierr = SNESGetSLES(snes,&sles);CHKERRQ(ierr); 98 ierr = SLESGetKSP(sles,&ksp);CHKERRQ(ierr); 99 ierr = KSPSetConvergenceTest(ksp,SNES_EQ_TR_KSPConverged_Private,(void *)snes);CHKERRQ(ierr); 100 101 for (i=0; i<maxits; i++) { 102 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 103 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 104 105 /* Solve J Y = F, where J is Jacobian matrix */ 106 ierr = SLESSolve(snes->sles,F,Ytmp,&lits);CHKERRQ(ierr); 107 snes->linear_its += lits; 108 PLogInfo(snes,"SNESSolve_EQ_TR: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 109 ierr = VecNorm(Ytmp,NORM_2,&norm);CHKERRQ(ierr); 110 norm1 = norm; 111 while(1) { 112 ierr = VecCopy(Ytmp,Y);CHKERRQ(ierr); 113 norm = norm1; 114 115 /* Scale Y if need be and predict new value of F norm */ 116 if (norm >= delta) { 117 norm = delta/norm; 118 gpnorm = (1.0 - norm)*fnorm; 119 cnorm = norm; 120 PLogInfo(snes,"SNESSolve_EQ_TR: Scaling direction by %g\n",norm); 121 ierr = VecScale(&cnorm,Y);CHKERRQ(ierr); 122 norm = gpnorm; 123 ynorm = delta; 124 } else { 125 gpnorm = 0.0; 126 PLogInfo(snes,"SNESSolve_EQ_TR: Direction is in Trust Region\n"); 127 ynorm = norm; 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 PLogInfo(snes,"SNESSolve_EQ_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 141 PLogInfo(snes,"SNESSolve_EQ_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta); 142 neP->delta = delta; 143 if (rho > neP->sigma) break; 144 PLogInfo(snes,"SNESSolve_EQ_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 breakout = 1; 151 break; 152 } 153 snes->nfailures++; 154 } 155 if (!breakout) { 156 fnorm = gnorm; 157 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 158 snes->iter = i+1; 159 snes->norm = fnorm; 160 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 161 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 162 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 163 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm = || X || */ 164 SNESLogConvHistory(snes,fnorm,lits); 165 SNESMonitor(snes,i+1,fnorm); 166 167 /* Test for convergence */ 168 neP->itflag = 1; 169 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 170 if (reason) { 171 break; 172 } 173 } else { 174 break; 175 } 176 } 177 if (X != snes->vec_sol) { 178 /* Verify solution is in corect location */ 179 ierr = VecCopy(X,snes->vec_sol);CHKERRQ(ierr); 180 snes->vec_sol_always = snes->vec_sol; 181 snes->vec_func_always = snes->vec_func; 182 } 183 if (i == maxits) { 184 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 185 i--; 186 reason = SNES_DIVERGED_MAX_IT; 187 } 188 *its = i+1; 189 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 190 snes->reason = reason; 191 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 /*------------------------------------------------------------*/ 195 #undef __FUNC__ 196 #define __FUNC__ /*<a name=""></a>*/"SNESSetUp_EQ_TR" 197 static int SNESSetUp_EQ_TR(SNES snes) 198 { 199 int ierr; 200 201 PetscFunctionBegin; 202 snes->nwork = 4; 203 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 204 PLogObjectParents(snes,snes->nwork,snes->work); 205 snes->vec_sol_update_always = snes->work[3]; 206 PetscFunctionReturn(0); 207 } 208 /*------------------------------------------------------------*/ 209 #undef __FUNC__ 210 #define __FUNC__ /*<a name=""></a>*/"SNESDestroy_EQ_TR" 211 static int SNESDestroy_EQ_TR(SNES snes) 212 { 213 int ierr; 214 215 PetscFunctionBegin; 216 if (snes->nwork) { 217 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 218 } 219 ierr = PetscFree(snes->data);CHKERRQ(ierr); 220 PetscFunctionReturn(0); 221 } 222 /*------------------------------------------------------------*/ 223 224 #undef __FUNC__ 225 #define __FUNC__ /*<a name=""></a>*/"SNESSetFromOptions_EQ_TR" 226 static int SNESSetFromOptions_EQ_TR(SNES snes) 227 { 228 SNES_EQ_TR *ctx = (SNES_EQ_TR *)snes->data; 229 double tmp; 230 int ierr; 231 PetscTruth flg; 232 233 PetscFunctionBegin; 234 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp,&flg);CHKERRQ(ierr); 235 if (flg) {ctx->mu = tmp;} 236 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp,&flg);CHKERRQ(ierr); 237 if (flg) {ctx->eta = tmp;} 238 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp,&flg);CHKERRQ(ierr); 239 if (flg) {ctx->sigma = tmp;} 240 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp,&flg);CHKERRQ(ierr); 241 if (flg) {ctx->delta0 = tmp;} 242 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp,&flg);CHKERRQ(ierr); 243 if (flg) {ctx->delta1 = tmp;} 244 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp,&flg);CHKERRQ(ierr); 245 if (flg) {ctx->delta2 = tmp;} 246 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp,&flg);CHKERRQ(ierr); 247 if (flg) {ctx->delta3 = tmp;} 248 PetscFunctionReturn(0); 249 } 250 251 #undef __FUNC__ 252 #define __FUNC__ /*<a name=""></a>*/"SNESPrintHelp_EQ_TR" 253 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 254 { 255 SNES_EQ_TR *ctx = (SNES_EQ_TR *)snes->data; 256 int ierr; 257 MPI_Comm comm = snes->comm; 258 259 PetscFunctionBegin; 260 ierr = (*PetscHelpPrintf)(comm," method SNESEQTR (tr) for systems of nonlinear equations:\n");CHKERRQ(ierr); 261 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu);CHKERRQ(ierr); 262 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta);CHKERRQ(ierr); 263 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma);CHKERRQ(ierr); 264 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0);CHKERRQ(ierr); 265 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1);CHKERRQ(ierr); 266 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2);CHKERRQ(ierr); 267 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3);CHKERRQ(ierr); 268 PetscFunctionReturn(0); 269 } 270 271 #undef __FUNC__ 272 #define __FUNC__ /*<a name=""></a>*/"SNESView_EQ_TR" 273 static int SNESView_EQ_TR(SNES snes,Viewer viewer) 274 { 275 SNES_EQ_TR *tr = (SNES_EQ_TR *)snes->data; 276 int ierr; 277 PetscTruth isascii; 278 279 PetscFunctionBegin; 280 ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr); 281 if (isascii) { 282 ierr = ViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 283 ierr = ViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 284 } else { 285 SETERRQ1(1,1,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name); 286 } 287 PetscFunctionReturn(0); 288 } 289 290 /* ---------------------------------------------------------------- */ 291 #undef __FUNC__ 292 #define __FUNC__ /*<a name=""></a>*/"SNESConverged_EQ_TR" 293 /*@C 294 SNESConverged_EQ_TR - Monitors the convergence of the trust region 295 method SNESEQTR for solving systems of nonlinear equations (default). 296 297 Collective on SNES 298 299 Input Parameters: 300 + snes - the SNES context 301 . xnorm - 2-norm of current iterate 302 . pnorm - 2-norm of current step 303 . fnorm - 2-norm of function 304 - dummy - unused context 305 306 Output Parameter: 307 . reason - one of 308 $ SNES_CONVERGED_FNORM_ABS - (fnorm < atol), 309 $ SNES_CONVERGED_PNORM_RELATIVE - (pnorm < xtol*xnorm), 310 $ SNES_CONVERGED_FNORM_RELATIVE - (fnorm < rtol*fnorm0), 311 $ SNES_DIVERGED_FUNCTION_COUNT - (nfct > maxf), 312 $ SNES_DIVERGED_FNORM_NAN - (fnorm == NaN), 313 $ SNES_CONVERGED_TR_DELTA - (delta < xnorm*deltatol), 314 $ SNES_CONVERGED_ITERATING - (otherwise) 315 316 where 317 + delta - trust region paramenter 318 . deltatol - trust region size tolerance, 319 set with SNESSetTrustRegionTolerance() 320 . maxf - maximum number of function evaluations, 321 set with SNESSetTolerances() 322 . nfct - number of function evaluations, 323 . atol - absolute function norm tolerance, 324 set with SNESSetTolerances() 325 - xtol - relative function norm tolerance, 326 set with SNESSetTolerances() 327 328 Level: intermediate 329 330 .keywords: SNES, nonlinear, default, converged, convergence 331 332 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 333 @*/ 334 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,SNESConvergedReason *reason,void *dummy) 335 { 336 SNES_EQ_TR *neP = (SNES_EQ_TR *)snes->data; 337 int ierr; 338 339 PetscFunctionBegin; 340 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 341 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 342 } 343 344 if (fnorm != fnorm) { 345 PLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 346 *reason = SNES_DIVERGED_FNORM_NAN; 347 } else if (neP->delta < xnorm * snes->deltatol) { 348 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 349 *reason = SNES_CONVERGED_TR_DELTA; 350 } else if (neP->itflag) { 351 ierr = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,reason,dummy);CHKERRQ(ierr); 352 } else if (snes->nfuncs > snes->max_funcs) { 353 PLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n",snes->nfuncs,snes->max_funcs); 354 *reason = SNES_DIVERGED_FUNCTION_COUNT; 355 } else { 356 *reason = SNES_CONVERGED_ITERATING; 357 } 358 PetscFunctionReturn(0); 359 } 360 /* ------------------------------------------------------------ */ 361 EXTERN_C_BEGIN 362 #undef __FUNC__ 363 #define __FUNC__ /*<a name=""></a>*/"SNESCreate_EQ_TR" 364 int SNESCreate_EQ_TR(SNES snes) 365 { 366 SNES_EQ_TR *neP; 367 368 PetscFunctionBegin; 369 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 370 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 371 } 372 snes->setup = SNESSetUp_EQ_TR; 373 snes->solve = SNESSolve_EQ_TR; 374 snes->destroy = SNESDestroy_EQ_TR; 375 snes->converged = SNESConverged_EQ_TR; 376 snes->printhelp = SNESPrintHelp_EQ_TR; 377 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 378 snes->view = SNESView_EQ_TR; 379 snes->nwork = 0; 380 381 neP = PetscNew(SNES_EQ_TR);CHKPTRQ(neP); 382 PLogObjectMemory(snes,sizeof(SNES_EQ_TR)); 383 snes->data = (void*)neP; 384 neP->mu = 0.25; 385 neP->eta = 0.75; 386 neP->delta = 0.0; 387 neP->delta0 = 0.2; 388 neP->delta1 = 0.3; 389 neP->delta2 = 0.75; 390 neP->delta3 = 2.0; 391 neP->sigma = 0.0001; 392 neP->itflag = 0; 393 neP->rnorm0 = 0; 394 neP->ttol = 0; 395 PetscFunctionReturn(0); 396 } 397 EXTERN_C_END 398 399