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