1 /*$Id: tr.c,v 1.128 2001/08/07 03:04:11 balay Exp $*/ 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,PetscReal 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 PetscReal 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 PetscReal rho,fnorm,gnorm,gpnorm,xnorm,delta,norm,ynorm,norm1; 63 PetscScalar 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 PetscLogInfo(snes,"SNESSolve_EQ_TR: Using Krylov convergence test SNES_EQ_TR_KSPConverged_Private\n"); 104 } 105 106 for (i=0; i<maxits; i++) { 107 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 108 ierr = SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 109 110 /* Solve J Y = F, where J is Jacobian matrix */ 111 ierr = SLESSolve(snes->sles,F,Ytmp,&lits);CHKERRQ(ierr); 112 snes->linear_its += lits; 113 PetscLogInfo(snes,"SNESSolve_EQ_TR: iter=%d, linear solve iterations=%d\n",snes->iter,lits); 114 ierr = VecNorm(Ytmp,NORM_2,&norm);CHKERRQ(ierr); 115 norm1 = norm; 116 while(1) { 117 ierr = VecCopy(Ytmp,Y);CHKERRQ(ierr); 118 norm = norm1; 119 120 /* Scale Y if need be and predict new value of F norm */ 121 if (norm >= delta) { 122 norm = delta/norm; 123 gpnorm = (1.0 - norm)*fnorm; 124 cnorm = norm; 125 PetscLogInfo(snes,"SNESSolve_EQ_TR: Scaling direction by %g\n",norm); 126 ierr = VecScale(&cnorm,Y);CHKERRQ(ierr); 127 norm = gpnorm; 128 ynorm = delta; 129 } else { 130 gpnorm = 0.0; 131 PetscLogInfo(snes,"SNESSolve_EQ_TR: Direction is in Trust Region\n"); 132 ynorm = norm; 133 } 134 ierr = VecAYPX(&mone,X,Y);CHKERRQ(ierr); /* Y <- X - Y */ 135 ierr = VecCopy(X,snes->vec_sol_update_always);CHKERRQ(ierr); 136 ierr = SNESComputeFunction(snes,Y,G);CHKERRQ(ierr); /* F(X) */ 137 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); /* gnorm <- || g || */ 138 if (fnorm == gpnorm) rho = 0.0; 139 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 140 141 /* Update size of trust region */ 142 if (rho < neP->mu) delta *= neP->delta1; 143 else if (rho < neP->eta) delta *= neP->delta2; 144 else delta *= neP->delta3; 145 PetscLogInfo(snes,"SNESSolve_EQ_TR: fnorm=%g, gnorm=%g, ynorm=%g\n",fnorm,gnorm,ynorm); 146 PetscLogInfo(snes,"SNESSolve_EQ_TR: gpred=%g, rho=%g, delta=%g\n",gpnorm,rho,delta); 147 neP->delta = delta; 148 if (rho > neP->sigma) break; 149 PetscLogInfo(snes,"SNESSolve_EQ_TR: Trying again in smaller region\n"); 150 /* check to see if progress is hopeless */ 151 neP->itflag = 0; 152 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 153 if (reason) { 154 /* We're not progressing, so return with the current iterate */ 155 breakout = 1; 156 break; 157 } 158 snes->nfailures++; 159 } 160 if (!breakout) { 161 fnorm = gnorm; 162 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 163 snes->iter = i+1; 164 snes->norm = fnorm; 165 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 166 TMP = F; F = G; snes->vec_func_always = F; G = TMP; 167 TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP; 168 ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); /* xnorm = || X || */ 169 SNESLogConvHistory(snes,fnorm,lits); 170 SNESMonitor(snes,i+1,fnorm); 171 172 /* Test for convergence */ 173 neP->itflag = 1; 174 ierr = (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 175 if (reason) { 176 break; 177 } 178 } else { 179 break; 180 } 181 } 182 /* Verify solution is in corect location */ 183 if (X != snes->vec_sol) { 184 ierr = VecCopy(X,snes->vec_sol);CHKERRQ(ierr); 185 } 186 if (F != snes->vec_func) { 187 ierr = VecCopy(F,snes->vec_func);CHKERRQ(ierr); 188 } 189 snes->vec_sol_always = snes->vec_sol; 190 snes->vec_func_always = snes->vec_func; 191 if (i == maxits) { 192 PetscLogInfo(snes,"SNESSolve_EQ_TR: Maximum number of iterations has been reached: %d\n",maxits); 193 i--; 194 reason = SNES_DIVERGED_MAX_IT; 195 } 196 *its = i+1; 197 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 198 snes->reason = reason; 199 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 200 PetscFunctionReturn(0); 201 } 202 /*------------------------------------------------------------*/ 203 #undef __FUNCT__ 204 #define __FUNCT__ "SNESSetUp_EQ_TR" 205 static int SNESSetUp_EQ_TR(SNES snes) 206 { 207 int ierr; 208 209 PetscFunctionBegin; 210 snes->nwork = 4; 211 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 212 PetscLogObjectParents(snes,snes->nwork,snes->work); 213 snes->vec_sol_update_always = snes->work[3]; 214 PetscFunctionReturn(0); 215 } 216 /*------------------------------------------------------------*/ 217 #undef __FUNCT__ 218 #define __FUNCT__ "SNESDestroy_EQ_TR" 219 static int SNESDestroy_EQ_TR(SNES snes) 220 { 221 int ierr; 222 223 PetscFunctionBegin; 224 if (snes->nwork) { 225 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 226 } 227 ierr = PetscFree(snes->data);CHKERRQ(ierr); 228 PetscFunctionReturn(0); 229 } 230 /*------------------------------------------------------------*/ 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "SNESSetFromOptions_EQ_TR" 234 static int SNESSetFromOptions_EQ_TR(SNES snes) 235 { 236 SNES_EQ_TR *ctx = (SNES_EQ_TR *)snes->data; 237 int ierr; 238 239 PetscFunctionBegin; 240 ierr = PetscOptionsHead("SNES trust region options for nonlinear equations");CHKERRQ(ierr); 241 ierr = PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,0);CHKERRQ(ierr); 242 ierr = PetscOptionsReal("-snes_eq_tr_mu","mu","None",ctx->mu,&ctx->mu,0);CHKERRQ(ierr); 243 ierr = PetscOptionsReal("-snes_eq_tr_eta","eta","None",ctx->eta,&ctx->eta,0);CHKERRQ(ierr); 244 ierr = PetscOptionsReal("-snes_eq_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,0);CHKERRQ(ierr); 245 ierr = PetscOptionsReal("-snes_eq_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,0);CHKERRQ(ierr); 246 ierr = PetscOptionsReal("-snes_eq_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,0);CHKERRQ(ierr); 247 ierr = PetscOptionsReal("-snes_eq_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,0);CHKERRQ(ierr); 248 ierr = PetscOptionsReal("-snes_eq_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,0);CHKERRQ(ierr); 249 ierr = PetscOptionsTail();CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "SNESView_EQ_TR" 255 static int SNESView_EQ_TR(SNES snes,PetscViewer viewer) 256 { 257 SNES_EQ_TR *tr = (SNES_EQ_TR *)snes->data; 258 int ierr; 259 PetscTruth isascii; 260 261 PetscFunctionBegin; 262 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 263 if (isascii) { 264 ierr = PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 265 ierr = PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 266 } else { 267 SETERRQ1(1,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name); 268 } 269 PetscFunctionReturn(0); 270 } 271 272 /* ---------------------------------------------------------------- */ 273 #undef __FUNCT__ 274 #define __FUNCT__ "SNESConverged_EQ_TR" 275 /*@C 276 SNESConverged_EQ_TR - Monitors the convergence of the trust region 277 method SNESEQTR for solving systems of nonlinear equations (default). 278 279 Collective on SNES 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 Output Parameter: 289 . reason - one of 290 $ SNES_CONVERGED_FNORM_ABS - (fnorm < atol), 291 $ SNES_CONVERGED_PNORM_RELATIVE - (pnorm < xtol*xnorm), 292 $ SNES_CONVERGED_FNORM_RELATIVE - (fnorm < rtol*fnorm0), 293 $ SNES_DIVERGED_FUNCTION_COUNT - (nfct > maxf), 294 $ SNES_DIVERGED_FNORM_NAN - (fnorm == NaN), 295 $ SNES_CONVERGED_TR_DELTA - (delta < xnorm*deltatol), 296 $ SNES_CONVERGED_ITERATING - (otherwise) 297 298 where 299 + delta - trust region paramenter 300 . deltatol - trust region size tolerance, 301 set with SNESSetTrustRegionTolerance() 302 . maxf - maximum number of function evaluations, 303 set with SNESSetTolerances() 304 . nfct - number of function evaluations, 305 . atol - absolute function norm tolerance, 306 set with SNESSetTolerances() 307 - xtol - relative function norm tolerance, 308 set with SNESSetTolerances() 309 310 Level: intermediate 311 312 .keywords: SNES, nonlinear, default, converged, convergence 313 314 .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged() 315 @*/ 316 int SNESConverged_EQ_TR(SNES snes,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) 317 { 318 SNES_EQ_TR *neP = (SNES_EQ_TR *)snes->data; 319 int ierr; 320 321 PetscFunctionBegin; 322 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 323 SETERRQ(PETSC_ERR_ARG_WRONG,"For SNES_NONLINEAR_EQUATIONS only"); 324 } 325 326 if (fnorm != fnorm) { 327 PetscLogInfo(snes,"SNESConverged_EQ_TR:Failed to converged, function norm is NaN\n"); 328 *reason = SNES_DIVERGED_FNORM_NAN; 329 } else if (neP->delta < xnorm * snes->deltatol) { 330 PetscLogInfo(snes,"SNESConverged_EQ_TR: Converged due to trust region param %g<%g*%g\n",neP->delta,xnorm,snes->deltatol); 331 *reason = SNES_CONVERGED_TR_DELTA; 332 } else if (neP->itflag) { 333 ierr = SNESConverged_EQ_LS(snes,xnorm,pnorm,fnorm,reason,dummy);CHKERRQ(ierr); 334 } else if (snes->nfuncs > snes->max_funcs) { 335 PetscLogInfo(snes,"SNESConverged_EQ_TR: Exceeded maximum number of function evaluations: %d > %d\n",snes->nfuncs,snes->max_funcs); 336 *reason = SNES_DIVERGED_FUNCTION_COUNT; 337 } else { 338 *reason = SNES_CONVERGED_ITERATING; 339 } 340 PetscFunctionReturn(0); 341 } 342 /* ------------------------------------------------------------ */ 343 EXTERN_C_BEGIN 344 #undef __FUNCT__ 345 #define __FUNCT__ "SNESCreate_EQ_TR" 346 int SNESCreate_EQ_TR(SNES snes) 347 { 348 SNES_EQ_TR *neP; 349 int ierr; 350 351 PetscFunctionBegin; 352 if (snes->method_class != SNES_NONLINEAR_EQUATIONS) { 353 SETERRQ(PETSC_ERR_ARG_WRONG,"For SNES_NONLINEAR_EQUATIONS only"); 354 } 355 snes->setup = SNESSetUp_EQ_TR; 356 snes->solve = SNESSolve_EQ_TR; 357 snes->destroy = SNESDestroy_EQ_TR; 358 snes->converged = SNESConverged_EQ_TR; 359 snes->setfromoptions = SNESSetFromOptions_EQ_TR; 360 snes->view = SNESView_EQ_TR; 361 snes->nwork = 0; 362 363 ierr = PetscNew(SNES_EQ_TR,&neP);CHKERRQ(ierr); 364 PetscLogObjectMemory(snes,sizeof(SNES_EQ_TR)); 365 snes->data = (void*)neP; 366 neP->mu = 0.25; 367 neP->eta = 0.75; 368 neP->delta = 0.0; 369 neP->delta0 = 0.2; 370 neP->delta1 = 0.3; 371 neP->delta2 = 0.75; 372 neP->delta3 = 2.0; 373 neP->sigma = 0.0001; 374 neP->itflag = 0; 375 neP->rnorm0 = 0; 376 neP->ttol = 0; 377 PetscFunctionReturn(0); 378 } 379 EXTERN_C_END 380 381