1 /*$Id: tr.c,v 1.123 2001/03/23 23:24:15 balay Exp curfman $*/ 2 3 #include "src/snes/impls/tr/tr.h" /*I "petscsnes.h" I*/ 4 5 /* 6 This convergence test determines if the two norm of the 7 solution lies outside the trust region, if so it halts. 8 */ 9 #undef __FUNCT__ 10 #define __FUNCT__ "SNES_EQ_TR_KSPConverged_Private" 11 int SNES_EQ_TR_KSPConverged_Private(KSP ksp,int n,double rnorm,KSPConvergedReason *reason,void *ctx) 12 { 13 SNES snes = (SNES) ctx; 14 SNES_KSP_EW_ConvCtx *kctx = (SNES_KSP_EW_ConvCtx*)snes->kspconvctx; 15 SNES_EQ_TR *neP = (SNES_EQ_TR*)snes->data; 16 Vec x; 17 double norm; 18 int ierr; 19 20 PetscFunctionBegin; 21 if (snes->ksp_ewconv) { 22 if (!kctx) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Eisenstat-Walker onvergence context not created"); 23 if (!n) {ierr = SNES_KSP_EW_ComputeRelativeTolerance_Private(snes,ksp);CHKERRQ(ierr);} 24 kctx->lresid_last = rnorm; 25 } 26 ierr = KSPDefaultConverged(ksp,n,rnorm,reason,ctx);CHKERRQ(ierr); 27 if (*reason) { 28 PetscLogInfo(snes,"SNES_EQ_TR_KSPConverged_Private: regular convergence test KSP iterations=%d, rnorm=%g\n",n,rnorm); 29 } 30 31 /* Determine norm of solution */ 32 ierr = KSPBuildSolution(ksp,0,&x);CHKERRQ(ierr); 33 ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); 34 if (norm >= neP->delta) { 35 PetscLogInfo(snes,"SNES_EQ_TR_KSPConverged_Private: KSP iterations=%d, rnorm=%g\n",n,rnorm); 36 PetscLogInfo(snes,"SNES_EQ_TR_KSPConverged_Private: Ending linear iteration early, delta=%g, length=%g\n",neP->delta,norm); 37 *reason = KSP_CONVERGED_STEP_LENGTH; 38 } 39 PetscFunctionReturn(0); 40 } 41 42 /* 43 SNESSolve_EQ_TR - Implements Newton's Method with a very simple trust 44 region approach for solving systems of nonlinear equations. 45 46 The basic algorithm is taken from "The Minpack Project", by More', 47 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 48 of Mathematical Software", Wayne Cowell, editor. 49 50 This is intended as a model implementation, since it does not 51 necessarily have many of the bells and whistles of other 52 implementations. 53 */ 54 #undef __FUNCT__ 55 #define __FUNCT__ "SNESSolve_EQ_TR" 56 static int SNESSolve_EQ_TR(SNES snes,int *its) 57 { 58 SNES_EQ_TR *neP = (SNES_EQ_TR*)snes->data; 59 Vec X,F,Y,G,TMP,Ytmp; 60 int maxits,i,ierr,lits,breakout = 0; 61 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 62 double rho,fnorm,gnorm,gpnorm,xnorm,delta,norm,ynorm,norm1; 63 Scalar mone = -1.0,cnorm; 64 KSP ksp; 65 SLES sles; 66 SNESConvergedReason reason; 67 PetscTruth conv; 68 69 PetscFunctionBegin; 70 maxits = snes->max_its; /* maximum number of iterations */ 71 X = snes->vec_sol; /* solution vector */ 72 F = snes->vec_func; /* residual vector */ 73 Y = snes->work[0]; /* work vectors */ 74 G = snes->work[1]; 75 Ytmp = snes->work[2]; 76 77 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 78 snes->iter = 0; 79 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 80 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm = || X || */ 81 82 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); /* F(X) */ 83 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- || F || */ 84 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 85 snes->norm = fnorm; 86 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 87 delta = neP->delta0*fnorm; 88 neP->delta = delta; 89 SNESLogConvHistory(snes,fnorm,0); 90 SNESMonitor(snes,0,fnorm); 91 92 if (fnorm < snes->atol) {*its = 0; snes->reason = SNES_CONVERGED_FNORM_ABS; PetscFunctionReturn(0);} 93 94 /* set parameter for default relative tolerance convergence test */ 95 snes->ttol = fnorm*snes->rtol; 96 97 /* Set the stopping criteria to use the More' trick. */ 98 ierr = PetscOptionsHasName(PETSC_NULL,"-snes_tr_ksp_regular_convergence_test",&conv);CHKERRQ(ierr); 99 if (!conv) { 100 ierr = SNESGetSLES(snes,&sles);CHKERRQ(ierr); 101 ierr = SLESGetKSP(sles,&ksp);CHKERRQ(ierr); 102 ierr = KSPSetConvergenceTest(ksp,SNES_EQ_TR_KSPConverged_Private,(void *)snes);CHKERRQ(ierr); 103 } 104 105 for (i=0; i<maxits; i++) { 106 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 107 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 108 109 /* Solve J Y = F, where J is Jacobian matrix */ 110 ierr = SLESSolve(snes->sles,F,Ytmp,&lits);CHKERRQ(ierr); 111 snes->linear_its += lits; 112 PetscLogInfo(snes,"SNESSolve_EQ_TR: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 113 ierr = VecNorm(Ytmp,NORM_2,&norm);CHKERRQ(ierr); 114 norm1 = norm; 115 while(1) { 116 ierr = VecCopy(Ytmp,Y);CHKERRQ(ierr); 117 norm = norm1; 118 119 /* Scale Y if need be and predict new value of F norm */ 120 if (norm >= delta) { 121 norm = delta/norm; 122 gpnorm = (1.0 - norm)*fnorm; 123 cnorm = norm; 124 PetscLogInfo(snes,"SNESSolve_EQ_TR: Scaling direction by %g\n",norm); 125 ierr = VecScale(&cnorm,Y);CHKERRQ(ierr); 126 norm = gpnorm; 127 ynorm = delta; 128 } else { 129 gpnorm = 0.0; 130 PetscLogInfo(snes,"SNESSolve_EQ_TR: Direction is in Trust Region\n"); 131 ynorm = norm; 132 } 133 ierr = VecAYPX(&mone,X,Y);CHKERRQ(ierr); /* Y <- X - Y */ 134 ierr = VecCopy(X,snes->vec_sol_update_always);CHKERRQ(ierr); 135 ierr = SNESComputeFunction(snes,Y,G);CHKERRQ(ierr); /* F(X) */ 136 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); /* gnorm <- || g || */ 137 if (fnorm == gpnorm) rho = 0.0; 138 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 139 140 /* Update size of trust region */ 141 if (rho < neP->mu) delta *= neP->delta1; 142 else if (rho < neP->eta) delta *= neP->delta2; 143 else delta *= neP->delta3; 144 PetscLogInfo(snes,"SNESSolve_EQ_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 145 PetscLogInfo(snes,"SNESSolve_EQ_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta); 146 neP->delta = delta; 147 if (rho > neP->sigma) break; 148 PetscLogInfo(snes,"SNESSolve_EQ_TR: Trying again in smaller region\n"); 149 /* check to see if progress is hopeless */ 150 neP->itflag = 0; 151 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 152 if (reason) { 153 /* We're not progressing, so return with the current iterate */ 154 breakout = 1; 155 break; 156 } 157 snes->nfailures++; 158 } 159 if (!breakout) { 160 fnorm = gnorm; 161 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 162 snes->iter = i+1; 163 snes->norm = fnorm; 164 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 165 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 166 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 167 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm = || X || */ 168 SNESLogConvHistory(snes,fnorm,lits); 169 SNESMonitor(snes,i+1,fnorm); 170 171 /* Test for convergence */ 172 neP->itflag = 1; 173 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 174 if (reason) { 175 break; 176 } 177 } else { 178 break; 179 } 180 } 181 /* Verify solution is in corect location */ 182 if (X != snes->vec_sol) { 183 ierr = VecCopy(X,snes->vec_sol);CHKERRQ(ierr); 184 } 185 if (F != snes->vec_func) { 186 ierr = VecCopy(F,snes->vec_func);CHKERRQ(ierr); 187 } 188 snes->vec_sol_always = snes->vec_sol; 189 snes->vec_func_always = snes->vec_func; 190 if (i == maxits) { 191 PetscLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 192 i--; 193 reason = SNES_DIVERGED_MAX_IT; 194 } 195 *its = i+1; 196 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 197 snes->reason = reason; 198 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 199 PetscFunctionReturn(0); 200 } 201 /*------------------------------------------------------------*/ 202 #undef __FUNCT__ 203 #define __FUNCT__ "SNESSetUp_EQ_TR" 204 static int SNESSetUp_EQ_TR(SNES snes) 205 { 206 int ierr; 207 208 PetscFunctionBegin; 209 snes->nwork = 4; 210 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 211 PetscLogObjectParents(snes,snes->nwork,snes->work); 212 snes->vec_sol_update_always = snes->work[3]; 213 PetscFunctionReturn(0); 214 } 215 /*------------------------------------------------------------*/ 216 #undef __FUNCT__ 217 #define __FUNCT__ "SNESDestroy_EQ_TR" 218 static int SNESDestroy_EQ_TR(SNES snes) 219 { 220 int ierr; 221 222 PetscFunctionBegin; 223 if (snes->nwork) { 224 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 225 } 226 ierr = PetscFree(snes->data);CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 /*------------------------------------------------------------*/ 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "SNESSetFromOptions_EQ_TR" 233 static int SNESSetFromOptions_EQ_TR(SNES snes) 234 { 235 SNES_EQ_TR *ctx = (SNES_EQ_TR *)snes->data; 236 int ierr; 237 238 PetscFunctionBegin; 239 ierr = PetscOptionsHead("SNES trust region options for nonlinear equations");CHKERRQ(ierr); 240 ierr = PetscOptionsDouble("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,0);CHKERRQ(ierr); 241 ierr = PetscOptionsDouble("-snes_eq_tr_mu","mu","None",ctx->mu,&ctx->mu,0);CHKERRQ(ierr); 242 ierr = PetscOptionsDouble("-snes_eq_tr_eta","eta","None",ctx->eta,&ctx->eta,0);CHKERRQ(ierr); 243 ierr = PetscOptionsDouble("-snes_eq_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,0);CHKERRQ(ierr); 244 ierr = PetscOptionsDouble("-snes_eq_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,0);CHKERRQ(ierr); 245 ierr = PetscOptionsDouble("-snes_eq_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,0);CHKERRQ(ierr); 246 ierr = PetscOptionsDouble("-snes_eq_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,0);CHKERRQ(ierr); 247 ierr = PetscOptionsDouble("-snes_eq_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,0);CHKERRQ(ierr); 248 ierr = PetscOptionsTail();CHKERRQ(ierr); 249 PetscFunctionReturn(0); 250 } 251 252 #undef __FUNCT__ 253 #define __FUNCT__ "SNESView_EQ_TR" 254 static int SNESView_EQ_TR(SNES snes,PetscViewer viewer) 255 { 256 SNES_EQ_TR *tr = (SNES_EQ_TR *)snes->data; 257 int ierr; 258 PetscTruth isascii; 259 260 PetscFunctionBegin; 261 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 262 if (isascii) { 263 ierr = PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 264 ierr = PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 265 } else { 266 SETERRQ1(1,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name); 267 } 268 PetscFunctionReturn(0); 269 } 270 271 /* ---------------------------------------------------------------- */ 272 #undef __FUNCT__ 273 #define __FUNCT__ "SNESConverged_EQ_TR" 274 /*@C 275 SNESConverged_EQ_TR - Monitors the convergence of the trust region 276 method SNESEQTR for solving systems of nonlinear equations (default). 277 278 Collective on SNES 279 280 Input Parameters: 281 + snes - the SNES context 282 . xnorm - 2-norm of current iterate 283 . pnorm - 2-norm of current step 284 . fnorm - 2-norm of function 285 - dummy - unused context 286 287 Output Parameter: 288 . reason - one of 289 $ SNES_CONVERGED_FNORM_ABS - (fnorm < atol), 290 $ SNES_CONVERGED_PNORM_RELATIVE - (pnorm < xtol*xnorm), 291 $ SNES_CONVERGED_FNORM_RELATIVE - (fnorm < rtol*fnorm0), 292 $ SNES_DIVERGED_FUNCTION_COUNT - (nfct > maxf), 293 $ SNES_DIVERGED_FNORM_NAN - (fnorm == NaN), 294 $ SNES_CONVERGED_TR_DELTA - (delta < xnorm*deltatol), 295 $ SNES_CONVERGED_ITERATING - (otherwise) 296 297 where 298 + delta - trust region paramenter 299 . deltatol - trust region size tolerance, 300 set with SNESSetTrustRegionTolerance() 301 . maxf - maximum number of function evaluations, 302 set with SNESSetTolerances() 303 . nfct - number of function evaluations, 304 . atol - absolute function norm tolerance, 305 set with SNESSetTolerances() 306 - xtol - relative function norm tolerance, 307 set with SNESSetTolerances() 308 309 Level: intermediate 310 311 .keywords: SNES, nonlinear, default, converged, convergence 312 313 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 314 @*/ 315 int SNESConverged_EQ_TR(SNES snes,double xnorm,double pnorm,double fnorm,SNESConvergedReason *reason,void *dummy) 316 { 317 SNES_EQ_TR *neP = (SNES_EQ_TR *)snes->data; 318 int ierr; 319 320 PetscFunctionBegin; 321 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 322 SETERRQ(PETSC_ERR_ARG_WRONG,"For SNES_NONLINEAR_EQUATIONS only"); 323 } 324 325 if (fnorm != fnorm) { 326 PetscLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 327 *reason = SNES_DIVERGED_FNORM_NAN; 328 } else if (neP->delta < xnorm * snes->deltatol) { 329 PetscLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 330 *reason = SNES_CONVERGED_TR_DELTA; 331 } else if (neP->itflag) { 332 ierr = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,reason,dummy);CHKERRQ(ierr); 333 } else if (snes->nfuncs > snes->max_funcs) { 334 PetscLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n",snes->nfuncs,snes->max_funcs); 335 *reason = SNES_DIVERGED_FUNCTION_COUNT; 336 } else { 337 *reason = SNES_CONVERGED_ITERATING; 338 } 339 PetscFunctionReturn(0); 340 } 341 /* ------------------------------------------------------------ */ 342 EXTERN_C_BEGIN 343 #undef __FUNCT__ 344 #define __FUNCT__ "SNESCreate_EQ_TR" 345 int SNESCreate_EQ_TR(SNES snes) 346 { 347 SNES_EQ_TR *neP; 348 int ierr; 349 350 PetscFunctionBegin; 351 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 352 SETERRQ(PETSC_ERR_ARG_WRONG,"For SNES_NONLINEAR_EQUATIONS only"); 353 } 354 snes->setup = SNESSetUp_EQ_TR; 355 snes->solve = SNESSolve_EQ_TR; 356 snes->destroy = SNESDestroy_EQ_TR; 357 snes->converged = SNESConverged_EQ_TR; 358 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 359 snes->view = SNESView_EQ_TR; 360 snes->nwork = 0; 361 362 ierr = PetscNew(SNES_EQ_TR,&neP);CHKERRQ(ierr); 363 PetscLogObjectMemory(snes,sizeof(SNES_EQ_TR)); 364 snes->data = (void*)neP; 365 neP->mu = 0.25; 366 neP->eta = 0.75; 367 neP->delta = 0.0; 368 neP->delta0 = 0.2; 369 neP->delta1 = 0.3; 370 neP->delta2 = 0.75; 371 neP->delta3 = 2.0; 372 neP->sigma = 0.0001; 373 neP->itflag = 0; 374 neP->rnorm0 = 0; 375 neP->ttol = 0; 376 PetscFunctionReturn(0); 377 } 378 EXTERN_C_END 379 380