1 #ifndef lint 2 static char vcid[] = "$Id: tr.c,v 1.67 1997/01/06 20:29:57 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; 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_len; /* 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 || */ snes->iter = 0; 81 82 ierr = SNESComputeFunction(snes,X,F); CHKERRQ(ierr); /* F(X) */ 83 ierr = VecNorm(F, NORM_2,&fnorm ); CHKERRQ(ierr); /* fnorm <- || F || */ 84 snes->norm = fnorm; 85 if (history && history_len > 0) history[0] = fnorm; 86 delta = neP->delta0*fnorm; 87 neP->delta = delta; 88 SNESMonitor(snes,0,fnorm); 89 90 /* set parameter for default relative tolerance convergence test */ 91 snes->ttol = fnorm*snes->rtol; 92 93 /* Set the stopping criteria to use the More' trick. */ 94 ierr = SNESGetSLES(snes,&sles); CHKERRQ(ierr); 95 ierr = SLESGetKSP(sles,&ksp); CHKERRQ(ierr); 96 ierr = KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,(void *)snes);CHKERRQ(ierr); 97 98 for ( i=0; i<maxits; i++ ) { 99 snes->iter = i+1; 100 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 101 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 102 103 /* Solve J Y = F, where J is Jacobian matrix */ 104 ierr = SLESSolve(snes->sles,F,Ytmp,&lits); CHKERRQ(ierr); 105 snes->linear_its += PetscAbsInt(lits); 106 PLogInfo(snes,"SNESSolve_EQ_TR: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 107 ierr = VecNorm(Ytmp,NORM_2,&norm); CHKERRQ(ierr); 108 norm1 = norm; 109 while(1) { 110 ierr = VecCopy(Ytmp,Y); CHKERRQ(ierr); 111 norm = norm1; 112 113 /* Scale Y if need be and predict new value of F norm */ 114 if (norm >= delta) { 115 norm = delta/norm; 116 gpnorm = (1.0 - norm)*fnorm; 117 cnorm = norm; 118 PLogInfo(snes,"SNESSolve_EQ_TR: Scaling direction by %g\n",norm ); 119 ierr = VecScale(&cnorm,Y); CHKERRQ(ierr); 120 norm = gpnorm; 121 ynorm = delta; 122 } else { 123 gpnorm = 0.0; 124 PLogInfo(snes,"SNESSolve_EQ_TR: Direction is in Trust Region\n" ); 125 ynorm = norm; 126 } 127 ierr = VecAYPX(&mone,X,Y); CHKERRQ(ierr); /* Y <- X - Y */ 128 ierr = VecCopy(X,snes->vec_sol_update_always); CHKERRQ(ierr); 129 ierr = SNESComputeFunction(snes,Y,G); CHKERRQ(ierr); /* F(X) */ 130 ierr = VecNorm(G,NORM_2,&gnorm); CHKERRQ(ierr); /* gnorm <- || g || */ 131 if (fnorm == gpnorm) rho = 0.0; 132 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 133 134 /* Update size of trust region */ 135 if (rho < neP->mu) delta *= neP->delta1; 136 else if (rho < neP->eta) delta *= neP->delta2; 137 else delta *= neP->delta3; 138 PLogInfo(snes,"SNESSolve_EQ_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 139 PLogInfo(snes,"SNESSolve_EQ_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta); 140 neP->delta = delta; 141 if (rho > neP->sigma) break; 142 PLogInfo(snes,"SNESSolve_EQ_TR: Trying again in smaller region\n"); 143 /* check to see if progress is hopeless */ 144 neP->itflag = 0; 145 if ((*snes->converged)(snes,xnorm,ynorm,fnorm,snes->cnvP)) { 146 /* We're not progressing, so return with the current iterate */ 147 if (X != snes->vec_sol) { 148 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 149 snes->vec_sol_always = snes->vec_sol; 150 snes->vec_func_always = snes->vec_func; 151 } 152 } 153 snes->nfailures++; 154 } 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 /* Verify solution is in corect location */ 167 if (X != snes->vec_sol) { 168 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 169 snes->vec_sol_always = snes->vec_sol; 170 snes->vec_func_always = snes->vec_func; 171 } 172 break; 173 } 174 } 175 if (i == maxits) { 176 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 177 i--; 178 } 179 *its = i+1; 180 return 0; 181 } 182 /*------------------------------------------------------------*/ 183 #undef __FUNC__ 184 #define __FUNC__ "SNESSetUp_EQ_TR" 185 static int SNESSetUp_EQ_TR( SNES snes ) 186 { 187 int ierr; 188 snes->nwork = 4; 189 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERRQ(ierr); 190 PLogObjectParents(snes,snes->nwork,snes->work); 191 snes->vec_sol_update_always = snes->work[3]; 192 return 0; 193 } 194 /*------------------------------------------------------------*/ 195 #undef __FUNC__ 196 #define __FUNC__ "SNESDestroy_EQ_TR" 197 static int SNESDestroy_EQ_TR(PetscObject obj ) 198 { 199 SNES snes = (SNES) obj; 200 int ierr; 201 202 if (snes->nwork) { 203 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 204 } 205 PetscFree(snes->data); 206 return 0; 207 } 208 /*------------------------------------------------------------*/ 209 210 #undef __FUNC__ 211 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 212 static int SNESSetFromOptions_EQ_TR(SNES snes) 213 { 214 SNES_TR *ctx = (SNES_TR *)snes->data; 215 double tmp; 216 int ierr,flg; 217 218 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp, &flg); CHKERRQ(ierr); 219 if (flg) {ctx->mu = tmp;} 220 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp, &flg); CHKERRQ(ierr); 221 if (flg) {ctx->eta = tmp;} 222 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp, &flg); CHKERRQ(ierr); 223 if (flg) {ctx->sigma = tmp;} 224 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp, &flg); CHKERRQ(ierr); 225 if (flg) {ctx->delta0 = tmp;} 226 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp, &flg); CHKERRQ(ierr); 227 if (flg) {ctx->delta1 = tmp;} 228 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp, &flg); CHKERRQ(ierr); 229 if (flg) {ctx->delta2 = tmp;} 230 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp, &flg); CHKERRQ(ierr); 231 if (flg) {ctx->delta3 = tmp;} 232 return 0; 233 } 234 235 #undef __FUNC__ 236 #define __FUNC__ "SNESPrintHelp_EQ_TR" 237 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 238 { 239 SNES_TR *ctx = (SNES_TR *)snes->data; 240 241 PetscFPrintf(snes->comm,stdout," method SNES_EQ_TR (tr) for systems of nonlinear equations:\n"); 242 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu); 243 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta); 244 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma); 245 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0); 246 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1); 247 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2); 248 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3); 249 return 0; 250 } 251 252 #undef __FUNC__ 253 #define __FUNC__ "SNESView_EQ_TR" 254 static int SNESView_EQ_TR(PetscObject obj,Viewer viewer) 255 { 256 SNES snes = (SNES)obj; 257 SNES_TR *tr = (SNES_TR *)snes->data; 258 FILE *fd; 259 int ierr; 260 ViewerType vtype; 261 262 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 263 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 264 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 265 PetscFPrintf(snes->comm,fd," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma); 266 PetscFPrintf(snes->comm,fd," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n", 267 tr->delta0,tr->delta1,tr->delta2,tr->delta3); 268 } 269 return 0; 270 } 271 272 /* ---------------------------------------------------------------- */ 273 #undef __FUNC__ 274 #define __FUNC__ "SNESConverged_EQ_TR" 275 /*@ 276 SNESConverged_EQ_TR - Monitors the convergence of the trust region 277 method SNES_EQ_TR for solving systems of nonlinear equations (default). 278 279 Input Parameters: 280 . snes - the SNES context 281 . xnorm - 2-norm of current iterate 282 . pnorm - 2-norm of current step 283 . fnorm - 2-norm of function 284 . dummy - unused context 285 286 Returns: 287 $ 1 if ( delta < xnorm*deltatol ), 288 $ 2 if ( fnorm < atol ), 289 $ 3 if ( pnorm < xtol*xnorm ), 290 $ -2 if ( nfct > maxf ), 291 $ -1 if ( delta < xnorm*epsmch ), 292 $ 0 otherwise, 293 294 where 295 $ delta - trust region paramenter 296 $ deltatol - trust region size tolerance, 297 $ set with SNESSetTrustRegionTolerance() 298 $ maxf - maximum number of function evaluations, 299 $ set with SNESSetTolerances() 300 $ nfct - number of function evaluations, 301 $ atol - absolute function norm tolerance, 302 $ set with SNESSetTolerances() 303 $ xtol - relative function norm tolerance, 304 $ set with SNESSetTolerances() 305 306 .keywords: SNES, nonlinear, default, converged, convergence 307 308 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 309 @*/ 310 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,void *dummy) 311 { 312 SNES_TR *neP = (SNES_TR *)snes->data; 313 double epsmch = 1.0e-14; /* This must be fixed */ 314 int info; 315 316 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 317 SETERRQ(1,0,"For SNES_NONLINEAR_EQUATIONS only"); 318 319 if (neP->delta < xnorm * snes->deltatol) { 320 PLogInfo(snes, 321 "SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 322 return 1; 323 } 324 if (neP->itflag) { 325 info = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,dummy); 326 if (info) return info; 327 } 328 if (neP->delta < xnorm * epsmch) { 329 PLogInfo(snes, 330 "SNESConverged_EQ_TR: Converged due to trust region param %g < %g * %g\n",neP->delta,xnorm, epsmch); 331 return -1; 332 } 333 return 0; 334 } 335 /* ------------------------------------------------------------ */ 336 #undef __FUNC__ 337 #define __FUNC__ "SNESCreate_EQ_TR" 338 int SNESCreate_EQ_TR(SNES snes ) 339 { 340 SNES_TR *neP; 341 342 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 343 SETERRQ(1,0,"For SNES_NONLINEAR_EQUATIONS only"); 344 snes->type = SNES_EQ_TR; 345 snes->setup = SNESSetUp_EQ_TR; 346 snes->solve = SNESSolve_EQ_TR; 347 snes->destroy = SNESDestroy_EQ_TR; 348 snes->converged = SNESConverged_EQ_TR; 349 snes->printhelp = SNESPrintHelp_EQ_TR; 350 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 351 snes->view = SNESView_EQ_TR; 352 snes->nwork = 0; 353 354 neP = PetscNew(SNES_TR); CHKPTRQ(neP); 355 PLogObjectMemory(snes,sizeof(SNES_TR)); 356 snes->data = (void *) neP; 357 neP->mu = 0.25; 358 neP->eta = 0.75; 359 neP->delta = 0.0; 360 neP->delta0 = 0.2; 361 neP->delta1 = 0.3; 362 neP->delta2 = 0.75; 363 neP->delta3 = 2.0; 364 neP->sigma = 0.0001; 365 neP->itflag = 0; 366 neP->rnorm0 = 0; 367 neP->ttol = 0; 368 return 0; 369 } 370