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