1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: tr.c,v 1.101 1999/09/27 21:31:45 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, 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 SNESConvergedReason reason; 71 72 PetscFunctionBegin; 73 maxits = snes->max_its; /* maximum number of iterations */ 74 X = snes->vec_sol; /* solution vector */ 75 F = snes->vec_func; /* residual vector */ 76 Y = snes->work[0]; /* work vectors */ 77 G = snes->work[1]; 78 Ytmp = snes->work[2]; 79 80 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm = || X || */ 81 82 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); /* F(X) */ 83 ierr = VecNorm(F, NORM_2,&fnorm );CHKERRQ(ierr); /* fnorm <- || F || */ 84 ierr = PetscAMSTakeAccess(snes);CHKERRQ(ierr); 85 snes->norm = fnorm; 86 snes->iter = 0; 87 ierr = PetscAMSGrantAccess(snes);CHKERRQ(ierr); 88 delta = neP->delta0*fnorm; 89 neP->delta = delta; 90 SNESLogConvHistory(snes,fnorm,0); 91 SNESMonitor(snes,0,fnorm); 92 93 if (fnorm < snes->atol) {*its = 0; snes->reason = SNES_CONVERGED_FNORM_ABS; 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 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 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 150 if (reason) { 151 /* We're not progressing, so return with the current iterate */ 152 breakout = 1; 153 break; 154 } 155 snes->nfailures++; 156 } 157 if (!breakout) { 158 fnorm = gnorm; 159 ierr = PetscAMSTakeAccess(snes);CHKERRQ(ierr); 160 snes->iter = i+1; 161 snes->norm = fnorm; 162 ierr = PetscAMSGrantAccess(snes);CHKERRQ(ierr); 163 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 164 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 165 VecNorm(X, NORM_2,&xnorm ); /* xnorm = || X || */ 166 SNESLogConvHistory(snes,fnorm,lits); 167 SNESMonitor(snes,i+1,fnorm); 168 169 /* Test for convergence */ 170 neP->itflag = 1; 171 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 172 if (reason) { 173 break; 174 } 175 } else { 176 break; 177 } 178 } 179 if (X != snes->vec_sol) { 180 /* Verify solution is in corect location */ 181 ierr = VecCopy(X,snes->vec_sol);CHKERRQ(ierr); 182 snes->vec_sol_always = snes->vec_sol; 183 snes->vec_func_always = snes->vec_func; 184 } 185 if (i == maxits) { 186 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 187 i--; 188 reason = SNES_DIVERGED_MAX_IT; 189 } 190 *its = i+1; 191 ierr = PetscAMSTakeAccess(snes);CHKERRQ(ierr); 192 snes->reason = reason; 193 ierr = PetscAMSGrantAccess(snes);CHKERRQ(ierr); 194 PetscFunctionReturn(0); 195 } 196 /*------------------------------------------------------------*/ 197 #undef __FUNC__ 198 #define __FUNC__ "SNESSetUp_EQ_TR" 199 static int SNESSetUp_EQ_TR(SNES snes) 200 { 201 int ierr; 202 203 PetscFunctionBegin; 204 snes->nwork = 4; 205 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work );CHKERRQ(ierr); 206 PLogObjectParents(snes,snes->nwork,snes->work); 207 snes->vec_sol_update_always = snes->work[3]; 208 PetscFunctionReturn(0); 209 } 210 /*------------------------------------------------------------*/ 211 #undef __FUNC__ 212 #define __FUNC__ "SNESDestroy_EQ_TR" 213 static int SNESDestroy_EQ_TR(SNES snes ) 214 { 215 int ierr; 216 217 PetscFunctionBegin; 218 if (snes->nwork) { 219 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 220 } 221 ierr = PetscFree(snes->data);CHKERRQ(ierr); 222 PetscFunctionReturn(0); 223 } 224 /*------------------------------------------------------------*/ 225 226 #undef __FUNC__ 227 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 228 static int SNESSetFromOptions_EQ_TR(SNES snes) 229 { 230 SNES_TR *ctx = (SNES_TR *)snes->data; 231 double tmp; 232 int ierr,flg; 233 234 PetscFunctionBegin; 235 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp, &flg);CHKERRQ(ierr); 236 if (flg) {ctx->mu = tmp;} 237 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp, &flg);CHKERRQ(ierr); 238 if (flg) {ctx->eta = tmp;} 239 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp, &flg);CHKERRQ(ierr); 240 if (flg) {ctx->sigma = tmp;} 241 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp, &flg);CHKERRQ(ierr); 242 if (flg) {ctx->delta0 = tmp;} 243 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp, &flg);CHKERRQ(ierr); 244 if (flg) {ctx->delta1 = tmp;} 245 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp, &flg);CHKERRQ(ierr); 246 if (flg) {ctx->delta2 = tmp;} 247 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp, &flg);CHKERRQ(ierr); 248 if (flg) {ctx->delta3 = tmp;} 249 PetscFunctionReturn(0); 250 } 251 252 #undef __FUNC__ 253 #define __FUNC__ "SNESPrintHelp_EQ_TR" 254 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 255 { 256 SNES_TR *ctx = (SNES_TR *)snes->data; 257 int ierr; 258 MPI_Comm comm = snes->comm; 259 260 PetscFunctionBegin; 261 ierr = (*PetscHelpPrintf)(comm," method SNES_EQ_TR (tr) for systems of nonlinear equations:\n");CHKERRQ(ierr); 262 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu);CHKERRQ(ierr); 263 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta);CHKERRQ(ierr); 264 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma);CHKERRQ(ierr); 265 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0);CHKERRQ(ierr); 266 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1);CHKERRQ(ierr); 267 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2);CHKERRQ(ierr); 268 ierr = (*PetscHelpPrintf)(comm," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3);CHKERRQ(ierr); 269 PetscFunctionReturn(0); 270 } 271 272 #undef __FUNC__ 273 #define __FUNC__ "SNESView_EQ_TR" 274 static int SNESView_EQ_TR(SNES snes,Viewer viewer) 275 { 276 SNES_TR *tr = (SNES_TR *)snes->data; 277 int ierr; 278 279 PetscFunctionBegin; 280 if (PetscTypeCompare(viewer,ASCII_VIEWER)) { 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 SETERRQ(1,1,"Viewer type not supported for this object"); 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 SNES_EQ_TR 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_TR *neP = (SNES_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_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 EXTERN_C_END 397 398