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