1 /*$Id: tr.c,v 1.116 2000/04/18 23:07:54 bsmith Exp balay $*/ 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 __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 /* Verify solution is in corect location */ 178 if (X != snes->vec_sol) { 179 ierr = VecCopy(X,snes->vec_sol);CHKERRQ(ierr); 180 } 181 if (F != snes->vec_func) { 182 ierr = VecCopy(F,snes->vec_func);CHKERRQ(ierr); 183 } 184 snes->vec_sol_always = snes->vec_sol; 185 snes->vec_func_always = snes->vec_func; 186 if (i == maxits) { 187 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 188 i--; 189 reason = SNES_DIVERGED_MAX_IT; 190 } 191 *its = i+1; 192 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 193 snes->reason = reason; 194 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 195 PetscFunctionReturn(0); 196 } 197 /*------------------------------------------------------------*/ 198 #undef __FUNC__ 199 #define __FUNC__ /*<a name=""></a>*/"SNESSetUp_EQ_TR" 200 static int SNESSetUp_EQ_TR(SNES snes) 201 { 202 int ierr; 203 204 PetscFunctionBegin; 205 snes->nwork = 4; 206 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 207 PLogObjectParents(snes,snes->nwork,snes->work); 208 snes->vec_sol_update_always = snes->work[3]; 209 PetscFunctionReturn(0); 210 } 211 /*------------------------------------------------------------*/ 212 #undef __FUNC__ 213 #define __FUNC__ /*<a name=""></a>*/"SNESDestroy_EQ_TR" 214 static int SNESDestroy_EQ_TR(SNES snes) 215 { 216 int ierr; 217 218 PetscFunctionBegin; 219 if (snes->nwork) { 220 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 221 } 222 ierr = PetscFree(snes->data);CHKERRQ(ierr); 223 PetscFunctionReturn(0); 224 } 225 /*------------------------------------------------------------*/ 226 227 #undef __FUNC__ 228 #define __FUNC__ /*<a name=""></a>*/"SNESSetFromOptions_EQ_TR" 229 static int SNESSetFromOptions_EQ_TR(SNES snes) 230 { 231 SNES_EQ_TR *ctx = (SNES_EQ_TR *)snes->data; 232 double tmp; 233 int ierr; 234 PetscTruth flg; 235 236 PetscFunctionBegin; 237 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp,&flg);CHKERRQ(ierr); 238 if (flg) {ctx->mu = tmp;} 239 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp,&flg);CHKERRQ(ierr); 240 if (flg) {ctx->eta = tmp;} 241 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp,&flg);CHKERRQ(ierr); 242 if (flg) {ctx->sigma = tmp;} 243 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp,&flg);CHKERRQ(ierr); 244 if (flg) {ctx->delta0 = tmp;} 245 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp,&flg);CHKERRQ(ierr); 246 if (flg) {ctx->delta1 = tmp;} 247 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp,&flg);CHKERRQ(ierr); 248 if (flg) {ctx->delta2 = tmp;} 249 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp,&flg);CHKERRQ(ierr); 250 if (flg) {ctx->delta3 = tmp;} 251 PetscFunctionReturn(0); 252 } 253 254 #undef __FUNC__ 255 #define __FUNC__ /*<a name=""></a>*/"SNESPrintHelp_EQ_TR" 256 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 257 { 258 SNES_EQ_TR *ctx = (SNES_EQ_TR *)snes->data; 259 int ierr; 260 MPI_Comm comm = snes->comm; 261 262 PetscFunctionBegin; 263 ierr = (*PetscHelpPrintf)(comm," method SNESEQTR (tr) for systems of nonlinear equations:\n");CHKERRQ(ierr); 264 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu);CHKERRQ(ierr); 265 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta);CHKERRQ(ierr); 266 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma);CHKERRQ(ierr); 267 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0);CHKERRQ(ierr); 268 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1);CHKERRQ(ierr); 269 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2);CHKERRQ(ierr); 270 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3);CHKERRQ(ierr); 271 PetscFunctionReturn(0); 272 } 273 274 #undef __FUNC__ 275 #define __FUNC__ /*<a name=""></a>*/"SNESView_EQ_TR" 276 static int SNESView_EQ_TR(SNES snes,Viewer viewer) 277 { 278 SNES_EQ_TR *tr = (SNES_EQ_TR *)snes->data; 279 int ierr; 280 PetscTruth isascii; 281 282 PetscFunctionBegin; 283 ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr); 284 if (isascii) { 285 ierr = ViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 286 ierr = ViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 287 } else { 288 SETERRQ1(1,1,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name); 289 } 290 PetscFunctionReturn(0); 291 } 292 293 /* ---------------------------------------------------------------- */ 294 #undef __FUNC__ 295 #define __FUNC__ /*<a name=""></a>*/"SNESConverged_EQ_TR" 296 /*@C 297 SNESConverged_EQ_TR - Monitors the convergence of the trust region 298 method SNESEQTR for solving systems of nonlinear equations (default). 299 300 Collective on SNES 301 302 Input Parameters: 303 + snes - the SNES context 304 . xnorm - 2-norm of current iterate 305 . pnorm - 2-norm of current step 306 . fnorm - 2-norm of function 307 - dummy - unused context 308 309 Output Parameter: 310 . reason - one of 311 $ SNES_CONVERGED_FNORM_ABS - (fnorm < atol), 312 $ SNES_CONVERGED_PNORM_RELATIVE - (pnorm < xtol*xnorm), 313 $ SNES_CONVERGED_FNORM_RELATIVE - (fnorm < rtol*fnorm0), 314 $ SNES_DIVERGED_FUNCTION_COUNT - (nfct > maxf), 315 $ SNES_DIVERGED_FNORM_NAN - (fnorm == NaN), 316 $ SNES_CONVERGED_TR_DELTA - (delta < xnorm*deltatol), 317 $ SNES_CONVERGED_ITERATING - (otherwise) 318 319 where 320 + delta - trust region paramenter 321 . deltatol - trust region size tolerance, 322 set with SNESSetTrustRegionTolerance() 323 . maxf - maximum number of function evaluations, 324 set with SNESSetTolerances() 325 . nfct - number of function evaluations, 326 . atol - absolute function norm tolerance, 327 set with SNESSetTolerances() 328 - xtol - relative function norm tolerance, 329 set with SNESSetTolerances() 330 331 Level: intermediate 332 333 .keywords: SNES, nonlinear, default, converged, convergence 334 335 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 336 @*/ 337 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,SNESConvergedReason *reason,void *dummy) 338 { 339 SNES_EQ_TR *neP = (SNES_EQ_TR *)snes->data; 340 int ierr; 341 342 PetscFunctionBegin; 343 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 344 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 345 } 346 347 if (fnorm != fnorm) { 348 PLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 349 *reason = SNES_DIVERGED_FNORM_NAN; 350 } else if (neP->delta < xnorm * snes->deltatol) { 351 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 352 *reason = SNES_CONVERGED_TR_DELTA; 353 } else if (neP->itflag) { 354 ierr = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,reason,dummy);CHKERRQ(ierr); 355 } else if (snes->nfuncs > snes->max_funcs) { 356 PLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n",snes->nfuncs,snes->max_funcs); 357 *reason = SNES_DIVERGED_FUNCTION_COUNT; 358 } else { 359 *reason = SNES_CONVERGED_ITERATING; 360 } 361 PetscFunctionReturn(0); 362 } 363 /* ------------------------------------------------------------ */ 364 EXTERN_C_BEGIN 365 #undef __FUNC__ 366 #define __FUNC__ /*<a name=""></a>*/"SNESCreate_EQ_TR" 367 int SNESCreate_EQ_TR(SNES snes) 368 { 369 SNES_EQ_TR *neP; 370 371 PetscFunctionBegin; 372 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 373 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 374 } 375 snes->setup = SNESSetUp_EQ_TR; 376 snes->solve = SNESSolve_EQ_TR; 377 snes->destroy = SNESDestroy_EQ_TR; 378 snes->converged = SNESConverged_EQ_TR; 379 snes->printhelp = SNESPrintHelp_EQ_TR; 380 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 381 snes->view = SNESView_EQ_TR; 382 snes->nwork = 0; 383 384 neP = PetscNew(SNES_EQ_TR);CHKPTRQ(neP); 385 PLogObjectMemory(snes,sizeof(SNES_EQ_TR)); 386 snes->data = (void*)neP; 387 neP->mu = 0.25; 388 neP->eta = 0.75; 389 neP->delta = 0.0; 390 neP->delta0 = 0.2; 391 neP->delta1 = 0.3; 392 neP->delta2 = 0.75; 393 neP->delta3 = 2.0; 394 neP->sigma = 0.0001; 395 neP->itflag = 0; 396 neP->rnorm0 = 0; 397 neP->ttol = 0; 398 PetscFunctionReturn(0); 399 } 400 EXTERN_C_END 401 402