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