1 #ifndef lint 2 static char vcid[] = "$Id: tr.c,v 1.42 1996/01/26 04:35:25 bsmith Exp bsmith $"; 3 #endif 4 5 #include <math.h> 6 #include "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((PetscObject)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((PetscObject)snes,"SNES: KSP iterations=%d, rnorm=%g\n",n,rnorm); 38 PLogInfo((PetscObject)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_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_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 if (snes->monitor) {ierr = (*snes->monitor)(snes,0,fnorm,snes->monP); CHKERRQ(ierr);} 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((PetscObject)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((PetscObject)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((PetscObject)snes,"%d: f_norm=%g, g_norm=%g, ynorm=%g\n", 129 i, fnorm, gnorm, ynorm ); 130 PLogInfo((PetscObject)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((PetscObject)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 if (snes->monitor) {ierr = (*snes->monitor)(snes,i+1,fnorm,snes->monP); CHKERRQ(ierr);} 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((PetscObject)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_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_TR(PetscObject obj ) 186 { 187 SNES snes = (SNES) obj; 188 int ierr; 189 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 190 PetscFree(snes->data); 191 return 0; 192 } 193 /*------------------------------------------------------------*/ 194 195 static int SNESSetFromOptions_TR(SNES snes) 196 { 197 SNES_TR *ctx = (SNES_TR *)snes->data; 198 double tmp; 199 int ierr,flg; 200 201 ierr = OptionsGetDouble(snes->prefix,"-mu",&tmp, &flg); CHKERRQ(ierr); 202 if (flg) {ctx->mu = tmp;} 203 ierr = OptionsGetDouble(snes->prefix,"-eta",&tmp, &flg); CHKERRQ(ierr); 204 if (flg) {ctx->eta = tmp;} 205 ierr = OptionsGetDouble(snes->prefix,"-sigma",&tmp, &flg); CHKERRQ(ierr); 206 if (flg) {ctx->sigma = tmp;} 207 ierr = OptionsGetDouble(snes->prefix,"-delta0",&tmp, &flg); CHKERRQ(ierr); 208 if (flg) {ctx->delta0 = tmp;} 209 ierr = OptionsGetDouble(snes->prefix,"-delta1",&tmp, &flg); CHKERRQ(ierr); 210 if (flg) {ctx->delta1 = tmp;} 211 ierr = OptionsGetDouble(snes->prefix,"-delta2",&tmp, &flg); CHKERRQ(ierr); 212 if (flg) {ctx->delta2 = tmp;} 213 ierr = OptionsGetDouble(snes->prefix,"-delta3",&tmp, &flg); CHKERRQ(ierr); 214 if (flg) {ctx->delta3 = tmp;} 215 return 0; 216 } 217 218 static int SNESPrintHelp_TR(SNES snes,char *p) 219 { 220 SNES_TR *ctx = (SNES_TR *)snes->data; 221 222 MPIU_fprintf(snes->comm,stdout," method tr (system of nonlinear equations):\n"); 223 MPIU_fprintf(snes->comm,stdout," %ssnes_trust_region_mu mu (default %g)\n",p,ctx->mu); 224 MPIU_fprintf(snes->comm,stdout," %ssnes_trust_region_eta eta (default %g)\n",p,ctx->eta); 225 MPIU_fprintf(snes->comm,stdout," %ssnes_trust_region_sigma sigma (default %g)\n",p,ctx->sigma); 226 MPIU_fprintf(snes->comm,stdout," %ssnes_trust_region_delta0 delta0 (default %g)\n",p,ctx->delta0); 227 MPIU_fprintf(snes->comm,stdout," %ssnes_trust_region_delta1 delta1 (default %g)\n",p,ctx->delta1); 228 MPIU_fprintf(snes->comm,stdout," %ssnes_trust_region_delta2 delta2 (default %g)\n",p,ctx->delta2); 229 MPIU_fprintf(snes->comm,stdout," %ssnes_trust_region_delta3 delta3 (default %g)\n",p,ctx->delta3); 230 return 0; 231 } 232 233 static int SNESView_TR(PetscObject obj,Viewer viewer) 234 { 235 SNES snes = (SNES)obj; 236 SNES_TR *tr = (SNES_TR *)snes->data; 237 FILE *fd; 238 int ierr; 239 240 ierr = ViewerFileGetPointer(viewer,&fd); CHKERRQ(ierr); 241 MPIU_fprintf(snes->comm,fd," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma); 242 MPIU_fprintf(snes->comm,fd," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n", 243 tr->delta0,tr->delta1,tr->delta2,tr->delta3); 244 return 0; 245 } 246 247 /* ---------------------------------------------------------------- */ 248 /*@ 249 SNESTrustRegionDefaultConverged - Default test for monitoring the 250 convergence of the trust region method SNES_EQ_TR for solving systems 251 of nonlinear equations. 252 253 Input Parameters: 254 . snes - the SNES context 255 . xnorm - 2-norm of current iterate 256 . pnorm - 2-norm of current step 257 . fnorm - 2-norm of function 258 . dummy - unused context 259 260 Returns: 261 $ 1 if ( delta < xnorm*deltatol ), 262 $ 2 if ( fnorm < atol ), 263 $ 3 if ( pnorm < xtol*xnorm ), 264 $ -2 if ( nfct > maxf ), 265 $ -1 if ( delta < xnorm*epsmch ), 266 $ 0 otherwise, 267 268 where 269 $ delta - trust region paramenter 270 $ deltatol - trust region size tolerance, 271 $ set with SNESSetTrustRegionTolerance() 272 $ maxf - maximum number of function evaluations, 273 $ set with SNESSetMaxFunctionEvaluations() 274 $ nfct - number of function evaluations, 275 $ atol - absolute function norm tolerance, 276 $ set with SNESSetAbsoluteTolerance() 277 $ xtol - relative function norm tolerance, 278 $ set with SNESSetRelativeTolerance() 279 280 .keywords: SNES, nonlinear, default, converged, convergence 281 282 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 283 @*/ 284 int SNESTrustRegionDefaultConverged(SNES snes,double xnorm,double pnorm, 285 double fnorm,void *dummy) 286 { 287 SNES_TR *neP = (SNES_TR *)snes->data; 288 double epsmch = 1.0e-14; /* This must be fixed */ 289 int info; 290 291 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 292 SETERRQ(1,"SNESDefaultConverged:For SNES_NONLINEAR_EQUATIONS only"); 293 294 if (neP->delta < xnorm * snes->deltatol) { 295 PLogInfo((PetscObject)snes, 296 "SNES:Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 297 return 1; 298 } 299 if (neP->itflag) { 300 info = SNESDefaultConverged(snes,xnorm,pnorm,fnorm,dummy); 301 if (info) return info; 302 } 303 if (neP->delta < xnorm * epsmch) { 304 PLogInfo((PetscObject)snes, 305 "SNES: Converged due to trust region param %g < %g * %g\n",neP->delta,xnorm, epsmch); 306 return -1; 307 } 308 return 0; 309 } 310 /* ------------------------------------------------------------ */ 311 int SNESCreate_TR(SNES snes ) 312 { 313 SNES_TR *neP; 314 315 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 316 SETERRQ(1,"SNESCreate_TR:For SNES_NONLINEAR_EQUATIONS only"); 317 snes->type = SNES_EQ_NTR; 318 snes->setup = SNESSetUp_TR; 319 snes->solve = SNESSolve_TR; 320 snes->destroy = SNESDestroy_TR; 321 snes->converged = SNESTrustRegionDefaultConverged; 322 snes->printhelp = SNESPrintHelp_TR; 323 snes->setfromoptions = SNESSetFromOptions_TR; 324 snes->view = SNESView_TR; 325 326 neP = PetscNew(SNES_TR); CHKPTRQ(neP); 327 PLogObjectMemory(snes,sizeof(SNES_TR)); 328 snes->data = (void *) neP; 329 neP->mu = 0.25; 330 neP->eta = 0.75; 331 neP->delta = 0.0; 332 neP->delta0 = 0.2; 333 neP->delta1 = 0.3; 334 neP->delta2 = 0.75; 335 neP->delta3 = 2.0; 336 neP->sigma = 0.0001; 337 neP->itflag = 0; 338 neP->rnorm0 = 0; 339 neP->ttol = 0; 340 return 0; 341 } 342