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