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