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