1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: tr.c,v 1.92 1999/02/09 23:25:07 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, history_len, ierr, lits, breakout = 0; 65 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 66 double rho, fnorm, gnorm, gpnorm, xnorm, delta,norm,*history, ynorm,norm1; 67 Scalar mone = -1.0,cnorm; 68 KSP ksp; 69 SLES sles; 70 71 PetscFunctionBegin; 72 history = snes->conv_hist; /* convergence history */ 73 history_len = snes->conv_hist_size; /* convergence history length */ 74 maxits = snes->max_its; /* maximum number of iterations */ 75 X = snes->vec_sol; /* solution vector */ 76 F = snes->vec_func; /* residual vector */ 77 Y = snes->work[0]; /* work vectors */ 78 G = snes->work[1]; 79 Ytmp = snes->work[2]; 80 81 ierr = VecNorm(X,NORM_2,&xnorm); CHKERRQ(ierr); /* xnorm = || X || */ 82 PetscAMSTakeAccess(snes); 83 snes->iter = 0; 84 PetscAMSGrantAccess(snes); 85 86 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* F(X) */ 87 ierr = VecNorm(F, NORM_2,&fnorm ); CHKERRQ(ierr); /* fnorm <- || F || */ 88 PetscAMSTakeAccess(snes); 89 snes->norm = fnorm; 90 PetscAMSGrantAccess(snes); 91 if (history) history[0] = fnorm; 92 delta = neP->delta0*fnorm; 93 neP->delta = delta; 94 SNESMonitor(snes,0,fnorm); 95 96 if (fnorm < snes->atol) {*its = 0; PetscFunctionReturn(0);} 97 98 /* set parameter for default relative tolerance convergence test */ 99 snes->ttol = fnorm*snes->rtol; 100 101 /* Set the stopping criteria to use the More' trick. */ 102 ierr = SNESGetSLES(snes,&sles); CHKERRQ(ierr); 103 ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr); 104 ierr = KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,(void *)snes);CHKERRQ(ierr); 105 106 for ( i=0; i<maxits; i++ ) { 107 PetscAMSTakeAccess(snes); 108 snes->iter = i+1; 109 PetscAMSGrantAccess(snes); 110 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 111 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 112 113 /* Solve J Y = F, where J is Jacobian matrix */ 114 ierr = SLESSolve(snes->sles,F,Ytmp,&lits); CHKERRQ(ierr); 115 snes->linear_its += PetscAbsInt(lits); 116 PLogInfo(snes,"SNESSolve_EQ_TR: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 117 ierr = VecNorm(Ytmp,NORM_2,&norm); CHKERRQ(ierr); 118 norm1 = norm; 119 while(1) { 120 ierr = VecCopy(Ytmp,Y); CHKERRQ(ierr); 121 norm = norm1; 122 123 /* Scale Y if need be and predict new value of F norm */ 124 if (norm >= delta) { 125 norm = delta/norm; 126 gpnorm = (1.0 - norm)*fnorm; 127 cnorm = norm; 128 PLogInfo(snes,"SNESSolve_EQ_TR: Scaling direction by %g\n",norm ); 129 ierr = VecScale(&cnorm,Y); CHKERRQ(ierr); 130 norm = gpnorm; 131 ynorm = delta; 132 } else { 133 gpnorm = 0.0; 134 PLogInfo(snes,"SNESSolve_EQ_TR: Direction is in Trust Region\n" ); 135 ynorm = norm; 136 } 137 ierr = VecAYPX(&mone,X,Y); CHKERRQ(ierr); /* Y <- X - Y */ 138 ierr = VecCopy(X,snes->vec_sol_update_always); CHKERRQ(ierr); 139 ierr = SNESComputeFunction(snes,Y,G); CHKERRQ(ierr); /* F(X) */ 140 ierr = VecNorm(G,NORM_2,&gnorm); CHKERRQ(ierr); /* gnorm <- || g || */ 141 if (fnorm == gpnorm) rho = 0.0; 142 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 143 144 /* Update size of trust region */ 145 if (rho < neP->mu) delta *= neP->delta1; 146 else if (rho < neP->eta) delta *= neP->delta2; 147 else delta *= neP->delta3; 148 PLogInfo(snes,"SNESSolve_EQ_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 149 PLogInfo(snes,"SNESSolve_EQ_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta); 150 neP->delta = delta; 151 if (rho > neP->sigma) break; 152 PLogInfo(snes,"SNESSolve_EQ_TR: Trying again in smaller region\n"); 153 /* check to see if progress is hopeless */ 154 neP->itflag = 0; 155 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 156 /* We're not progressing, so return with the current iterate */ 157 breakout = 1; break; 158 } 159 snes->nfailures++; 160 } 161 if (!breakout) { 162 fnorm = gnorm; 163 PetscAMSTakeAccess(snes); 164 snes->norm = fnorm; 165 PetscAMSGrantAccess(snes); 166 if (history && history_len > i+1) history[i+1] = fnorm; 167 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 168 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 169 VecNorm(X, NORM_2,&xnorm ); /* xnorm = || X || */ 170 SNESMonitor(snes,i+1,fnorm); 171 172 /* Test for convergence */ 173 neP->itflag = 1; 174 if ((*snes->converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 175 break; 176 } 177 } else { 178 break; 179 } 180 } 181 if (X != snes->vec_sol) { 182 /* Verify solution is in corect location */ 183 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 184 snes->vec_sol_always = snes->vec_sol; 185 snes->vec_func_always = snes->vec_func; 186 } 187 if (i == maxits) { 188 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 189 i--; 190 } 191 if (history) snes->conv_act_size = (history_len < i+1) ? history_len : i+1; 192 *its = i+1; 193 PetscFunctionReturn(0); 194 } 195 /*------------------------------------------------------------*/ 196 #undef __FUNC__ 197 #define __FUNC__ "SNESSetUp_EQ_TR" 198 static int SNESSetUp_EQ_TR(SNES snes) 199 { 200 int ierr; 201 202 PetscFunctionBegin; 203 snes->nwork = 4; 204 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERRQ(ierr); 205 PLogObjectParents(snes,snes->nwork,snes->work); 206 snes->vec_sol_update_always = snes->work[3]; 207 PetscFunctionReturn(0); 208 } 209 /*------------------------------------------------------------*/ 210 #undef __FUNC__ 211 #define __FUNC__ "SNESDestroy_EQ_TR" 212 static int SNESDestroy_EQ_TR(SNES snes ) 213 { 214 int ierr; 215 216 PetscFunctionBegin; 217 if (snes->nwork) { 218 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 219 } 220 PetscFree(snes->data); 221 PetscFunctionReturn(0); 222 } 223 /*------------------------------------------------------------*/ 224 225 #undef __FUNC__ 226 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 227 static int SNESSetFromOptions_EQ_TR(SNES snes) 228 { 229 SNES_TR *ctx = (SNES_TR *)snes->data; 230 double tmp; 231 int ierr,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__ "SNESPrintHelp_EQ_TR" 253 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 254 { 255 SNES_TR *ctx = (SNES_TR *)snes->data; 256 257 PetscFunctionBegin; 258 PetscFPrintf(snes->comm,stdout," method SNES_EQ_TR (tr) for systems of nonlinear equations:\n"); 259 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu); 260 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta); 261 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma); 262 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0); 263 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1); 264 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2); 265 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3); 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 ViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma); 281 ViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3); 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 Returns: 305 + 1 if ( delta < xnorm*deltatol ), 306 . 2 if ( fnorm < atol ), 307 . 3 if ( pnorm < xtol*xnorm ), 308 . -2 if ( nfct > maxf ), 309 . -1 if ( delta < xnorm*epsmch ), 310 - 0 otherwise 311 312 where 313 + delta - trust region paramenter 314 . deltatol - trust region size tolerance, 315 set with SNESSetTrustRegionTolerance() 316 . maxf - maximum number of function evaluations, 317 set with SNESSetTolerances() 318 . nfct - number of function evaluations, 319 . atol - absolute function norm tolerance, 320 set with SNESSetTolerances() 321 - xtol - relative function norm tolerance, 322 set with SNESSetTolerances() 323 324 .keywords: SNES, nonlinear, default, converged, convergence 325 326 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 327 @*/ 328 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,void *dummy) 329 { 330 SNES_TR *neP = (SNES_TR *)snes->data; 331 double epsmch = 1.0e-14; /* This must be fixed */ 332 int info; 333 334 PetscFunctionBegin; 335 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 336 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 337 } 338 339 if (fnorm != fnorm) { 340 PLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 341 PetscFunctionReturn(-3); 342 } 343 if (neP->delta < xnorm * snes->deltatol) { 344 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n", 345 neP->delta,xnorm,snes->deltatol); 346 PetscFunctionReturn(1); 347 } 348 if (neP->itflag) { 349 info = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,dummy); 350 if (info) PetscFunctionReturn(info); 351 } else if (snes->nfuncs > snes->max_funcs) { 352 PLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n", 353 snes->nfuncs, snes->max_funcs ); 354 PetscFunctionReturn(-2); 355 } 356 if (neP->delta < xnorm * epsmch) { 357 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g < %g * %g\n", 358 neP->delta,xnorm, epsmch); 359 PetscFunctionReturn(-1); 360 } 361 PetscFunctionReturn(0); 362 } 363 /* ------------------------------------------------------------ */ 364 EXTERN_C_BEGIN 365 #undef __FUNC__ 366 #define __FUNC__ "SNESCreate_EQ_TR" 367 int SNESCreate_EQ_TR(SNES snes ) 368 { 369 SNES_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_TR); CHKPTRQ(neP); 385 PLogObjectMemory(snes,sizeof(SNES_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