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