1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: tr.c,v 1.103 1999/10/13 20:38:30 bsmith Exp bsmith $"; 3 #endif 4 5 #include "src/snes/impls/tr/tr.h" /*I "snes.h" I*/ 6 7 /* 8 This convergence test determines if the two norm of the 9 solution lies outside the trust region, if so it halts. 10 */ 11 #undef __FUNC__ 12 #define __FUNC__ "SNES_EQ_TR_KSPConverged_Private" 13 int SNES_EQ_TR_KSPConverged_Private(KSP ksp,int n, double rnorm, void *ctx) 14 { 15 SNES snes = (SNES) ctx; 16 SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx; 17 SNES_EQ_TR *neP = (SNES_EQ_TR*)snes->data; 18 Vec x; 19 double norm; 20 int ierr, convinfo; 21 22 PetscFunctionBegin; 23 if (snes->ksp_ewconv) { 24 if (!kctx) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Eisenstat-Walker onvergence context not created"); 25 if (n == 0) {ierr = SNES_KSP_EW_ComputeRelativeTolerance_Private(snes,ksp);CHKERRQ(ierr);} 26 kctx->lresid_last = rnorm; 27 } 28 convinfo = KSPDefaultConverged(ksp,n,rnorm,ctx); 29 if (convinfo) { 30 PLogInfo(snes,"SNES_EQ_TR_KSPConverged_Private: KSP iterations=%d, rnorm=%g\n",n,rnorm); 31 PetscFunctionReturn(convinfo); 32 } 33 34 /* Determine norm of solution */ 35 ierr = KSPBuildSolution(ksp,0,&x);CHKERRQ(ierr); 36 ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); 37 if (norm >= neP->delta) { 38 PLogInfo(snes,"SNES_EQ_TR_KSPConverged_Private: KSP iterations=%d, rnorm=%g\n",n,rnorm); 39 PLogInfo(snes,"SNES_EQ_TR_KSPConverged_Private: Ending linear iteration early, delta=%g, length=%g\n", 40 neP->delta,norm); 41 PetscFunctionReturn(1); 42 } 43 PetscFunctionReturn(0); 44 } 45 46 /* 47 SNESSolve_EQ_TR - Implements Newton's Method with a very simple trust 48 region approach for solving systems of nonlinear equations. 49 50 The basic algorithm is taken from "The Minpack Project", by More', 51 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 52 of Mathematical Software", Wayne Cowell, editor. 53 54 This is intended as a model implementation, since it does not 55 necessarily have many of the bells and whistles of other 56 implementations. 57 */ 58 #undef __FUNC__ 59 #define __FUNC__ "SNESSolve_EQ_TR" 60 static int SNESSolve_EQ_TR(SNES snes,int *its) 61 { 62 SNES_EQ_TR *neP = (SNES_EQ_TR *) snes->data; 63 Vec X, F, Y, G, TMP, Ytmp; 64 int maxits, i, ierr, lits, breakout = 0; 65 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 66 double rho, fnorm, gnorm, gpnorm, xnorm, delta,norm,ynorm,norm1; 67 Scalar mone = -1.0,cnorm; 68 KSP ksp; 69 SLES sles; 70 SNESConvergedReason reason; 71 72 PetscFunctionBegin; 73 maxits = snes->max_its; /* maximum number of iterations */ 74 X = snes->vec_sol; /* solution vector */ 75 F = snes->vec_func; /* residual vector */ 76 Y = snes->work[0]; /* work vectors */ 77 G = snes->work[1]; 78 Ytmp = snes->work[2]; 79 80 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm = || X || */ 81 82 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); /* F(X) */ 83 ierr = VecNorm(F, NORM_2,&fnorm );CHKERRQ(ierr); /* fnorm <- || F || */ 84 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 85 snes->norm = fnorm; 86 snes->iter = 0; 87 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 88 delta = neP->delta0*fnorm; 89 neP->delta = delta; 90 SNESLogConvHistory(snes,fnorm,0); 91 SNESMonitor(snes,0,fnorm); 92 93 if (fnorm < snes->atol) {*its = 0; snes->reason = SNES_CONVERGED_FNORM_ABS; PetscFunctionReturn(0);} 94 95 /* set parameter for default relative tolerance convergence test */ 96 snes->ttol = fnorm*snes->rtol; 97 98 /* Set the stopping criteria to use the More' trick. */ 99 ierr = SNESGetSLES(snes,&sles);CHKERRQ(ierr); 100 ierr = SLESGetKSP(sles,&ksp);CHKERRQ(ierr); 101 ierr = KSPSetConvergenceTest(ksp,SNES_EQ_TR_KSPConverged_Private,(void *)snes);CHKERRQ(ierr); 102 103 for ( i=0; i<maxits; i++ ) { 104 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 105 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 106 107 /* Solve J Y = F, where J is Jacobian matrix */ 108 ierr = SLESSolve(snes->sles,F,Ytmp,&lits);CHKERRQ(ierr); 109 snes->linear_its += PetscAbsInt(lits); 110 PLogInfo(snes,"SNESSolve_EQ_TR: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 111 ierr = VecNorm(Ytmp,NORM_2,&norm);CHKERRQ(ierr); 112 norm1 = norm; 113 while(1) { 114 ierr = VecCopy(Ytmp,Y);CHKERRQ(ierr); 115 norm = norm1; 116 117 /* Scale Y if need be and predict new value of F norm */ 118 if (norm >= delta) { 119 norm = delta/norm; 120 gpnorm = (1.0 - norm)*fnorm; 121 cnorm = norm; 122 PLogInfo(snes,"SNESSolve_EQ_TR: Scaling direction by %g\n",norm ); 123 ierr = VecScale(&cnorm,Y);CHKERRQ(ierr); 124 norm = gpnorm; 125 ynorm = delta; 126 } else { 127 gpnorm = 0.0; 128 PLogInfo(snes,"SNESSolve_EQ_TR: Direction is in Trust Region\n" ); 129 ynorm = norm; 130 } 131 ierr = VecAYPX(&mone,X,Y);CHKERRQ(ierr); /* Y <- X - Y */ 132 ierr = VecCopy(X,snes->vec_sol_update_always);CHKERRQ(ierr); 133 ierr = SNESComputeFunction(snes,Y,G);CHKERRQ(ierr); /* F(X) */ 134 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); /* gnorm <- || g || */ 135 if (fnorm == gpnorm) rho = 0.0; 136 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 137 138 /* Update size of trust region */ 139 if (rho < neP->mu) delta *= neP->delta1; 140 else if (rho < neP->eta) delta *= neP->delta2; 141 else delta *= neP->delta3; 142 PLogInfo(snes,"SNESSolve_EQ_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 143 PLogInfo(snes,"SNESSolve_EQ_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta); 144 neP->delta = delta; 145 if (rho > neP->sigma) break; 146 PLogInfo(snes,"SNESSolve_EQ_TR: Trying again in smaller region\n"); 147 /* check to see if progress is hopeless */ 148 neP->itflag = 0; 149 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 150 if (reason) { 151 /* We're not progressing, so return with the current iterate */ 152 breakout = 1; 153 break; 154 } 155 snes->nfailures++; 156 } 157 if (!breakout) { 158 fnorm = gnorm; 159 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 160 snes->iter = i+1; 161 snes->norm = fnorm; 162 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 163 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 164 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 165 VecNorm(X, NORM_2,&xnorm ); /* xnorm = || X || */ 166 SNESLogConvHistory(snes,fnorm,lits); 167 SNESMonitor(snes,i+1,fnorm); 168 169 /* Test for convergence */ 170 neP->itflag = 1; 171 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 172 if (reason) { 173 break; 174 } 175 } else { 176 break; 177 } 178 } 179 if (X != snes->vec_sol) { 180 /* Verify solution is in corect location */ 181 ierr = VecCopy(X,snes->vec_sol);CHKERRQ(ierr); 182 snes->vec_sol_always = snes->vec_sol; 183 snes->vec_func_always = snes->vec_func; 184 } 185 if (i == maxits) { 186 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 187 i--; 188 reason = SNES_DIVERGED_MAX_IT; 189 } 190 *its = i+1; 191 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 192 snes->reason = reason; 193 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196 /*------------------------------------------------------------*/ 197 #undef __FUNC__ 198 #define __FUNC__ "SNESSetUp_EQ_TR" 199 static int SNESSetUp_EQ_TR(SNES snes) 200 { 201 int ierr; 202 203 PetscFunctionBegin; 204 snes->nwork = 4; 205 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work );CHKERRQ(ierr); 206 PLogObjectParents(snes,snes->nwork,snes->work); 207 snes->vec_sol_update_always = snes->work[3]; 208 PetscFunctionReturn(0); 209 } 210 /*------------------------------------------------------------*/ 211 #undef __FUNC__ 212 #define __FUNC__ "SNESDestroy_EQ_TR" 213 static int SNESDestroy_EQ_TR(SNES snes ) 214 { 215 int ierr; 216 217 PetscFunctionBegin; 218 if (snes->nwork) { 219 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 220 } 221 ierr = PetscFree(snes->data);CHKERRQ(ierr); 222 PetscFunctionReturn(0); 223 } 224 /*------------------------------------------------------------*/ 225 226 #undef __FUNC__ 227 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 228 static int SNESSetFromOptions_EQ_TR(SNES snes) 229 { 230 SNES_EQ_TR *ctx = (SNES_EQ_TR *)snes->data; 231 double tmp; 232 int ierr,flg; 233 234 PetscFunctionBegin; 235 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp, &flg);CHKERRQ(ierr); 236 if (flg) {ctx->mu = tmp;} 237 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp, &flg);CHKERRQ(ierr); 238 if (flg) {ctx->eta = tmp;} 239 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp, &flg);CHKERRQ(ierr); 240 if (flg) {ctx->sigma = tmp;} 241 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp, &flg);CHKERRQ(ierr); 242 if (flg) {ctx->delta0 = tmp;} 243 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp, &flg);CHKERRQ(ierr); 244 if (flg) {ctx->delta1 = tmp;} 245 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp, &flg);CHKERRQ(ierr); 246 if (flg) {ctx->delta2 = tmp;} 247 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp, &flg);CHKERRQ(ierr); 248 if (flg) {ctx->delta3 = tmp;} 249 PetscFunctionReturn(0); 250 } 251 252 #undef __FUNC__ 253 #define __FUNC__ "SNESPrintHelp_EQ_TR" 254 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 255 { 256 SNES_EQ_TR *ctx = (SNES_EQ_TR *)snes->data; 257 int ierr; 258 MPI_Comm comm = snes->comm; 259 260 PetscFunctionBegin; 261 ierr = (*PetscHelpPrintf)(comm," method SNESEQTR (tr) for systems of nonlinear equations:\n");CHKERRQ(ierr); 262 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu);CHKERRQ(ierr); 263 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta);CHKERRQ(ierr); 264 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma);CHKERRQ(ierr); 265 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0);CHKERRQ(ierr); 266 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1);CHKERRQ(ierr); 267 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2);CHKERRQ(ierr); 268 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3);CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNC__ 273 #define __FUNC__ "SNESView_EQ_TR" 274 static int SNESView_EQ_TR(SNES snes,Viewer viewer) 275 { 276 SNES_EQ_TR *tr = (SNES_EQ_TR *)snes->data; 277 int ierr; 278 PetscTruth isascii; 279 280 PetscFunctionBegin; 281 ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr); 282 if (isascii) { 283 ierr = ViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 284 ierr = ViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 285 } else { 286 SETERRQ1(1,1,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name); 287 } 288 PetscFunctionReturn(0); 289 } 290 291 /* ---------------------------------------------------------------- */ 292 #undef __FUNC__ 293 #define __FUNC__ "SNESConverged_EQ_TR" 294 /*@ 295 SNESConverged_EQ_TR - Monitors the convergence of the trust region 296 method SNESEQTR for solving systems of nonlinear equations (default). 297 298 Collective on SNES 299 300 Input Parameters: 301 + snes - the SNES context 302 . xnorm - 2-norm of current iterate 303 . pnorm - 2-norm of current step 304 . fnorm - 2-norm of function 305 - dummy - unused context 306 307 Output Parameter: 308 . reason - one of 309 $ SNES_CONVERGED_FNORM_ABS - ( fnorm < atol ), 310 $ SNES_CONVERGED_PNORM_RELATIVE - ( pnorm < xtol*xnorm ), 311 $ SNES_CONVERGED_FNORM_RELATIVE - ( fnorm < rtol*fnorm0 ), 312 $ SNES_DIVERGED_FUNCTION_COUNT - ( nfct > maxf ), 313 $ SNES_DIVERGED_FNORM_NAN - ( fnorm == NaN ), 314 $ SNES_CONVERGED_TR_DELTA - ( delta < xnorm*deltatol ), 315 $ SNES_CONVERGED_ITERATING - ( otherwise ) 316 317 where 318 + delta - trust region paramenter 319 . deltatol - trust region size tolerance, 320 set with SNESSetTrustRegionTolerance() 321 . maxf - maximum number of function evaluations, 322 set with SNESSetTolerances() 323 . nfct - number of function evaluations, 324 . atol - absolute function norm tolerance, 325 set with SNESSetTolerances() 326 - xtol - relative function norm tolerance, 327 set with SNESSetTolerances() 328 329 Level: intermediate 330 331 .keywords: SNES, nonlinear, default, converged, convergence 332 333 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 334 @*/ 335 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,SNESConvergedReason *reason,void *dummy) 336 { 337 SNES_EQ_TR *neP = (SNES_EQ_TR *)snes->data; 338 int ierr; 339 340 PetscFunctionBegin; 341 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 342 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 343 } 344 345 if (fnorm != fnorm) { 346 PLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 347 *reason = SNES_DIVERGED_FNORM_NAN; 348 } else if (neP->delta < xnorm * snes->deltatol) { 349 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 350 *reason = SNES_CONVERGED_TR_DELTA; 351 } else if (neP->itflag) { 352 ierr = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,reason,dummy);CHKERRQ(ierr); 353 } else if (snes->nfuncs > snes->max_funcs) { 354 PLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n",snes->nfuncs, snes->max_funcs ); 355 *reason = SNES_DIVERGED_FUNCTION_COUNT; 356 } else { 357 *reason = SNES_CONVERGED_ITERATING; 358 } 359 PetscFunctionReturn(0); 360 } 361 /* ------------------------------------------------------------ */ 362 EXTERN_C_BEGIN 363 #undef __FUNC__ 364 #define __FUNC__ "SNESCreate_EQ_TR" 365 int SNESCreate_EQ_TR(SNES snes ) 366 { 367 SNES_EQ_TR *neP; 368 369 PetscFunctionBegin; 370 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 371 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 372 } 373 snes->setup = SNESSetUp_EQ_TR; 374 snes->solve = SNESSolve_EQ_TR; 375 snes->destroy = SNESDestroy_EQ_TR; 376 snes->converged = SNESConverged_EQ_TR; 377 snes->printhelp = SNESPrintHelp_EQ_TR; 378 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 379 snes->view = SNESView_EQ_TR; 380 snes->nwork = 0; 381 382 neP = PetscNew(SNES_EQ_TR);CHKPTRQ(neP); 383 PLogObjectMemory(snes,sizeof(SNES_EQ_TR)); 384 snes->data = (void *) neP; 385 neP->mu = 0.25; 386 neP->eta = 0.75; 387 neP->delta = 0.0; 388 neP->delta0 = 0.2; 389 neP->delta1 = 0.3; 390 neP->delta2 = 0.75; 391 neP->delta3 = 2.0; 392 neP->sigma = 0.0001; 393 neP->itflag = 0; 394 neP->rnorm0 = 0; 395 neP->ttol = 0; 396 PetscFunctionReturn(0); 397 } 398 EXTERN_C_END 399 400