1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: tr.c,v 1.75 1997/07/09 20:59:53 balay Exp curfman $"; 3 #endif 4 5 #include <math.h> 6 #include "src/snes/impls/tr/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 #undef __FUNC__ 14 #define __FUNC__ "SNES_TR_KSPConverged_Private" 15 int SNES_TR_KSPConverged_Private(KSP ksp,int n, double rnorm, void *ctx) 16 { 17 SNES snes = (SNES) ctx; 18 SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx; 19 SNES_TR *neP = (SNES_TR*)snes->data; 20 Vec x; 21 double norm; 22 int ierr, convinfo; 23 24 if (snes->ksp_ewconv) { 25 if (!kctx) SETERRQ(1,0,"Convergence context does not exist"); 26 if (n == 0) SNES_KSP_EW_ComputeRelativeTolerance_Private(snes,ksp); 27 kctx->lresid_last = rnorm; 28 } 29 convinfo = KSPDefaultConverged(ksp,n,rnorm,ctx); 30 if (convinfo) { 31 PLogInfo(snes,"SNES: KSP iterations=%d, rnorm=%g\n",n,rnorm); 32 return convinfo; 33 } 34 35 /* Determine norm of solution */ 36 ierr = KSPBuildSolution(ksp,0,&x); CHKERRQ(ierr); 37 ierr = VecNorm(x,NORM_2,&norm); CHKERRQ(ierr); 38 if (norm >= neP->delta) { 39 PLogInfo(snes,"SNES: KSP iterations=%d, rnorm=%g\n",n,rnorm); 40 PLogInfo(snes, 41 "SNES: Ending linear iteration early, delta=%g, length=%g\n",neP->delta,norm); 42 return 1; 43 } 44 return(0); 45 } 46 /* 47 SNESSolve_EQ_TR - Implements Newton's Method with a very simple trust 48 region approach for solving systems of nonlinear equations. 49 50 The basic algorithm is taken from "The Minpack Project", by More', 51 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 52 of Mathematical Software", Wayne Cowell, editor. 53 54 This is intended as a model implementation, since it does not 55 necessarily have many of the bells and whistles of other 56 implementations. 57 */ 58 #undef __FUNC__ 59 #define __FUNC__ "SNESSolve_EQ_TR" 60 static int SNESSolve_EQ_TR(SNES snes,int *its) 61 { 62 SNES_TR *neP = (SNES_TR *) snes->data; 63 Vec X, F, Y, G, TMP, Ytmp; 64 int maxits, i, history_len, ierr, lits, breakout = 0; 65 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 66 double rho, fnorm, gnorm, gpnorm, xnorm, delta,norm,*history, ynorm,norm1; 67 Scalar mone = -1.0,cnorm; 68 KSP ksp; 69 SLES sles; 70 71 history = snes->conv_hist; /* convergence history */ 72 history_len = snes->conv_hist_size; /* convergence history length */ 73 maxits = snes->max_its; /* maximum number of iterations */ 74 X = snes->vec_sol; /* solution vector */ 75 F = snes->vec_func; /* residual vector */ 76 Y = snes->work[0]; /* work vectors */ 77 G = snes->work[1]; 78 Ytmp = snes->work[2]; 79 80 ierr = VecNorm(X,NORM_2,&xnorm); CHKERRQ(ierr); /* xnorm = || X || */ 81 snes->iter = 0; 82 83 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* F(X) */ 84 ierr = VecNorm(F, NORM_2,&fnorm ); CHKERRQ(ierr); /* fnorm <- || F || */ 85 snes->norm = fnorm; 86 if (history) history[0] = fnorm; 87 delta = neP->delta0*fnorm; 88 neP->delta = delta; 89 SNESMonitor(snes,0,fnorm); 90 91 if (fnorm < snes->atol) {*its = 0; return 0;} 92 93 /* set parameter for default relative tolerance convergence test */ 94 snes->ttol = fnorm*snes->rtol; 95 96 /* Set the stopping criteria to use the More' trick. */ 97 ierr = SNESGetSLES(snes,&sles); CHKERRQ(ierr); 98 ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr); 99 ierr = KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,(void *)snes);CHKERRQ(ierr); 100 101 for ( i=0; i<maxits; i++ ) { 102 snes->iter = i+1; 103 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 104 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 105 106 /* Solve J Y = F, where J is Jacobian matrix */ 107 ierr = SLESSolve(snes->sles,F,Ytmp,&lits); CHKERRQ(ierr); 108 snes->linear_its += PetscAbsInt(lits); 109 PLogInfo(snes,"SNESSolve_EQ_TR: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 110 ierr = VecNorm(Ytmp,NORM_2,&norm); CHKERRQ(ierr); 111 norm1 = norm; 112 while(1) { 113 ierr = VecCopy(Ytmp,Y); CHKERRQ(ierr); 114 norm = norm1; 115 116 /* Scale Y if need be and predict new value of F norm */ 117 if (norm >= delta) { 118 norm = delta/norm; 119 gpnorm = (1.0 - norm)*fnorm; 120 cnorm = norm; 121 PLogInfo(snes,"SNESSolve_EQ_TR: Scaling direction by %g\n",norm ); 122 ierr = VecScale(&cnorm,Y); CHKERRQ(ierr); 123 norm = gpnorm; 124 ynorm = delta; 125 } else { 126 gpnorm = 0.0; 127 PLogInfo(snes,"SNESSolve_EQ_TR: Direction is in Trust Region\n" ); 128 ynorm = norm; 129 } 130 ierr = VecAYPX(&mone,X,Y); CHKERRQ(ierr); /* Y <- X - Y */ 131 ierr = VecCopy(X,snes->vec_sol_update_always); CHKERRQ(ierr); 132 ierr = SNESComputeFunction(snes,Y,G); CHKERRQ(ierr); /* F(X) */ 133 ierr = VecNorm(G,NORM_2,&gnorm); CHKERRQ(ierr); /* gnorm <- || g || */ 134 if (fnorm == gpnorm) rho = 0.0; 135 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 136 137 /* Update size of trust region */ 138 if (rho < neP->mu) delta *= neP->delta1; 139 else if (rho < neP->eta) delta *= neP->delta2; 140 else delta *= neP->delta3; 141 PLogInfo(snes,"SNESSolve_EQ_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 142 PLogInfo(snes,"SNESSolve_EQ_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta); 143 neP->delta = delta; 144 if (rho > neP->sigma) break; 145 PLogInfo(snes,"SNESSolve_EQ_TR: Trying again in smaller region\n"); 146 /* check to see if progress is hopeless */ 147 neP->itflag = 0; 148 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 149 /* We're not progressing, so return with the current iterate */ 150 breakout = 1; break; 151 } 152 snes->nfailures++; 153 } 154 if (!breakout) { 155 fnorm = gnorm; 156 snes->norm = fnorm; 157 if (history && history_len > i+1) history[i+1] = fnorm; 158 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 159 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 160 VecNorm(X, NORM_2,&xnorm ); /* xnorm = || X || */ 161 SNESMonitor(snes,i+1,fnorm); 162 163 /* Test for convergence */ 164 neP->itflag = 1; 165 if ((*snes->converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 166 break; 167 } 168 } else { 169 break; 170 } 171 } 172 if (X != snes->vec_sol) { 173 /* Verify solution is in corect location */ 174 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 175 snes->vec_sol_always = snes->vec_sol; 176 snes->vec_func_always = snes->vec_func; 177 } 178 if (i == maxits) { 179 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 180 i--; 181 } 182 if (history) snes->conv_act_size = (history_len < i+1) ? history_len : i+1; 183 *its = i+1; 184 return 0; 185 } 186 /*------------------------------------------------------------*/ 187 #undef __FUNC__ 188 #define __FUNC__ "SNESSetUp_EQ_TR" 189 static int SNESSetUp_EQ_TR( SNES snes ) 190 { 191 int ierr; 192 snes->nwork = 4; 193 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERRQ(ierr); 194 PLogObjectParents(snes,snes->nwork,snes->work); 195 snes->vec_sol_update_always = snes->work[3]; 196 return 0; 197 } 198 /*------------------------------------------------------------*/ 199 #undef __FUNC__ 200 #define __FUNC__ "SNESDestroy_EQ_TR" /* ADIC Ignore */ 201 static int SNESDestroy_EQ_TR(PetscObject obj ) 202 { 203 SNES snes = (SNES) obj; 204 int ierr; 205 206 if (snes->nwork) { 207 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 208 } 209 PetscFree(snes->data); 210 return 0; 211 } 212 /*------------------------------------------------------------*/ 213 214 #undef __FUNC__ 215 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 216 static int SNESSetFromOptions_EQ_TR(SNES snes) 217 { 218 SNES_TR *ctx = (SNES_TR *)snes->data; 219 double tmp; 220 int ierr,flg; 221 222 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp, &flg); CHKERRQ(ierr); 223 if (flg) {ctx->mu = tmp;} 224 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp, &flg); CHKERRQ(ierr); 225 if (flg) {ctx->eta = tmp;} 226 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp, &flg); CHKERRQ(ierr); 227 if (flg) {ctx->sigma = tmp;} 228 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp, &flg); CHKERRQ(ierr); 229 if (flg) {ctx->delta0 = tmp;} 230 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp, &flg); CHKERRQ(ierr); 231 if (flg) {ctx->delta1 = tmp;} 232 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp, &flg); CHKERRQ(ierr); 233 if (flg) {ctx->delta2 = tmp;} 234 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp, &flg); CHKERRQ(ierr); 235 if (flg) {ctx->delta3 = tmp;} 236 return 0; 237 } 238 239 #undef __FUNC__ 240 #define __FUNC__ "SNESPrintHelp_EQ_TR" /* ADIC Ignore */ 241 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 242 { 243 SNES_TR *ctx = (SNES_TR *)snes->data; 244 245 PetscFPrintf(snes->comm,stdout," method SNES_EQ_TR (tr) for systems of nonlinear equations:\n"); 246 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu); 247 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta); 248 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma); 249 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0); 250 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1); 251 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2); 252 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3); 253 return 0; 254 } 255 256 #undef __FUNC__ 257 #define __FUNC__ "SNESView_EQ_TR" /* ADIC Ignore */ 258 static int SNESView_EQ_TR(PetscObject obj,Viewer viewer) 259 { 260 SNES snes = (SNES)obj; 261 SNES_TR *tr = (SNES_TR *)snes->data; 262 FILE *fd; 263 int ierr; 264 ViewerType vtype; 265 266 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 267 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 268 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 269 PetscFPrintf(snes->comm,fd," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma); 270 PetscFPrintf(snes->comm,fd," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n", 271 tr->delta0,tr->delta1,tr->delta2,tr->delta3); 272 } 273 return 0; 274 } 275 276 /* ---------------------------------------------------------------- */ 277 #undef __FUNC__ 278 #define __FUNC__ "SNESConverged_EQ_TR" 279 /*@ 280 SNESConverged_EQ_TR - Monitors the convergence of the trust region 281 method SNES_EQ_TR for solving systems of nonlinear equations (default). 282 283 Input Parameters: 284 . snes - the SNES context 285 . xnorm - 2-norm of current iterate 286 . pnorm - 2-norm of current step 287 . fnorm - 2-norm of function 288 . dummy - unused context 289 290 Returns: 291 $ 1 if ( delta < xnorm*deltatol ), 292 $ 2 if ( fnorm < atol ), 293 $ 3 if ( pnorm < xtol*xnorm ), 294 $ -2 if ( nfct > maxf ), 295 $ -1 if ( delta < xnorm*epsmch ), 296 $ 0 otherwise, 297 298 where 299 $ delta - trust region paramenter 300 $ deltatol - trust region size tolerance, 301 $ set with SNESSetTrustRegionTolerance() 302 $ maxf - maximum number of function evaluations, 303 $ set with SNESSetTolerances() 304 $ nfct - number of function evaluations, 305 $ atol - absolute function norm tolerance, 306 $ set with SNESSetTolerances() 307 $ xtol - relative function norm tolerance, 308 $ set with SNESSetTolerances() 309 310 .keywords: SNES, nonlinear, default, converged, convergence 311 312 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 313 @*/ 314 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,void *dummy) 315 { 316 SNES_TR *neP = (SNES_TR *)snes->data; 317 double epsmch = 1.0e-14; /* This must be fixed */ 318 int info; 319 320 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 321 SETERRQ(1,0,"For SNES_NONLINEAR_EQUATIONS only"); 322 323 if (fnorm != fnorm) { 324 PLogInfo(snes,"SNES:Failed to converged, function norm is NaN\n"); 325 return -3; 326 } 327 if (neP->delta < xnorm * snes->deltatol) { 328 PLogInfo(snes, 329 "SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 330 return 1; 331 } 332 if (neP->itflag) { 333 info = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,dummy); 334 if (info) return info; 335 } else if (snes->nfuncs > snes->max_funcs) { 336 PLogInfo(snes, 337 "SNES: Exceeded maximum number of function evaluations: %d > %d\n", 338 snes->nfuncs, snes->max_funcs ); 339 return -2; 340 } 341 if (neP->delta < xnorm * epsmch) { 342 PLogInfo(snes, 343 "SNESConverged_EQ_TR: Converged due to trust region param %g < %g * %g\n",neP->delta,xnorm, epsmch); 344 return -1; 345 } 346 return 0; 347 } 348 /* ------------------------------------------------------------ */ 349 #undef __FUNC__ 350 #define __FUNC__ "SNESCreate_EQ_TR" 351 int SNESCreate_EQ_TR(SNES snes ) 352 { 353 SNES_TR *neP; 354 355 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 356 SETERRQ(1,0,"For SNES_NONLINEAR_EQUATIONS only"); 357 snes->type = SNES_EQ_TR; 358 snes->setup = SNESSetUp_EQ_TR; 359 snes->solve = SNESSolve_EQ_TR; 360 snes->destroy = SNESDestroy_EQ_TR; 361 snes->converged = SNESConverged_EQ_TR; 362 snes->printhelp = SNESPrintHelp_EQ_TR; 363 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 364 snes->view = SNESView_EQ_TR; 365 snes->nwork = 0; 366 367 neP = PetscNew(SNES_TR); CHKPTRQ(neP); 368 PLogObjectMemory(snes,sizeof(SNES_TR)); 369 snes->data = (void *) neP; 370 neP->mu = 0.25; 371 neP->eta = 0.75; 372 neP->delta = 0.0; 373 neP->delta0 = 0.2; 374 neP->delta1 = 0.3; 375 neP->delta2 = 0.75; 376 neP->delta3 = 2.0; 377 neP->sigma = 0.0001; 378 neP->itflag = 0; 379 neP->rnorm0 = 0; 380 neP->ttol = 0; 381 return 0; 382 } 383