1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: tr.c,v 1.100 1999/09/02 14:54:04 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_TR_KSPConverged_Private" 13 int SNES_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_TR *neP = (SNES_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_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_TR_KSPConverged_Private: KSP iterations=%d, rnorm=%g\n",n,rnorm); 39 PLogInfo(snes,"SNES_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_TR *neP = (SNES_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 = PetscAMSTakeAccess(snes);CHKERRQ(ierr); 85 snes->norm = fnorm; 86 snes->iter = 0; 87 ierr = PetscAMSGrantAccess(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_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; break; 153 } 154 snes->nfailures++; 155 } 156 if (!breakout) { 157 fnorm = gnorm; 158 ierr = PetscAMSTakeAccess(snes);CHKERRQ(ierr); 159 snes->iter = i+1; 160 snes->norm = fnorm; 161 ierr = PetscAMSGrantAccess(snes);CHKERRQ(ierr); 162 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 163 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 164 VecNorm(X, NORM_2,&xnorm ); /* xnorm = || X || */ 165 SNESLogConvHistory(snes,fnorm,lits); 166 SNESMonitor(snes,i+1,fnorm); 167 168 /* Test for convergence */ 169 neP->itflag = 1; 170 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 171 if (reason) { 172 break; 173 } 174 } else { 175 break; 176 } 177 } 178 if (X != snes->vec_sol) { 179 /* Verify solution is in corect location */ 180 ierr = VecCopy(X,snes->vec_sol);CHKERRQ(ierr); 181 snes->vec_sol_always = snes->vec_sol; 182 snes->vec_func_always = snes->vec_func; 183 } 184 if (i == maxits) { 185 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 186 i--; 187 reason = SNES_DIVERGED_MAX_IT; 188 } 189 *its = i+1; 190 snes->reason = reason; 191 PetscFunctionReturn(0); 192 } 193 /*------------------------------------------------------------*/ 194 #undef __FUNC__ 195 #define __FUNC__ "SNESSetUp_EQ_TR" 196 static int SNESSetUp_EQ_TR(SNES snes) 197 { 198 int ierr; 199 200 PetscFunctionBegin; 201 snes->nwork = 4; 202 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work );CHKERRQ(ierr); 203 PLogObjectParents(snes,snes->nwork,snes->work); 204 snes->vec_sol_update_always = snes->work[3]; 205 PetscFunctionReturn(0); 206 } 207 /*------------------------------------------------------------*/ 208 #undef __FUNC__ 209 #define __FUNC__ "SNESDestroy_EQ_TR" 210 static int SNESDestroy_EQ_TR(SNES snes ) 211 { 212 int ierr; 213 214 PetscFunctionBegin; 215 if (snes->nwork) { 216 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 217 } 218 ierr = PetscFree(snes->data);CHKERRQ(ierr); 219 PetscFunctionReturn(0); 220 } 221 /*------------------------------------------------------------*/ 222 223 #undef __FUNC__ 224 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 225 static int SNESSetFromOptions_EQ_TR(SNES snes) 226 { 227 SNES_TR *ctx = (SNES_TR *)snes->data; 228 double tmp; 229 int ierr,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_TR *ctx = (SNES_TR *)snes->data; 254 int ierr; 255 MPI_Comm comm = snes->comm; 256 257 PetscFunctionBegin; 258 ierr = (*PetscHelpPrintf)(comm," method SNES_EQ_TR (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_TR *tr = (SNES_TR *)snes->data; 274 int ierr; 275 ViewerType vtype; 276 277 PetscFunctionBegin; 278 ierr = ViewerGetType(viewer,&vtype);CHKERRQ(ierr); 279 if (PetscTypeCompare(vtype,ASCII_VIEWER)) { 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 SETERRQ(1,1,"Viewer type not supported for this object"); 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 SNES_EQ_TR 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_TR *neP = (SNES_TR *)snes->data; 335 double epsmch = 1.0e-14; /* This must be fixed */ 336 int ierr; 337 338 PetscFunctionBegin; 339 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 340 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 341 } 342 343 if (fnorm != fnorm) { 344 PLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 345 *reason = SNES_DIVERGED_FNORM_NAN; 346 } else if (neP->delta < xnorm * snes->deltatol) { 347 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 348 *reason = SNES_CONVERGED_TR_DELTA; 349 } else if (neP->itflag) { 350 ierr = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,reason,dummy);CHKERRQ(ierr); 351 } else if (snes->nfuncs > snes->max_funcs) { 352 PLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n",snes->nfuncs, snes->max_funcs ); 353 *reason = SNES_DIVERGED_FUNCTION_COUNT; 354 } else { 355 *reason = SNES_CONVERGED_ITERATING; 356 } 357 PetscFunctionReturn(0); 358 } 359 /* ------------------------------------------------------------ */ 360 EXTERN_C_BEGIN 361 #undef __FUNC__ 362 #define __FUNC__ "SNESCreate_EQ_TR" 363 int SNESCreate_EQ_TR(SNES snes ) 364 { 365 SNES_TR *neP; 366 367 PetscFunctionBegin; 368 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 369 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 370 } 371 snes->setup = SNESSetUp_EQ_TR; 372 snes->solve = SNESSolve_EQ_TR; 373 snes->destroy = SNESDestroy_EQ_TR; 374 snes->converged = SNESConverged_EQ_TR; 375 snes->printhelp = SNESPrintHelp_EQ_TR; 376 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 377 snes->view = SNESView_EQ_TR; 378 snes->nwork = 0; 379 380 neP = PetscNew(SNES_TR);CHKPTRQ(neP); 381 PLogObjectMemory(snes,sizeof(SNES_TR)); 382 snes->data = (void *) neP; 383 neP->mu = 0.25; 384 neP->eta = 0.75; 385 neP->delta = 0.0; 386 neP->delta0 = 0.2; 387 neP->delta1 = 0.3; 388 neP->delta2 = 0.75; 389 neP->delta3 = 2.0; 390 neP->sigma = 0.0001; 391 neP->itflag = 0; 392 neP->rnorm0 = 0; 393 neP->ttol = 0; 394 PetscFunctionReturn(0); 395 } 396 EXTERN_C_END 397 398