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