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