1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: tr.c,v 1.99 1999/06/30 23:54:13 balay 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 71 PetscFunctionBegin; 72 maxits = snes->max_its; /* maximum number of iterations */ 73 X = snes->vec_sol; /* solution vector */ 74 F = snes->vec_func; /* residual vector */ 75 Y = snes->work[0]; /* work vectors */ 76 G = snes->work[1]; 77 Ytmp = snes->work[2]; 78 79 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm = || X || */ 80 81 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); /* F(X) */ 82 ierr = VecNorm(F, NORM_2,&fnorm );CHKERRQ(ierr); /* fnorm <- || F || */ 83 ierr = PetscAMSTakeAccess(snes);CHKERRQ(ierr); 84 snes->norm = fnorm; 85 snes->iter = 0; 86 ierr = PetscAMSGrantAccess(snes);CHKERRQ(ierr); 87 delta = neP->delta0*fnorm; 88 neP->delta = delta; 89 SNESLogConvHistory(snes,fnorm,0); 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 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 104 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 105 106 /* Solve J Y = F, where J is Jacobian matrix */ 107 ierr = SLESSolve(snes->sles,F,Ytmp,&lits);CHKERRQ(ierr); 108 snes->linear_its += PetscAbsInt(lits); 109 PLogInfo(snes,"SNESSolve_EQ_TR: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 110 ierr = VecNorm(Ytmp,NORM_2,&norm);CHKERRQ(ierr); 111 norm1 = norm; 112 while(1) { 113 ierr = VecCopy(Ytmp,Y);CHKERRQ(ierr); 114 norm = norm1; 115 116 /* Scale Y if need be and predict new value of F norm */ 117 if (norm >= delta) { 118 norm = delta/norm; 119 gpnorm = (1.0 - norm)*fnorm; 120 cnorm = norm; 121 PLogInfo(snes,"SNESSolve_EQ_TR: Scaling direction by %g\n",norm ); 122 ierr = VecScale(&cnorm,Y);CHKERRQ(ierr); 123 norm = gpnorm; 124 ynorm = delta; 125 } else { 126 gpnorm = 0.0; 127 PLogInfo(snes,"SNESSolve_EQ_TR: Direction is in Trust Region\n" ); 128 ynorm = norm; 129 } 130 ierr = VecAYPX(&mone,X,Y);CHKERRQ(ierr); /* Y <- X - Y */ 131 ierr = VecCopy(X,snes->vec_sol_update_always);CHKERRQ(ierr); 132 ierr = SNESComputeFunction(snes,Y,G);CHKERRQ(ierr); /* F(X) */ 133 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); /* gnorm <- || g || */ 134 if (fnorm == gpnorm) rho = 0.0; 135 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 136 137 /* Update size of trust region */ 138 if (rho < neP->mu) delta *= neP->delta1; 139 else if (rho < neP->eta) delta *= neP->delta2; 140 else delta *= neP->delta3; 141 PLogInfo(snes,"SNESSolve_EQ_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 142 PLogInfo(snes,"SNESSolve_EQ_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta); 143 neP->delta = delta; 144 if (rho > neP->sigma) break; 145 PLogInfo(snes,"SNESSolve_EQ_TR: Trying again in smaller region\n"); 146 /* check to see if progress is hopeless */ 147 neP->itflag = 0; 148 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 149 /* We're not progressing, so return with the current iterate */ 150 breakout = 1; break; 151 } 152 snes->nfailures++; 153 } 154 if (!breakout) { 155 fnorm = gnorm; 156 ierr = PetscAMSTakeAccess(snes);CHKERRQ(ierr); 157 snes->iter = i+1; 158 snes->norm = fnorm; 159 ierr = PetscAMSGrantAccess(snes);CHKERRQ(ierr); 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 SNESLogConvHistory(snes,fnorm,lits); 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 *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(SNES snes ) 206 { 207 int ierr; 208 209 PetscFunctionBegin; 210 if (snes->nwork) { 211 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 212 } 213 ierr = PetscFree(snes->data);CHKERRQ(ierr); 214 PetscFunctionReturn(0); 215 } 216 /*------------------------------------------------------------*/ 217 218 #undef __FUNC__ 219 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 220 static int SNESSetFromOptions_EQ_TR(SNES snes) 221 { 222 SNES_TR *ctx = (SNES_TR *)snes->data; 223 double tmp; 224 int ierr,flg; 225 226 PetscFunctionBegin; 227 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp, &flg);CHKERRQ(ierr); 228 if (flg) {ctx->mu = tmp;} 229 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp, &flg);CHKERRQ(ierr); 230 if (flg) {ctx->eta = tmp;} 231 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp, &flg);CHKERRQ(ierr); 232 if (flg) {ctx->sigma = tmp;} 233 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp, &flg);CHKERRQ(ierr); 234 if (flg) {ctx->delta0 = tmp;} 235 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp, &flg);CHKERRQ(ierr); 236 if (flg) {ctx->delta1 = tmp;} 237 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp, &flg);CHKERRQ(ierr); 238 if (flg) {ctx->delta2 = tmp;} 239 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp, &flg);CHKERRQ(ierr); 240 if (flg) {ctx->delta3 = tmp;} 241 PetscFunctionReturn(0); 242 } 243 244 #undef __FUNC__ 245 #define __FUNC__ "SNESPrintHelp_EQ_TR" 246 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 247 { 248 SNES_TR *ctx = (SNES_TR *)snes->data; 249 int ierr; 250 MPI_Comm comm = snes->comm; 251 252 PetscFunctionBegin; 253 ierr = (*PetscHelpPrintf)(comm," method SNES_EQ_TR (tr) for systems of nonlinear equations:\n");CHKERRQ(ierr); 254 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu);CHKERRQ(ierr); 255 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta);CHKERRQ(ierr); 256 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma);CHKERRQ(ierr); 257 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0);CHKERRQ(ierr); 258 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1);CHKERRQ(ierr); 259 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2);CHKERRQ(ierr); 260 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3);CHKERRQ(ierr); 261 PetscFunctionReturn(0); 262 } 263 264 #undef __FUNC__ 265 #define __FUNC__ "SNESView_EQ_TR" 266 static int SNESView_EQ_TR(SNES snes,Viewer viewer) 267 { 268 SNES_TR *tr = (SNES_TR *)snes->data; 269 int ierr; 270 ViewerType vtype; 271 272 PetscFunctionBegin; 273 ierr = ViewerGetType(viewer,&vtype);CHKERRQ(ierr); 274 if (PetscTypeCompare(vtype,ASCII_VIEWER)) { 275 ierr = ViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 276 ierr = ViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 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 Level: intermediate 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 EXTERN_C_BEGIN 362 #undef __FUNC__ 363 #define __FUNC__ "SNESCreate_EQ_TR" 364 int SNESCreate_EQ_TR(SNES snes ) 365 { 366 SNES_TR *neP; 367 368 PetscFunctionBegin; 369 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 370 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 371 } 372 snes->setup = SNESSetUp_EQ_TR; 373 snes->solve = SNESSolve_EQ_TR; 374 snes->destroy = SNESDestroy_EQ_TR; 375 snes->converged = SNESConverged_EQ_TR; 376 snes->printhelp = SNESPrintHelp_EQ_TR; 377 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 378 snes->view = SNESView_EQ_TR; 379 snes->nwork = 0; 380 381 neP = PetscNew(SNES_TR);CHKPTRQ(neP); 382 PLogObjectMemory(snes,sizeof(SNES_TR)); 383 snes->data = (void *) neP; 384 neP->mu = 0.25; 385 neP->eta = 0.75; 386 neP->delta = 0.0; 387 neP->delta0 = 0.2; 388 neP->delta1 = 0.3; 389 neP->delta2 = 0.75; 390 neP->delta3 = 2.0; 391 neP->sigma = 0.0001; 392 neP->itflag = 0; 393 neP->rnorm0 = 0; 394 neP->ttol = 0; 395 PetscFunctionReturn(0); 396 } 397 EXTERN_C_END 398 399