1 #ifndef lint 2 static char vcid[] = "$Id: tr.c,v 1.70 1997/01/21 16:56:39 curfman 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 if (snes->ksp_ewconv) { 25 if (!kctx) SETERRQ(1,0,"Convergence context does not exist"); 26 if (n == 0) SNES_KSP_EW_ComputeRelativeTolerance_Private(snes,ksp); 27 kctx->lresid_last = rnorm; 28 } 29 convinfo = KSPDefaultConverged(ksp,n,rnorm,ctx); 30 if (convinfo) { 31 PLogInfo(snes,"SNES: KSP iterations=%d, rnorm=%g\n",n,rnorm); 32 return convinfo; 33 } 34 35 /* Determine norm of solution */ 36 ierr = KSPBuildSolution(ksp,0,&x); CHKERRQ(ierr); 37 ierr = VecNorm(x,NORM_2,&norm); CHKERRQ(ierr); 38 if (norm >= neP->delta) { 39 PLogInfo(snes,"SNES: KSP iterations=%d, rnorm=%g\n",n,rnorm); 40 PLogInfo(snes, 41 "SNES: Ending linear iteration early, delta=%g, length=%g\n",neP->delta,norm); 42 return 1; 43 } 44 return(0); 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 history = snes->conv_hist; /* convergence history */ 72 history_len = snes->conv_hist_len; /* convergence history length */ 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 snes->iter = 0; 82 83 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* F(X) */ 84 ierr = VecNorm(F, NORM_2,&fnorm ); CHKERRQ(ierr); /* fnorm <- || F || */ 85 snes->norm = fnorm; 86 if (history && history_len > 0) history[0] = fnorm; 87 delta = neP->delta0*fnorm; 88 neP->delta = delta; 89 SNESMonitor(snes,0,fnorm); 90 91 /* set parameter for default relative tolerance convergence test */ 92 snes->ttol = fnorm*snes->rtol; 93 94 /* Set the stopping criteria to use the More' trick. */ 95 ierr = SNESGetSLES(snes,&sles); CHKERRQ(ierr); 96 ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr); 97 ierr = KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,(void *)snes);CHKERRQ(ierr); 98 99 for ( i=0; i<maxits; i++ ) { 100 snes->iter = i+1; 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 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 147 /* We're not progressing, so return with the current iterate */ 148 breakout = 1; break; 149 } 150 snes->nfailures++; 151 } 152 if (!breakout) { 153 fnorm = gnorm; 154 snes->norm = fnorm; 155 if (history && history_len > i+1) history[i+1] = fnorm; 156 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 157 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 158 VecNorm(X, NORM_2,&xnorm ); /* xnorm = || X || */ 159 SNESMonitor(snes,i+1,fnorm); 160 161 /* Test for convergence */ 162 neP->itflag = 1; 163 if ((*snes->converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 164 break; 165 } 166 } else { 167 break; 168 } 169 } 170 if (X != snes->vec_sol) { 171 /* Verify solution is in corect location */ 172 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 173 snes->vec_sol_always = snes->vec_sol; 174 snes->vec_func_always = snes->vec_func; 175 } 176 if (i == maxits) { 177 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 178 i--; 179 } 180 *its = i+1; 181 return 0; 182 } 183 /*------------------------------------------------------------*/ 184 #undef __FUNC__ 185 #define __FUNC__ "SNESSetUp_EQ_TR" 186 static int SNESSetUp_EQ_TR( SNES snes ) 187 { 188 int ierr; 189 snes->nwork = 4; 190 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERRQ(ierr); 191 PLogObjectParents(snes,snes->nwork,snes->work); 192 snes->vec_sol_update_always = snes->work[3]; 193 return 0; 194 } 195 /*------------------------------------------------------------*/ 196 #undef __FUNC__ 197 #define __FUNC__ "SNESDestroy_EQ_TR" 198 static int SNESDestroy_EQ_TR(PetscObject obj ) 199 { 200 SNES snes = (SNES) obj; 201 int ierr; 202 203 if (snes->nwork) { 204 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 205 } 206 PetscFree(snes->data); 207 return 0; 208 } 209 /*------------------------------------------------------------*/ 210 211 #undef __FUNC__ 212 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 213 static int SNESSetFromOptions_EQ_TR(SNES snes) 214 { 215 SNES_TR *ctx = (SNES_TR *)snes->data; 216 double tmp; 217 int ierr,flg; 218 219 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp, &flg); CHKERRQ(ierr); 220 if (flg) {ctx->mu = tmp;} 221 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp, &flg); CHKERRQ(ierr); 222 if (flg) {ctx->eta = tmp;} 223 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp, &flg); CHKERRQ(ierr); 224 if (flg) {ctx->sigma = tmp;} 225 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp, &flg); CHKERRQ(ierr); 226 if (flg) {ctx->delta0 = tmp;} 227 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp, &flg); CHKERRQ(ierr); 228 if (flg) {ctx->delta1 = tmp;} 229 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp, &flg); CHKERRQ(ierr); 230 if (flg) {ctx->delta2 = tmp;} 231 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp, &flg); CHKERRQ(ierr); 232 if (flg) {ctx->delta3 = tmp;} 233 return 0; 234 } 235 236 #undef __FUNC__ 237 #define __FUNC__ "SNESPrintHelp_EQ_TR" 238 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 239 { 240 SNES_TR *ctx = (SNES_TR *)snes->data; 241 242 PetscFPrintf(snes->comm,stdout," method SNES_EQ_TR (tr) for systems of nonlinear equations:\n"); 243 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu); 244 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta); 245 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma); 246 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0); 247 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1); 248 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2); 249 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3); 250 return 0; 251 } 252 253 #undef __FUNC__ 254 #define __FUNC__ "SNESView_EQ_TR" 255 static int SNESView_EQ_TR(PetscObject obj,Viewer viewer) 256 { 257 SNES snes = (SNES)obj; 258 SNES_TR *tr = (SNES_TR *)snes->data; 259 FILE *fd; 260 int ierr; 261 ViewerType vtype; 262 263 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 264 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 265 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 266 PetscFPrintf(snes->comm,fd," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma); 267 PetscFPrintf(snes->comm,fd," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n", 268 tr->delta0,tr->delta1,tr->delta2,tr->delta3); 269 } 270 return 0; 271 } 272 273 /* ---------------------------------------------------------------- */ 274 #undef __FUNC__ 275 #define __FUNC__ "SNESConverged_EQ_TR" 276 /*@ 277 SNESConverged_EQ_TR - Monitors the convergence of the trust region 278 method SNES_EQ_TR for solving systems of nonlinear equations (default). 279 280 Input Parameters: 281 . snes - the SNES context 282 . xnorm - 2-norm of current iterate 283 . pnorm - 2-norm of current step 284 . fnorm - 2-norm of function 285 . dummy - unused context 286 287 Returns: 288 $ 1 if ( delta < xnorm*deltatol ), 289 $ 2 if ( fnorm < atol ), 290 $ 3 if ( pnorm < xtol*xnorm ), 291 $ -2 if ( nfct > maxf ), 292 $ -1 if ( delta < xnorm*epsmch ), 293 $ 0 otherwise, 294 295 where 296 $ delta - trust region paramenter 297 $ deltatol - trust region size tolerance, 298 $ set with SNESSetTrustRegionTolerance() 299 $ maxf - maximum number of function evaluations, 300 $ set with SNESSetTolerances() 301 $ nfct - number of function evaluations, 302 $ atol - absolute function norm tolerance, 303 $ set with SNESSetTolerances() 304 $ xtol - relative function norm tolerance, 305 $ set with SNESSetTolerances() 306 307 .keywords: SNES, nonlinear, default, converged, convergence 308 309 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 310 @*/ 311 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,void *dummy) 312 { 313 SNES_TR *neP = (SNES_TR *)snes->data; 314 double epsmch = 1.0e-14; /* This must be fixed */ 315 int info; 316 317 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 318 SETERRQ(1,0,"For SNES_NONLINEAR_EQUATIONS only"); 319 320 if (neP->delta < xnorm * snes->deltatol) { 321 PLogInfo(snes, 322 "SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 323 return 1; 324 } 325 if (neP->itflag) { 326 info = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,dummy); 327 if (info) return info; 328 } else if (snes->nfuncs > snes->max_funcs) { 329 PLogInfo(snes, 330 "SNES: Exceeded maximum number of function evaluations: %d > %d\n", 331 snes->nfuncs, snes->max_funcs ); 332 return -2; 333 } 334 if (neP->delta < xnorm * epsmch) { 335 PLogInfo(snes, 336 "SNESConverged_EQ_TR: Converged due to trust region param %g < %g * %g\n",neP->delta,xnorm, epsmch); 337 return -1; 338 } 339 return 0; 340 } 341 /* ------------------------------------------------------------ */ 342 #undef __FUNC__ 343 #define __FUNC__ "SNESCreate_EQ_TR" 344 int SNESCreate_EQ_TR(SNES snes ) 345 { 346 SNES_TR *neP; 347 348 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 349 SETERRQ(1,0,"For SNES_NONLINEAR_EQUATIONS only"); 350 snes->type = SNES_EQ_TR; 351 snes->setup = SNESSetUp_EQ_TR; 352 snes->solve = SNESSolve_EQ_TR; 353 snes->destroy = SNESDestroy_EQ_TR; 354 snes->converged = SNESConverged_EQ_TR; 355 snes->printhelp = SNESPrintHelp_EQ_TR; 356 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 357 snes->view = SNESView_EQ_TR; 358 snes->nwork = 0; 359 360 neP = PetscNew(SNES_TR); CHKPTRQ(neP); 361 PLogObjectMemory(snes,sizeof(SNES_TR)); 362 snes->data = (void *) neP; 363 neP->mu = 0.25; 364 neP->eta = 0.75; 365 neP->delta = 0.0; 366 neP->delta0 = 0.2; 367 neP->delta1 = 0.3; 368 neP->delta2 = 0.75; 369 neP->delta3 = 2.0; 370 neP->sigma = 0.0001; 371 neP->itflag = 0; 372 neP->rnorm0 = 0; 373 neP->ttol = 0; 374 return 0; 375 } 376