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