1 #ifndef lint 2 static char vcid[] = "$Id: tr.c,v 1.68 1997/01/14 22:58:25 curfman Exp bsmith $"; 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 } 148 snes->nfailures++; 149 } 150 fnorm = gnorm; 151 snes->norm = fnorm; 152 if (history && history_len > i+1) history[i+1] = fnorm; 153 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 154 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 155 VecNorm(X, NORM_2,&xnorm ); /* xnorm = || X || */ 156 SNESMonitor(snes,i+1,fnorm); 157 158 /* Test for convergence */ 159 neP->itflag = 1; 160 if ((*snes->converged)( snes, xnorm, ynorm, fnorm,snes->cnvP )) { 161 break; 162 } 163 } 164 if (X != snes->vec_sol) { 165 /* Verify solution is in corect location */ 166 ierr = VecCopy(X,snes->vec_sol); CHKERRQ(ierr); 167 snes->vec_sol_always = snes->vec_sol; 168 snes->vec_func_always = snes->vec_func; 169 } 170 if (i == maxits) { 171 PLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 172 i--; 173 } 174 *its = i+1; 175 return 0; 176 } 177 /*------------------------------------------------------------*/ 178 #undef __FUNC__ 179 #define __FUNC__ "SNESSetUp_EQ_TR" 180 static int SNESSetUp_EQ_TR( SNES snes ) 181 { 182 int ierr; 183 snes->nwork = 4; 184 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work ); CHKERRQ(ierr); 185 PLogObjectParents(snes,snes->nwork,snes->work); 186 snes->vec_sol_update_always = snes->work[3]; 187 return 0; 188 } 189 /*------------------------------------------------------------*/ 190 #undef __FUNC__ 191 #define __FUNC__ "SNESDestroy_EQ_TR" 192 static int SNESDestroy_EQ_TR(PetscObject obj ) 193 { 194 SNES snes = (SNES) obj; 195 int ierr; 196 197 if (snes->nwork) { 198 ierr = VecDestroyVecs(snes->work,snes->nwork); CHKERRQ(ierr); 199 } 200 PetscFree(snes->data); 201 return 0; 202 } 203 /*------------------------------------------------------------*/ 204 205 #undef __FUNC__ 206 #define __FUNC__ "SNESSetFromOptions_EQ_TR" 207 static int SNESSetFromOptions_EQ_TR(SNES snes) 208 { 209 SNES_TR *ctx = (SNES_TR *)snes->data; 210 double tmp; 211 int ierr,flg; 212 213 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_mu",&tmp, &flg); CHKERRQ(ierr); 214 if (flg) {ctx->mu = tmp;} 215 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_eta",&tmp, &flg); CHKERRQ(ierr); 216 if (flg) {ctx->eta = tmp;} 217 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_sigma",&tmp, &flg); CHKERRQ(ierr); 218 if (flg) {ctx->sigma = tmp;} 219 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta0",&tmp, &flg); CHKERRQ(ierr); 220 if (flg) {ctx->delta0 = tmp;} 221 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta1",&tmp, &flg); CHKERRQ(ierr); 222 if (flg) {ctx->delta1 = tmp;} 223 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta2",&tmp, &flg); CHKERRQ(ierr); 224 if (flg) {ctx->delta2 = tmp;} 225 ierr = OptionsGetDouble(snes->prefix,"-snes_eq_tr_delta3",&tmp, &flg); CHKERRQ(ierr); 226 if (flg) {ctx->delta3 = tmp;} 227 return 0; 228 } 229 230 #undef __FUNC__ 231 #define __FUNC__ "SNESPrintHelp_EQ_TR" 232 static int SNESPrintHelp_EQ_TR(SNES snes,char *p) 233 { 234 SNES_TR *ctx = (SNES_TR *)snes->data; 235 236 PetscFPrintf(snes->comm,stdout," method SNES_EQ_TR (tr) for systems of nonlinear equations:\n"); 237 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_mu <mu> (default %g)\n",p,ctx->mu); 238 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_eta <eta> (default %g)\n",p,ctx->eta); 239 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_sigma <sigma> (default %g)\n",p,ctx->sigma); 240 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta0 <delta0> (default %g)\n",p,ctx->delta0); 241 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta1 <delta1> (default %g)\n",p,ctx->delta1); 242 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta2 <delta2> (default %g)\n",p,ctx->delta2); 243 PetscFPrintf(snes->comm,stdout," %ssnes_eq_tr_delta3 <delta3> (default %g)\n",p,ctx->delta3); 244 return 0; 245 } 246 247 #undef __FUNC__ 248 #define __FUNC__ "SNESView_EQ_TR" 249 static int SNESView_EQ_TR(PetscObject obj,Viewer viewer) 250 { 251 SNES snes = (SNES)obj; 252 SNES_TR *tr = (SNES_TR *)snes->data; 253 FILE *fd; 254 int ierr; 255 ViewerType vtype; 256 257 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 258 if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 259 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 260 PetscFPrintf(snes->comm,fd," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma); 261 PetscFPrintf(snes->comm,fd," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n", 262 tr->delta0,tr->delta1,tr->delta2,tr->delta3); 263 } 264 return 0; 265 } 266 267 /* ---------------------------------------------------------------- */ 268 #undef __FUNC__ 269 #define __FUNC__ "SNESConverged_EQ_TR" 270 /*@ 271 SNESConverged_EQ_TR - Monitors the convergence of the trust region 272 method SNES_EQ_TR for solving systems of nonlinear equations (default). 273 274 Input Parameters: 275 . snes - the SNES context 276 . xnorm - 2-norm of current iterate 277 . pnorm - 2-norm of current step 278 . fnorm - 2-norm of function 279 . dummy - unused context 280 281 Returns: 282 $ 1 if ( delta < xnorm*deltatol ), 283 $ 2 if ( fnorm < atol ), 284 $ 3 if ( pnorm < xtol*xnorm ), 285 $ -2 if ( nfct > maxf ), 286 $ -1 if ( delta < xnorm*epsmch ), 287 $ 0 otherwise, 288 289 where 290 $ delta - trust region paramenter 291 $ deltatol - trust region size tolerance, 292 $ set with SNESSetTrustRegionTolerance() 293 $ maxf - maximum number of function evaluations, 294 $ set with SNESSetTolerances() 295 $ nfct - number of function evaluations, 296 $ atol - absolute function norm tolerance, 297 $ set with SNESSetTolerances() 298 $ xtol - relative function norm tolerance, 299 $ set with SNESSetTolerances() 300 301 .keywords: SNES, nonlinear, default, converged, convergence 302 303 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 304 @*/ 305 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,void *dummy) 306 { 307 SNES_TR *neP = (SNES_TR *)snes->data; 308 double epsmch = 1.0e-14; /* This must be fixed */ 309 int info; 310 311 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 312 SETERRQ(1,0,"For SNES_NONLINEAR_EQUATIONS only"); 313 314 if (neP->delta < xnorm * snes->deltatol) { 315 PLogInfo(snes, 316 "SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 317 return 1; 318 } 319 if (neP->itflag) { 320 info = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,dummy); 321 if (info) return info; 322 } 323 if (neP->delta < xnorm * epsmch) { 324 PLogInfo(snes, 325 "SNESConverged_EQ_TR: Converged due to trust region param %g < %g * %g\n",neP->delta,xnorm, epsmch); 326 return -1; 327 } 328 return 0; 329 } 330 /* ------------------------------------------------------------ */ 331 #undef __FUNC__ 332 #define __FUNC__ "SNESCreate_EQ_TR" 333 int SNESCreate_EQ_TR(SNES snes ) 334 { 335 SNES_TR *neP; 336 337 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) 338 SETERRQ(1,0,"For SNES_NONLINEAR_EQUATIONS only"); 339 snes->type = SNES_EQ_TR; 340 snes->setup = SNESSetUp_EQ_TR; 341 snes->solve = SNESSolve_EQ_TR; 342 snes->destroy = SNESDestroy_EQ_TR; 343 snes->converged = SNESConverged_EQ_TR; 344 snes->printhelp = SNESPrintHelp_EQ_TR; 345 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 346 snes->view = SNESView_EQ_TR; 347 snes->nwork = 0; 348 349 neP = PetscNew(SNES_TR); CHKPTRQ(neP); 350 PLogObjectMemory(snes,sizeof(SNES_TR)); 351 snes->data = (void *) neP; 352 neP->mu = 0.25; 353 neP->eta = 0.75; 354 neP->delta = 0.0; 355 neP->delta0 = 0.2; 356 neP->delta1 = 0.3; 357 neP->delta2 = 0.75; 358 neP->delta3 = 2.0; 359 neP->sigma = 0.0001; 360 neP->itflag = 0; 361 neP->rnorm0 = 0; 362 neP->ttol = 0; 363 return 0; 364 } 365