1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: tr.c,v 1.88 1998/12/03 04:05:34 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 snes->iter = 0; 83 84 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* F(X) */ 85 ierr = VecNorm(F, NORM_2,&fnorm ); CHKERRQ(ierr); /* fnorm <- || F || */ 86 snes->norm = fnorm; 87 if (history) history[0] = fnorm; 88 delta = neP->delta0*fnorm; 89 neP->delta = delta; 90 SNESMonitor(snes,0,fnorm); 91 92 if (fnorm < snes->atol) {*its = 0; PetscFunctionReturn(0);} 93 94 /* set parameter for default relative tolerance convergence test */ 95 snes->ttol = fnorm*snes->rtol; 96 97 /* Set the stopping criteria to use the More' trick. */ 98 ierr = SNESGetSLES(snes,&sles); CHKERRQ(ierr); 99 ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr); 100 ierr = KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,(void *)snes);CHKERRQ(ierr); 101 102 for ( i=0; i<maxits; i++ ) { 103 snes->iter = i+1; 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 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 150 /* We're not progressing, so return with the current iterate */ 151 breakout = 1; break; 152 } 153 snes->nfailures++; 154 } 155 if (!breakout) { 156 fnorm = gnorm; 157 snes->norm = fnorm; 158 if (history && history_len > i+1) history[i+1] = fnorm; 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 VecNorm(X, NORM_2,&xnorm ); /* xnorm = || X || */ 162 SNESMonitor(snes,i+1,fnorm); 163 164 /* Test for convergence */ 165 neP->itflag = 1; 166 if ((*snes->converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 167 break; 168 } 169 } else { 170 break; 171 } 172 } 173 if (X != snes->vec_sol) { 174 /* Verify solution is in corect location */ 175 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 176 snes->vec_sol_always = snes->vec_sol; 177 snes->vec_func_always = snes->vec_func; 178 } 179 if (i == maxits) { 180 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 181 i--; 182 } 183 if (history) snes->conv_act_size = (history_len < i+1) ? history_len : i+1; 184 *its = i+1; 185 PetscFunctionReturn(0); 186 } 187 /*------------------------------------------------------------*/ 188 #undef __FUNC__ 189 #define __FUNC__ "SNESSetUp_EQ_TR" 190 static int SNESSetUp_EQ_TR(SNES snes) 191 { 192 int ierr; 193 194 PetscFunctionBegin; 195 snes->nwork = 4; 196 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERRQ(ierr); 197 PLogObjectParents(snes,snes->nwork,snes->work); 198 snes->vec_sol_update_always = snes->work[3]; 199 PetscFunctionReturn(0); 200 } 201 /*------------------------------------------------------------*/ 202 #undef __FUNC__ 203 #define __FUNC__ "SNESDestroy_EQ_TR" 204 static int SNESDestroy_EQ_TR(SNES snes ) 205 { 206 int ierr; 207 208 PetscFunctionBegin; 209 if (snes->nwork) { 210 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 211 } 212 PetscFree(snes->data); 213 PetscFunctionReturn(0); 214 } 215 /*------------------------------------------------------------*/ 216 217 #undef __FUNC__ 218 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 219 static int SNESSetFromOptions_EQ_TR(SNES snes) 220 { 221 SNES_TR *ctx = (SNES_TR *)snes->data; 222 double tmp; 223 int ierr,flg; 224 225 PetscFunctionBegin; 226 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp, &flg); CHKERRQ(ierr); 227 if (flg) {ctx->mu = tmp;} 228 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp, &flg); CHKERRQ(ierr); 229 if (flg) {ctx->eta = tmp;} 230 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp, &flg); CHKERRQ(ierr); 231 if (flg) {ctx->sigma = tmp;} 232 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp, &flg); CHKERRQ(ierr); 233 if (flg) {ctx->delta0 = tmp;} 234 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp, &flg); CHKERRQ(ierr); 235 if (flg) {ctx->delta1 = tmp;} 236 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp, &flg); CHKERRQ(ierr); 237 if (flg) {ctx->delta2 = tmp;} 238 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp, &flg); CHKERRQ(ierr); 239 if (flg) {ctx->delta3 = tmp;} 240 PetscFunctionReturn(0); 241 } 242 243 #undef __FUNC__ 244 #define __FUNC__ "SNESPrintHelp_EQ_TR" 245 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 246 { 247 SNES_TR *ctx = (SNES_TR *)snes->data; 248 249 PetscFunctionBegin; 250 PetscFPrintf(snes->comm,stdout," method SNES_EQ_TR (tr) for systems of nonlinear equations:\n"); 251 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu); 252 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta); 253 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma); 254 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0); 255 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1); 256 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2); 257 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3); 258 PetscFunctionReturn(0); 259 } 260 261 #undef __FUNC__ 262 #define __FUNC__ "SNESView_EQ_TR" 263 static int SNESView_EQ_TR(SNES snes,Viewer viewer) 264 { 265 SNES_TR *tr = (SNES_TR *)snes->data; 266 FILE *fd; 267 int ierr; 268 ViewerType vtype; 269 270 PetscFunctionBegin; 271 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 272 if (PetscTypeCompare(vtype,ASCII_VIEWER)) { 273 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 274 PetscFPrintf(snes->comm,fd," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma); 275 PetscFPrintf(snes->comm,fd," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n", 276 tr->delta0,tr->delta1,tr->delta2,tr->delta3); 277 } else { 278 SETERRQ(1,1,"Viewer type not supported for this object"); 279 } 280 PetscFunctionReturn(0); 281 } 282 283 /* ---------------------------------------------------------------- */ 284 #undef __FUNC__ 285 #define __FUNC__ "SNESConverged_EQ_TR" 286 /*@ 287 SNESConverged_EQ_TR - Monitors the convergence of the trust region 288 method SNES_EQ_TR for solving systems of nonlinear equations (default). 289 290 Collective on SNES 291 292 Input Parameters: 293 + snes - the SNES context 294 . xnorm - 2-norm of current iterate 295 . pnorm - 2-norm of current step 296 . fnorm - 2-norm of function 297 - dummy - unused context 298 299 Returns: 300 + 1 if ( delta < xnorm*deltatol ), 301 . 2 if ( fnorm < atol ), 302 . 3 if ( pnorm < xtol*xnorm ), 303 . -2 if ( nfct > maxf ), 304 . -1 if ( delta < xnorm*epsmch ), 305 - 0 otherwise 306 307 where 308 + delta - trust region paramenter 309 . deltatol - trust region size tolerance, 310 set with SNESSetTrustRegionTolerance() 311 . maxf - maximum number of function evaluations, 312 set with SNESSetTolerances() 313 . nfct - number of function evaluations, 314 . atol - absolute function norm tolerance, 315 set with SNESSetTolerances() 316 - xtol - relative function norm tolerance, 317 set with SNESSetTolerances() 318 319 .keywords: SNES, nonlinear, default, converged, convergence 320 321 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 322 @*/ 323 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,void *dummy) 324 { 325 SNES_TR *neP = (SNES_TR *)snes->data; 326 double epsmch = 1.0e-14; /* This must be fixed */ 327 int info; 328 329 PetscFunctionBegin; 330 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 331 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 332 } 333 334 if (fnorm != fnorm) { 335 PLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 336 PetscFunctionReturn(-3); 337 } 338 if (neP->delta < xnorm * snes->deltatol) { 339 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n", 340 neP->delta,xnorm,snes->deltatol); 341 PetscFunctionReturn(1); 342 } 343 if (neP->itflag) { 344 info = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,dummy); 345 if (info) PetscFunctionReturn(info); 346 } else if (snes->nfuncs > snes->max_funcs) { 347 PLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n", 348 snes->nfuncs, snes->max_funcs ); 349 PetscFunctionReturn(-2); 350 } 351 if (neP->delta < xnorm * epsmch) { 352 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g < %g * %g\n", 353 neP->delta,xnorm, epsmch); 354 PetscFunctionReturn(-1); 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_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_TR); CHKPTRQ(neP); 380 PLogObjectMemory(snes,sizeof(SNES_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