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