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