1 #ifndef lint 2 static char vcid[] = "$Id: tr.c,v 1.55 1996/08/12 03:43:01 bsmith 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 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 if (snes->ksp_ewconv) { 23 if (!kctx) SETERRQ(1,"SNES_KSP_EW_Converged_Private:Convergence context does not exist"); 24 if (n == 0) SNES_KSP_EW_ComputeRelativeTolerance_Private(snes,ksp); 25 kctx->lresid_last = rnorm; 26 } 27 convinfo = KSPDefaultConverged(ksp,n,rnorm,ctx); 28 if (convinfo) { 29 PLogInfo(snes,"SNES: KSP iterations=%d, rnorm=%g\n",n,rnorm); 30 return convinfo; 31 } 32 33 /* Determine norm of solution */ 34 ierr = KSPBuildSolution(ksp,0,&x); CHKERRQ(ierr); 35 ierr = VecNorm(x,NORM_2,&norm); CHKERRQ(ierr); 36 if (norm >= neP->delta) { 37 PLogInfo(snes,"SNES: KSP iterations=%d, rnorm=%g\n",n,rnorm); 38 PLogInfo(snes, 39 "SNES: Ending linear iteration early, delta %g length %g\n",neP->delta,norm); 40 return 1; 41 } 42 return(0); 43 } 44 /* 45 SNESSolve_EQ_TR - Implements Newton's Method with a very simple trust 46 region approach for solving systems of nonlinear equations. 47 48 The basic algorithm is taken from "The Minpack Project", by More', 49 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 50 of Mathematical Software", Wayne Cowell, editor. 51 52 This is intended as a model implementation, since it does not 53 necessarily have many of the bells and whistles of other 54 implementations. 55 */ 56 static int SNESSolve_EQ_TR(SNES snes,int *its) 57 { 58 SNES_TR *neP = (SNES_TR *) snes->data; 59 Vec X, F, Y, G, TMP, Ytmp; 60 int maxits, i, history_len, ierr, lits; 61 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 62 double rho, fnorm, gnorm, gpnorm, xnorm, delta,norm,*history, ynorm; 63 Scalar mone = -1.0,cnorm; 64 KSP ksp; 65 SLES sles; 66 67 history = snes->conv_hist; /* convergence history */ 68 history_len = snes->conv_hist_len; /* convergence history length */ 69 maxits = snes->max_its; /* maximum number of iterations */ 70 X = snes->vec_sol; /* solution vector */ 71 F = snes->vec_func; /* residual vector */ 72 Y = snes->work[0]; /* work vectors */ 73 G = snes->work[1]; 74 Ytmp = snes->work[2]; 75 76 ierr = VecNorm(X,NORM_2,&xnorm); CHKERRQ(ierr); /* xnorm = || X || */ 77 snes->iter = 0; 78 79 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* F(X) */ 80 ierr = VecNorm(F, NORM_2,&fnorm ); CHKERRQ(ierr); /* fnorm <- || F || */ 81 snes->norm = fnorm; 82 if (history && history_len > 0) history[0] = fnorm; 83 delta = neP->delta0*fnorm; 84 neP->delta = delta; 85 SNESMonitor(snes,0,fnorm); 86 87 /* set parameter for default relative tolerance convergence test */ 88 snes->ttol = fnorm*snes->rtol; 89 90 /* Set the stopping criteria to use the More' trick. */ 91 ierr = SNESGetSLES(snes,&sles); CHKERRQ(ierr); 92 ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr); 93 ierr = KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,(void *)snes);CHKERRQ(ierr); 94 95 for ( i=0; i<maxits; i++ ) { 96 snes->iter = i+1; 97 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 98 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 99 ierr = SLESSolve(snes->sles,F,Ytmp,&lits); CHKERRQ(ierr); 100 ierr = VecNorm(Ytmp,NORM_2,&norm); CHKERRQ(ierr); 101 while(1) { 102 ierr = VecCopy(Ytmp,Y); CHKERRQ(ierr); 103 /* Scale Y if need be and predict new value of F norm */ 104 105 if (norm >= delta) { 106 norm = delta/norm; 107 gpnorm = (1.0 - norm)*fnorm; 108 cnorm = norm; 109 PLogInfo(snes, "Scaling direction by %g \n",norm ); 110 ierr = VecScale(&cnorm,Y); CHKERRQ(ierr); 111 norm = gpnorm; 112 ynorm = delta; 113 } else { 114 gpnorm = 0.0; 115 PLogInfo(snes,"Direction is in Trust Region \n" ); 116 ynorm = norm; 117 } 118 ierr = VecAYPX(&mone,X,Y); CHKERRQ(ierr); /* Y <- X + Y */ 119 ierr = VecCopy(X,snes->vec_sol_update_always); CHKERRQ(ierr); 120 ierr = SNESComputeFunction(snes,Y,G); CHKERRQ(ierr); /* F(X) */ 121 ierr = VecNorm(G,NORM_2,&gnorm); CHKERRQ(ierr); /* gnorm <- || g || */ 122 if (fnorm == gpnorm) rho = 0.0; 123 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 124 125 /* Update size of trust region */ 126 if (rho < neP->mu) delta *= neP->delta1; 127 else if (rho < neP->eta) delta *= neP->delta2; 128 else delta *= neP->delta3; 129 PLogInfo(snes,"%d: f_norm=%g, g_norm=%g, ynorm=%g\n", 130 i, fnorm, gnorm, ynorm ); 131 PLogInfo(snes,"gpred=%g, rho=%g, delta=%g,iters=%d\n", 132 gpnorm, rho, delta, lits); 133 neP->delta = delta; 134 if (rho > neP->sigma) break; 135 PLogInfo(snes,"Trying again in smaller region\n"); 136 /* check to see if progress is hopeless */ 137 neP->itflag = 0; 138 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 139 /* We're not progressing, so return with the current iterate */ 140 if (X != snes->vec_sol) { 141 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 142 snes->vec_sol_always = snes->vec_sol; 143 snes->vec_func_always = snes->vec_func; 144 } 145 } 146 snes->nfailures++; 147 } 148 fnorm = gnorm; 149 snes->norm = fnorm; 150 if (history && history_len > i+1) history[i+1] = fnorm; 151 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 152 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 153 VecNorm(X, NORM_2,&xnorm ); /* xnorm = || X || */ 154 SNESMonitor(snes,i+1,fnorm); 155 156 /* Test for convergence */ 157 neP->itflag = 1; 158 if ((*snes->converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 159 /* Verify solution is in corect location */ 160 if (X != snes->vec_sol) { 161 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 162 snes->vec_sol_always = snes->vec_sol; 163 snes->vec_func_always = snes->vec_func; 164 } 165 break; 166 } 167 } 168 if (i == maxits) { 169 PLogInfo(snes,"Maximum number of iterations has been reached: %d\n",maxits); 170 i--; 171 } 172 *its = i+1; 173 return 0; 174 } 175 /*------------------------------------------------------------*/ 176 static int SNESSetUp_EQ_TR( SNES snes ) 177 { 178 int ierr; 179 snes->nwork = 4; 180 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERRQ(ierr); 181 PLogObjectParents(snes,snes->nwork,snes->work); 182 snes->vec_sol_update_always = snes->work[3]; 183 return 0; 184 } 185 /*------------------------------------------------------------*/ 186 static int SNESDestroy_EQ_TR(PetscObject obj ) 187 { 188 SNES snes = (SNES) obj; 189 int ierr; 190 191 if (snes->nwork) { 192 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 193 } 194 PetscFree(snes->data); 195 return 0; 196 } 197 /*------------------------------------------------------------*/ 198 199 static int SNESSetFromOptions_EQ_TR(SNES snes) 200 { 201 SNES_TR *ctx = (SNES_TR *)snes->data; 202 double tmp; 203 int ierr,flg; 204 205 ierr = OptionsGetDouble(snes->prefix,"-mu",&tmp, &flg); CHKERRQ(ierr); 206 if (flg) {ctx->mu = tmp;} 207 ierr = OptionsGetDouble(snes->prefix,"-eta",&tmp, &flg); CHKERRQ(ierr); 208 if (flg) {ctx->eta = tmp;} 209 ierr = OptionsGetDouble(snes->prefix,"-sigma",&tmp, &flg); CHKERRQ(ierr); 210 if (flg) {ctx->sigma = tmp;} 211 ierr = OptionsGetDouble(snes->prefix,"-delta0",&tmp, &flg); CHKERRQ(ierr); 212 if (flg) {ctx->delta0 = tmp;} 213 ierr = OptionsGetDouble(snes->prefix,"-delta1",&tmp, &flg); CHKERRQ(ierr); 214 if (flg) {ctx->delta1 = tmp;} 215 ierr = OptionsGetDouble(snes->prefix,"-delta2",&tmp, &flg); CHKERRQ(ierr); 216 if (flg) {ctx->delta2 = tmp;} 217 ierr = OptionsGetDouble(snes->prefix,"-delta3",&tmp, &flg); CHKERRQ(ierr); 218 if (flg) {ctx->delta3 = tmp;} 219 return 0; 220 } 221 222 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 223 { 224 SNES_TR *ctx = (SNES_TR *)snes->data; 225 226 PetscFPrintf(snes->comm,stdout," method SNES_EQ_TR (tr) for systems of nonlinear equations:\n"); 227 PetscFPrintf(snes->comm,stdout," %ssnes_trust_region_mu <mu> (default %g)\n",p,ctx->mu); 228 PetscFPrintf(snes->comm,stdout," %ssnes_trust_region_eta <eta> (default %g)\n",p,ctx->eta); 229 PetscFPrintf(snes->comm,stdout," %ssnes_trust_region_sigma <sigma> (default %g)\n",p,ctx->sigma); 230 PetscFPrintf(snes->comm,stdout," %ssnes_trust_region_delta0 <delta0> (default %g)\n",p,ctx->delta0); 231 PetscFPrintf(snes->comm,stdout," %ssnes_trust_region_delta1 <delta1> (default %g)\n",p,ctx->delta1); 232 PetscFPrintf(snes->comm,stdout," %ssnes_trust_region_delta2 <delta2> (default %g)\n",p,ctx->delta2); 233 PetscFPrintf(snes->comm,stdout," %ssnes_trust_region_delta3 <delta3> (default %g)\n",p,ctx->delta3); 234 return 0; 235 } 236 237 static int SNESView_EQ_TR(PetscObject obj,Viewer viewer) 238 { 239 SNES snes = (SNES)obj; 240 SNES_TR *tr = (SNES_TR *)snes->data; 241 FILE *fd; 242 int ierr; 243 ViewerType vtype; 244 245 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 246 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 247 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 248 PetscFPrintf(snes->comm,fd," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma); 249 PetscFPrintf(snes->comm,fd," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n", 250 tr->delta0,tr->delta1,tr->delta2,tr->delta3); 251 } 252 return 0; 253 } 254 255 /* ---------------------------------------------------------------- */ 256 /*@ 257 SNESConverged_EQ_TR - Default test for monitoring the convergence of the 258 trust region method SNES_EQ_TR for solving systems of nonlinear equations. 259 260 Input Parameters: 261 . snes - the SNES context 262 . xnorm - 2-norm of current iterate 263 . pnorm - 2-norm of current step 264 . fnorm - 2-norm of function 265 . dummy - unused context 266 267 Returns: 268 $ 1 if ( delta < xnorm*deltatol ), 269 $ 2 if ( fnorm < atol ), 270 $ 3 if ( pnorm < xtol*xnorm ), 271 $ -2 if ( nfct > maxf ), 272 $ -1 if ( delta < xnorm*epsmch ), 273 $ 0 otherwise, 274 275 where 276 $ delta - trust region paramenter 277 $ deltatol - trust region size tolerance, 278 $ set with SNESSetTrustRegionTolerance() 279 $ maxf - maximum number of function evaluations, 280 $ set with SNESSetTolerances() 281 $ nfct - number of function evaluations, 282 $ atol - absolute function norm tolerance, 283 $ set with SNESSetTolerances() 284 $ xtol - relative function norm tolerance, 285 $ set with SNESSetTolerances() 286 287 .keywords: SNES, nonlinear, default, converged, convergence 288 289 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 290 @*/ 291 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,void *dummy) 292 { 293 SNES_TR *neP = (SNES_TR *)snes->data; 294 double epsmch = 1.0e-14; /* This must be fixed */ 295 int info; 296 297 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 298 SETERRQ(1,"SNESConverged_EQ_LS:For SNES_NONLINEAR_EQUATIONS only"); 299 300 if (neP->delta < xnorm * snes->deltatol) { 301 PLogInfo(snes, 302 "SNES: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 303 return 1; 304 } 305 if (neP->itflag) { 306 info = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,dummy); 307 if (info) return info; 308 } 309 if (neP->delta < xnorm * epsmch) { 310 PLogInfo(snes, 311 "SNES: Converged due to trust region param %g < %g * %g\n",neP->delta,xnorm, epsmch); 312 return -1; 313 } 314 return 0; 315 } 316 /* ------------------------------------------------------------ */ 317 int SNESCreate_EQ_TR(SNES snes ) 318 { 319 SNES_TR *neP; 320 321 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 322 SETERRQ(1,"SNESCreate_EQ_TR:For SNES_NONLINEAR_EQUATIONS only"); 323 snes->type = SNES_EQ_TR; 324 snes->setup = SNESSetUp_EQ_TR; 325 snes->solve = SNESSolve_EQ_TR; 326 snes->destroy = SNESDestroy_EQ_TR; 327 snes->converged = SNESConverged_EQ_TR; 328 snes->printhelp = SNESPrintHelp_EQ_TR; 329 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 330 snes->view = SNESView_EQ_TR; 331 snes->nwork = 0; 332 333 neP = PetscNew(SNES_TR); CHKPTRQ(neP); 334 PLogObjectMemory(snes,sizeof(SNES_TR)); 335 snes->data = (void *) neP; 336 neP->mu = 0.25; 337 neP->eta = 0.75; 338 neP->delta = 0.0; 339 neP->delta0 = 0.2; 340 neP->delta1 = 0.3; 341 neP->delta2 = 0.75; 342 neP->delta3 = 2.0; 343 neP->sigma = 0.0001; 344 neP->itflag = 0; 345 neP->rnorm0 = 0; 346 neP->ttol = 0; 347 return 0; 348 } 349