1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: tr.c,v 1.93 1999/02/09 23:31:40 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, history_len, ierr, lits, breakout = 0; 65 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 66 double rho, fnorm, gnorm, gpnorm, xnorm, delta,norm,*history, ynorm,norm1; 67 Scalar mone = -1.0,cnorm; 68 KSP ksp; 69 SLES sles; 70 71 PetscFunctionBegin; 72 history = snes->conv_hist; /* convergence history */ 73 history_len = snes->conv_hist_size; /* convergence history length */ 74 maxits = snes->max_its; /* maximum number of iterations */ 75 X = snes->vec_sol; /* solution vector */ 76 F = snes->vec_func; /* residual vector */ 77 Y = snes->work[0]; /* work vectors */ 78 G = snes->work[1]; 79 Ytmp = snes->work[2]; 80 81 ierr = VecNorm(X,NORM_2,&xnorm); CHKERRQ(ierr); /* xnorm = || X || */ 82 83 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* F(X) */ 84 ierr = VecNorm(F, NORM_2,&fnorm ); CHKERRQ(ierr); /* fnorm <- || F || */ 85 PetscAMSTakeAccess(snes); 86 snes->norm = fnorm; 87 snes->iter = 0; 88 PetscAMSGrantAccess(snes); 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 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 PetscAMSTakeAccess(snes); 159 snes->iter = i+1; 160 snes->norm = fnorm; 161 PetscAMSGrantAccess(snes); 162 if (history && history_len > i+1) history[i+1] = fnorm; 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 SNESMonitor(snes,i+1,fnorm); 167 168 /* Test for convergence */ 169 neP->itflag = 1; 170 if ((*snes->converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 171 break; 172 } 173 } else { 174 break; 175 } 176 } 177 if (X != snes->vec_sol) { 178 /* Verify solution is in corect location */ 179 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 180 snes->vec_sol_always = snes->vec_sol; 181 snes->vec_func_always = snes->vec_func; 182 } 183 if (i == maxits) { 184 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 185 i--; 186 } 187 if (history) snes->conv_act_size = (history_len < i+1) ? history_len : i+1; 188 *its = i+1; 189 PetscFunctionReturn(0); 190 } 191 /*------------------------------------------------------------*/ 192 #undef __FUNC__ 193 #define __FUNC__ "SNESSetUp_EQ_TR" 194 static int SNESSetUp_EQ_TR(SNES snes) 195 { 196 int ierr; 197 198 PetscFunctionBegin; 199 snes->nwork = 4; 200 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERRQ(ierr); 201 PLogObjectParents(snes,snes->nwork,snes->work); 202 snes->vec_sol_update_always = snes->work[3]; 203 PetscFunctionReturn(0); 204 } 205 /*------------------------------------------------------------*/ 206 #undef __FUNC__ 207 #define __FUNC__ "SNESDestroy_EQ_TR" 208 static int SNESDestroy_EQ_TR(SNES snes ) 209 { 210 int ierr; 211 212 PetscFunctionBegin; 213 if (snes->nwork) { 214 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 215 } 216 PetscFree(snes->data); 217 PetscFunctionReturn(0); 218 } 219 /*------------------------------------------------------------*/ 220 221 #undef __FUNC__ 222 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 223 static int SNESSetFromOptions_EQ_TR(SNES snes) 224 { 225 SNES_TR *ctx = (SNES_TR *)snes->data; 226 double tmp; 227 int ierr,flg; 228 229 PetscFunctionBegin; 230 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp, &flg); CHKERRQ(ierr); 231 if (flg) {ctx->mu = tmp;} 232 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp, &flg); CHKERRQ(ierr); 233 if (flg) {ctx->eta = tmp;} 234 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp, &flg); CHKERRQ(ierr); 235 if (flg) {ctx->sigma = tmp;} 236 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp, &flg); CHKERRQ(ierr); 237 if (flg) {ctx->delta0 = tmp;} 238 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp, &flg); CHKERRQ(ierr); 239 if (flg) {ctx->delta1 = tmp;} 240 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp, &flg); CHKERRQ(ierr); 241 if (flg) {ctx->delta2 = tmp;} 242 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp, &flg); CHKERRQ(ierr); 243 if (flg) {ctx->delta3 = tmp;} 244 PetscFunctionReturn(0); 245 } 246 247 #undef __FUNC__ 248 #define __FUNC__ "SNESPrintHelp_EQ_TR" 249 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 250 { 251 SNES_TR *ctx = (SNES_TR *)snes->data; 252 253 PetscFunctionBegin; 254 PetscFPrintf(snes->comm,stdout," method SNES_EQ_TR (tr) for systems of nonlinear equations:\n"); 255 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu); 256 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta); 257 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma); 258 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0); 259 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1); 260 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2); 261 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3); 262 PetscFunctionReturn(0); 263 } 264 265 #undef __FUNC__ 266 #define __FUNC__ "SNESView_EQ_TR" 267 static int SNESView_EQ_TR(SNES snes,Viewer viewer) 268 { 269 SNES_TR *tr = (SNES_TR *)snes->data; 270 int ierr; 271 ViewerType vtype; 272 273 PetscFunctionBegin; 274 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 275 if (PetscTypeCompare(vtype,ASCII_VIEWER)) { 276 ViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma); 277 ViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3); 278 } else { 279 SETERRQ(1,1,"Viewer type not supported for this object"); 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 Collective on SNES 292 293 Input Parameters: 294 + snes - the SNES context 295 . xnorm - 2-norm of current iterate 296 . pnorm - 2-norm of current step 297 . fnorm - 2-norm of function 298 - dummy - unused context 299 300 Returns: 301 + 1 if ( delta < xnorm*deltatol ), 302 . 2 if ( fnorm < atol ), 303 . 3 if ( pnorm < xtol*xnorm ), 304 . -2 if ( nfct > maxf ), 305 . -1 if ( delta < xnorm*epsmch ), 306 - 0 otherwise 307 308 where 309 + delta - trust region paramenter 310 . deltatol - trust region size tolerance, 311 set with SNESSetTrustRegionTolerance() 312 . maxf - maximum number of function evaluations, 313 set with SNESSetTolerances() 314 . nfct - number of function evaluations, 315 . atol - absolute function norm tolerance, 316 set with SNESSetTolerances() 317 - xtol - relative function norm tolerance, 318 set with SNESSetTolerances() 319 320 .keywords: SNES, nonlinear, default, converged, convergence 321 322 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 323 @*/ 324 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,void *dummy) 325 { 326 SNES_TR *neP = (SNES_TR *)snes->data; 327 double epsmch = 1.0e-14; /* This must be fixed */ 328 int info; 329 330 PetscFunctionBegin; 331 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 332 SETERRQ(PETSC_ERR_ARG_WRONG,0,"For SNES_NONLINEAR_EQUATIONS only"); 333 } 334 335 if (fnorm != fnorm) { 336 PLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 337 PetscFunctionReturn(-3); 338 } 339 if (neP->delta < xnorm * snes->deltatol) { 340 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n", 341 neP->delta,xnorm,snes->deltatol); 342 PetscFunctionReturn(1); 343 } 344 if (neP->itflag) { 345 info = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,dummy); 346 if (info) PetscFunctionReturn(info); 347 } else if (snes->nfuncs > snes->max_funcs) { 348 PLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n", 349 snes->nfuncs, snes->max_funcs ); 350 PetscFunctionReturn(-2); 351 } 352 if (neP->delta < xnorm * epsmch) { 353 PLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g < %g * %g\n", 354 neP->delta,xnorm, epsmch); 355 PetscFunctionReturn(-1); 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