1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: tr.c,v 1.79 1997/12/01 01:56:54 bsmith Exp bsmith $"; 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 SNESSolve_EQ_TR - Implements Newton's Method with a very simple trust 49 region approach for solving systems of nonlinear equations. 50 51 The basic algorithm is taken from "The Minpack Project", by More', 52 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 53 of Mathematical Software", Wayne Cowell, editor. 54 55 This is intended as a model implementation, since it does not 56 necessarily have many of the bells and whistles of other 57 implementations. 58 */ 59 #undef __FUNC__ 60 #define __FUNC__ "SNESSolve_EQ_TR" 61 static int SNESSolve_EQ_TR(SNES snes,int *its) 62 { 63 SNES_TR *neP = (SNES_TR *) snes->data; 64 Vec X, F, Y, G, TMP, Ytmp; 65 int maxits, i, history_len, ierr, lits, breakout = 0; 66 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 67 double rho, fnorm, gnorm, gpnorm, xnorm, delta,norm,*history, ynorm,norm1; 68 Scalar mone = -1.0,cnorm; 69 KSP ksp; 70 SLES sles; 71 72 PetscFunctionBegin; 73 history = snes->conv_hist; /* convergence history */ 74 history_len = snes->conv_hist_size; /* convergence history length */ 75 maxits = snes->max_its; /* maximum number of iterations */ 76 X = snes->vec_sol; /* solution vector */ 77 F = snes->vec_func; /* residual vector */ 78 Y = snes->work[0]; /* work vectors */ 79 G = snes->work[1]; 80 Ytmp = snes->work[2]; 81 82 ierr = VecNorm(X,NORM_2,&xnorm); CHKERRQ(ierr); /* xnorm = || X || */ 83 snes->iter = 0; 84 85 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* F(X) */ 86 ierr = VecNorm(F, NORM_2,&fnorm ); CHKERRQ(ierr); /* fnorm <- || F || */ 87 snes->norm = fnorm; 88 if (history) history[0] = fnorm; 89 delta = neP->delta0*fnorm; 90 neP->delta = delta; 91 SNESMonitor(snes,0,fnorm); 92 93 if (fnorm < snes->atol) {*its = 0; 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 snes->iter = i+1; 105 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 106 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 107 108 /* Solve J Y = F, where J is Jacobian matrix */ 109 ierr = SLESSolve(snes->sles,F,Ytmp,&lits); CHKERRQ(ierr); 110 snes->linear_its += PetscAbsInt(lits); 111 PLogInfo(snes,"SNESSolve_EQ_TR: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 112 ierr = VecNorm(Ytmp,NORM_2,&norm); CHKERRQ(ierr); 113 norm1 = norm; 114 while(1) { 115 ierr = VecCopy(Ytmp,Y); CHKERRQ(ierr); 116 norm = norm1; 117 118 /* Scale Y if need be and predict new value of F norm */ 119 if (norm >= delta) { 120 norm = delta/norm; 121 gpnorm = (1.0 - norm)*fnorm; 122 cnorm = norm; 123 PLogInfo(snes,"SNESSolve_EQ_TR: Scaling direction by %g\n",norm ); 124 ierr = VecScale(&cnorm,Y); CHKERRQ(ierr); 125 norm = gpnorm; 126 ynorm = delta; 127 } else { 128 gpnorm = 0.0; 129 PLogInfo(snes,"SNESSolve_EQ_TR: Direction is in Trust Region\n" ); 130 ynorm = norm; 131 } 132 ierr = VecAYPX(&mone,X,Y); CHKERRQ(ierr); /* Y <- X - Y */ 133 ierr = VecCopy(X,snes->vec_sol_update_always); CHKERRQ(ierr); 134 ierr = SNESComputeFunction(snes,Y,G); CHKERRQ(ierr); /* F(X) */ 135 ierr = VecNorm(G,NORM_2,&gnorm); CHKERRQ(ierr); /* gnorm <- || g || */ 136 if (fnorm == gpnorm) rho = 0.0; 137 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 138 139 /* Update size of trust region */ 140 if (rho < neP->mu) delta *= neP->delta1; 141 else if (rho < neP->eta) delta *= neP->delta2; 142 else delta *= neP->delta3; 143 PLogInfo(snes,"SNESSolve_EQ_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 144 PLogInfo(snes,"SNESSolve_EQ_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta); 145 neP->delta = delta; 146 if (rho > neP->sigma) break; 147 PLogInfo(snes,"SNESSolve_EQ_TR: Trying again in smaller region\n"); 148 /* check to see if progress is hopeless */ 149 neP->itflag = 0; 150 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 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 snes->norm = fnorm; 159 if (history && history_len > i+1) history[i+1] = fnorm; 160 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 161 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 162 VecNorm(X, NORM_2,&xnorm ); /* xnorm = || X || */ 163 SNESMonitor(snes,i+1,fnorm); 164 165 /* Test for convergence */ 166 neP->itflag = 1; 167 if ((*snes->converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 168 break; 169 } 170 } else { 171 break; 172 } 173 } 174 if (X != snes->vec_sol) { 175 /* Verify solution is in corect location */ 176 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 177 snes->vec_sol_always = snes->vec_sol; 178 snes->vec_func_always = snes->vec_func; 179 } 180 if (i == maxits) { 181 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 182 i--; 183 } 184 if (history) snes->conv_act_size = (history_len < i+1) ? history_len : i+1; 185 *its = i+1; 186 PetscFunctionReturn(0); 187 } 188 /*------------------------------------------------------------*/ 189 #undef __FUNC__ 190 #define __FUNC__ "SNESSetUp_EQ_TR" 191 static int SNESSetUp_EQ_TR( SNES snes ) 192 { 193 int ierr; 194 195 PetscFunctionBegin; 196 snes->nwork = 4; 197 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERRQ(ierr); 198 PLogObjectParents(snes,snes->nwork,snes->work); 199 snes->vec_sol_update_always = snes->work[3]; 200 PetscFunctionReturn(0); 201 } 202 /*------------------------------------------------------------*/ 203 #undef __FUNC__ 204 #define __FUNC__ "SNESDestroy_EQ_TR" 205 static int SNESDestroy_EQ_TR(PetscObject obj ) 206 { 207 SNES snes = (SNES) obj; 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(PetscObject obj,Viewer viewer) 266 { 267 SNES snes = (SNES)obj; 268 SNES_TR *tr = (SNES_TR *)snes->data; 269 FILE *fd; 270 int ierr; 271 ViewerType vtype; 272 273 PetscFunctionBegin; 274 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 275 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 276 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 277 PetscFPrintf(snes->comm,fd," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma); 278 PetscFPrintf(snes->comm,fd," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n", 279 tr->delta0,tr->delta1,tr->delta2,tr->delta3); 280 } 281 PetscFunctionReturn(0); 282 } 283 284 /* ---------------------------------------------------------------- */ 285 #undef __FUNC__ 286 #define __FUNC__ "SNESConverged_EQ_TR" 287 /*@ 288 SNESConverged_EQ_TR - Monitors the convergence of the trust region 289 method SNES_EQ_TR for solving systems of nonlinear equations (default). 290 291 Input Parameters: 292 . snes - the SNES context 293 . xnorm - 2-norm of current iterate 294 . pnorm - 2-norm of current step 295 . fnorm - 2-norm of function 296 . dummy - unused context 297 298 Returns: 299 $ 1 if ( delta < xnorm*deltatol ), 300 $ 2 if ( fnorm < atol ), 301 $ 3 if ( pnorm < xtol*xnorm ), 302 $ -2 if ( nfct > maxf ), 303 $ -1 if ( delta < xnorm*epsmch ), 304 $ 0 otherwise, 305 306 where 307 $ delta - trust region paramenter 308 $ deltatol - trust region size tolerance, 309 $ set with SNESSetTrustRegionTolerance() 310 $ maxf - maximum number of function evaluations, 311 $ set with SNESSetTolerances() 312 $ nfct - number of function evaluations, 313 $ atol - absolute function norm tolerance, 314 $ set with SNESSetTolerances() 315 $ xtol - relative function norm tolerance, 316 $ set with SNESSetTolerances() 317 318 .keywords: SNES, nonlinear, default, converged, convergence 319 320 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 321 @*/ 322 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,void *dummy) 323 { 324 SNES_TR *neP = (SNES_TR *)snes->data; 325 double epsmch = 1.0e-14; /* This must be fixed */ 326 int info; 327 328 PetscFunctionBegin; 329 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 330 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 331 } 332 333 if (fnorm != fnorm) { 334 PLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 335 PetscFunctionReturn(-3); 336 } 337 if (neP->delta < xnorm * snes->deltatol) { 338 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n", 339 neP->delta,xnorm,snes->deltatol); 340 PetscFunctionReturn(1); 341 } 342 if (neP->itflag) { 343 info = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,dummy); 344 if (info) PetscFunctionReturn(info); 345 } else if (snes->nfuncs > snes->max_funcs) { 346 PLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n", 347 snes->nfuncs, snes->max_funcs ); 348 PetscFunctionReturn(-2); 349 } 350 if (neP->delta < xnorm * epsmch) { 351 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g < %g * %g\n", 352 neP->delta,xnorm, epsmch); 353 PetscFunctionReturn(-1); 354 } 355 PetscFunctionReturn(0); 356 } 357 /* ------------------------------------------------------------ */ 358 #undef __FUNC__ 359 #define __FUNC__ "SNESCreate_EQ_TR" 360 int SNESCreate_EQ_TR(SNES snes ) 361 { 362 SNES_TR *neP; 363 364 PetscFunctionBegin; 365 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 366 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 367 } 368 snes->type = SNES_EQ_TR; 369 snes->setup = SNESSetUp_EQ_TR; 370 snes->solve = SNESSolve_EQ_TR; 371 snes->destroy = SNESDestroy_EQ_TR; 372 snes->converged = SNESConverged_EQ_TR; 373 snes->printhelp = SNESPrintHelp_EQ_TR; 374 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 375 snes->view = SNESView_EQ_TR; 376 snes->nwork = 0; 377 378 neP = PetscNew(SNES_TR); CHKPTRQ(neP); 379 PLogObjectMemory(snes,sizeof(SNES_TR)); 380 snes->data = (void *) neP; 381 neP->mu = 0.25; 382 neP->eta = 0.75; 383 neP->delta = 0.0; 384 neP->delta0 = 0.2; 385 neP->delta1 = 0.3; 386 neP->delta2 = 0.75; 387 neP->delta3 = 2.0; 388 neP->sigma = 0.0001; 389 neP->itflag = 0; 390 neP->rnorm0 = 0; 391 neP->ttol = 0; 392 PetscFunctionReturn(0); 393 } 394