1 #define PETSCSNES_DLL 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_TR_KSPConverged_Private" 11 PetscErrorCode SNES_TR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *ctx) 12 { 13 SNES snes = (SNES) ctx; 14 SNES_TR *neP = (SNES_TR*)snes->data; 15 Vec x; 16 PetscReal nrm; 17 PetscErrorCode ierr; 18 19 PetscFunctionBegin; 20 ierr = KSPDefaultConverged(ksp,n,rnorm,reason,ctx);CHKERRQ(ierr); 21 if (*reason) { 22 ierr = PetscInfo2(snes,"default convergence test KSP iterations=%D, rnorm=%G\n",n,rnorm);CHKERRQ(ierr); 23 } 24 /* Determine norm of solution */ 25 ierr = KSPBuildSolution(ksp,0,&x);CHKERRQ(ierr); 26 ierr = VecNorm(x,NORM_2,&nrm);CHKERRQ(ierr); 27 if (nrm >= neP->delta) { 28 ierr = PetscInfo2(snes,"Ending linear iteration early, delta=%G, length=%G\n",neP->delta,nrm);CHKERRQ(ierr); 29 *reason = KSP_CONVERGED_STEP_LENGTH; 30 } 31 PetscFunctionReturn(0); 32 } 33 34 /* ---------------------------------------------------------------- */ 35 #undef __FUNCT__ 36 #define __FUNCT__ "SNES_TR_Converged_Private" 37 /* 38 SNES_TR_Converged_Private -test convergence JUST for 39 the trust region tolerance. 40 41 */ 42 static PetscErrorCode SNES_TR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy) 43 { 44 SNES_TR *neP = (SNES_TR *)snes->data; 45 PetscErrorCode ierr; 46 47 PetscFunctionBegin; 48 *reason = SNES_CONVERGED_ITERATING; 49 if (neP->delta < xnorm * snes->deltatol) { 50 ierr = PetscInfo3(snes,"Converged due to trust region param %G<%G*%G\n",neP->delta,xnorm,snes->deltatol);CHKERRQ(ierr); 51 *reason = SNES_CONVERGED_TR_DELTA; 52 } else if (snes->nfuncs >= snes->max_funcs) { 53 ierr = PetscInfo1(snes,"Exceeded maximum number of function evaluations: %D\n",snes->max_funcs);CHKERRQ(ierr); 54 *reason = SNES_DIVERGED_FUNCTION_COUNT; 55 } 56 PetscFunctionReturn(0); 57 } 58 59 60 /* 61 SNESSolve_TR - Implements Newton's Method with a very simple trust 62 region approach for solving systems of nonlinear equations. 63 64 65 */ 66 #undef __FUNCT__ 67 #define __FUNCT__ "SNESSolve_TR" 68 static PetscErrorCode SNESSolve_TR(SNES snes) 69 { 70 SNES_TR *neP = (SNES_TR*)snes->data; 71 Vec X,F,Y,G,Ytmp; 72 PetscErrorCode ierr; 73 PetscInt maxits,i,lits; 74 MatStructure flg = DIFFERENT_NONZERO_PATTERN; 75 PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1; 76 PetscScalar cnorm; 77 KSP ksp; 78 SNESConvergedReason reason = SNES_CONVERGED_ITERATING; 79 PetscTruth conv,breakout = PETSC_FALSE; 80 81 PetscFunctionBegin; 82 maxits = snes->max_its; /* maximum number of iterations */ 83 X = snes->vec_sol; /* solution vector */ 84 F = snes->vec_func; /* residual vector */ 85 Y = snes->work[0]; /* work vectors */ 86 G = snes->work[1]; 87 Ytmp = snes->work[2]; 88 89 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 90 snes->iter = 0; 91 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 92 93 ierr = SNESComputeFunction(snes,X,F);CHKERRQ(ierr); /* F(X) */ 94 ierr = VecNorm(F,NORM_2,&fnorm);CHKERRQ(ierr); /* fnorm <- || F || */ 95 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 96 snes->norm = fnorm; 97 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 98 delta = neP->delta0*fnorm; 99 neP->delta = delta; 100 SNESLogConvHistory(snes,fnorm,0); 101 SNESMonitor(snes,0,fnorm); 102 103 /* set parameter for default relative tolerance convergence test */ 104 snes->ttol = fnorm*snes->rtol; 105 /* test convergence */ 106 ierr = (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);CHKERRQ(ierr); 107 if (snes->reason) PetscFunctionReturn(0); 108 109 /* Set the stopping criteria to use the More' trick. */ 110 ierr = PetscOptionsHasName(PETSC_NULL,"-snes_tr_ksp_regular_convergence_test",&conv);CHKERRQ(ierr); 111 if (!conv) { 112 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr); 113 ierr = KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,(void*)snes);CHKERRQ(ierr); 114 ierr = PetscInfo(snes,"Using Krylov convergence test SNES_TR_KSPConverged_Private\n");CHKERRQ(ierr); 115 } 116 117 for (i=0; i<maxits; i++) { 118 119 /* Call general purpose update function */ 120 if (snes->ops->update) { 121 ierr = (*snes->ops->update)(snes, snes->iter);CHKERRQ(ierr); 122 } 123 124 /* Solve J Y = F, where J is Jacobian matrix */ 125 ierr = SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);CHKERRQ(ierr); 126 ierr = KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre,flg);CHKERRQ(ierr); 127 ierr = SNES_KSPSolve(snes,snes->ksp,F,Ytmp);CHKERRQ(ierr); 128 ierr = KSPGetIterationNumber(snes->ksp,&lits);CHKERRQ(ierr); 129 snes->linear_its += lits; 130 ierr = PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);CHKERRQ(ierr); 131 ierr = VecNorm(Ytmp,NORM_2,&nrm);CHKERRQ(ierr); 132 norm1 = nrm; 133 while(1) { 134 ierr = VecCopy(Ytmp,Y);CHKERRQ(ierr); 135 nrm = norm1; 136 137 /* Scale Y if need be and predict new value of F norm */ 138 if (nrm >= delta) { 139 nrm = delta/nrm; 140 gpnorm = (1.0 - nrm)*fnorm; 141 cnorm = nrm; 142 ierr = PetscInfo1(snes,"Scaling direction by %G\n",nrm);CHKERRQ(ierr); 143 ierr = VecScale(Y,cnorm);CHKERRQ(ierr); 144 nrm = gpnorm; 145 ynorm = delta; 146 } else { 147 gpnorm = 0.0; 148 ierr = PetscInfo(snes,"Direction is in Trust Region\n");CHKERRQ(ierr); 149 ynorm = nrm; 150 } 151 ierr = VecAYPX(Y,-1.0,X);CHKERRQ(ierr); /* Y <- X - Y */ 152 ierr = VecCopy(X,snes->vec_sol_update);CHKERRQ(ierr); 153 ierr = SNESComputeFunction(snes,Y,G);CHKERRQ(ierr); /* F(X) */ 154 ierr = VecNorm(G,NORM_2,&gnorm);CHKERRQ(ierr); /* gnorm <- || g || */ 155 if (fnorm == gpnorm) rho = 0.0; 156 else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm); 157 158 /* Update size of trust region */ 159 if (rho < neP->mu) delta *= neP->delta1; 160 else if (rho < neP->eta) delta *= neP->delta2; 161 else delta *= neP->delta3; 162 ierr = PetscInfo3(snes,"fnorm=%G, gnorm=%G, ynorm=%G\n",fnorm,gnorm,ynorm);CHKERRQ(ierr); 163 ierr = PetscInfo3(snes,"gpred=%G, rho=%G, delta=%G\n",gpnorm,rho,delta);CHKERRQ(ierr); 164 neP->delta = delta; 165 if (rho > neP->sigma) break; 166 ierr = PetscInfo(snes,"Trying again in smaller region\n");CHKERRQ(ierr); 167 /* check to see if progress is hopeless */ 168 neP->itflag = PETSC_FALSE; 169 ierr = SNES_TR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 170 if (!reason) { ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); } 171 if (reason) { 172 /* We're not progressing, so return with the current iterate */ 173 SNESMonitor(snes,i+1,fnorm); 174 breakout = PETSC_TRUE; 175 break; 176 } 177 snes->numFailures++; 178 } 179 if (!breakout) { 180 /* Update function and solution vectors */ 181 fnorm = gnorm; 182 ierr = VecCopy(G,F);CHKERRQ(ierr); 183 ierr = VecCopy(Y,X);CHKERRQ(ierr); 184 /* Monitor convergence */ 185 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 186 snes->iter = i+1; 187 snes->norm = fnorm; 188 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 189 SNESLogConvHistory(snes,snes->norm,lits); 190 SNESMonitor(snes,snes->iter,snes->norm); 191 /* Test for convergence, xnorm = || X || */ 192 neP->itflag = PETSC_TRUE; 193 if (snes->ops->converged != SNESSkipConverged) { ierr = VecNorm(X,NORM_2,&xnorm);CHKERRQ(ierr); } 194 ierr = (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);CHKERRQ(ierr); 195 if (reason) break; 196 } else { 197 break; 198 } 199 } 200 if (i == maxits) { 201 ierr = PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);CHKERRQ(ierr); 202 if (!reason) reason = SNES_DIVERGED_MAX_IT; 203 } 204 ierr = PetscObjectTakeAccess(snes);CHKERRQ(ierr); 205 snes->reason = reason; 206 ierr = PetscObjectGrantAccess(snes);CHKERRQ(ierr); 207 PetscFunctionReturn(0); 208 } 209 /*------------------------------------------------------------*/ 210 #undef __FUNCT__ 211 #define __FUNCT__ "SNESSetUp_TR" 212 static PetscErrorCode SNESSetUp_TR(SNES snes) 213 { 214 PetscErrorCode ierr; 215 216 PetscFunctionBegin; 217 if (!snes->vec_sol_update) { 218 ierr = VecDuplicate(snes->vec_sol,&snes->vec_sol_update);CHKERRQ(ierr); 219 ierr = PetscLogObjectParent(snes,snes->vec_sol_update);CHKERRQ(ierr); 220 } 221 if (!snes->work) { 222 snes->nwork = 3; 223 ierr = VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);CHKERRQ(ierr); 224 ierr = PetscLogObjectParents(snes,snes->nwork,snes->work);CHKERRQ(ierr); 225 } 226 PetscFunctionReturn(0); 227 } 228 /*------------------------------------------------------------*/ 229 #undef __FUNCT__ 230 #define __FUNCT__ "SNESDestroy_TR" 231 static PetscErrorCode SNESDestroy_TR(SNES snes) 232 { 233 PetscErrorCode ierr; 234 235 PetscFunctionBegin; 236 if (snes->vec_sol_update) { 237 ierr = VecDestroy(snes->vec_sol_update);CHKERRQ(ierr); 238 snes->vec_sol_update = PETSC_NULL; 239 } 240 if (snes->nwork) { 241 ierr = VecDestroyVecs(snes->work,snes->nwork);CHKERRQ(ierr); 242 snes->nwork = 0; 243 snes->work = PETSC_NULL; 244 } 245 ierr = PetscFree(snes->data);CHKERRQ(ierr); 246 PetscFunctionReturn(0); 247 } 248 /*------------------------------------------------------------*/ 249 250 #undef __FUNCT__ 251 #define __FUNCT__ "SNESSetFromOptions_TR" 252 static PetscErrorCode SNESSetFromOptions_TR(SNES snes) 253 { 254 SNES_TR *ctx = (SNES_TR *)snes->data; 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 ierr = PetscOptionsHead("SNES trust region options for nonlinear equations");CHKERRQ(ierr); 259 ierr = PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,0);CHKERRQ(ierr); 260 ierr = PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,0);CHKERRQ(ierr); 261 ierr = PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,0);CHKERRQ(ierr); 262 ierr = PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,0);CHKERRQ(ierr); 263 ierr = PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,0);CHKERRQ(ierr); 264 ierr = PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,0);CHKERRQ(ierr); 265 ierr = PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,0);CHKERRQ(ierr); 266 ierr = PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,0);CHKERRQ(ierr); 267 ierr = PetscOptionsTail();CHKERRQ(ierr); 268 PetscFunctionReturn(0); 269 } 270 271 #undef __FUNCT__ 272 #define __FUNCT__ "SNESView_TR" 273 static PetscErrorCode SNESView_TR(SNES snes,PetscViewer viewer) 274 { 275 SNES_TR *tr = (SNES_TR *)snes->data; 276 PetscErrorCode ierr; 277 PetscTruth iascii; 278 279 PetscFunctionBegin; 280 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 281 if (iascii) { 282 ierr = PetscViewerASCIIPrintf(viewer," mu=%G, eta=%G, sigma=%G\n",tr->mu,tr->eta,tr->sigma);CHKERRQ(ierr); 283 ierr = PetscViewerASCIIPrintf(viewer," delta0=%G, delta1=%G, delta2=%G, delta3=%G\n",tr->delta0,tr->delta1,tr->delta2,tr->delta3);CHKERRQ(ierr); 284 } else { 285 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name); 286 } 287 PetscFunctionReturn(0); 288 } 289 /* ------------------------------------------------------------ */ 290 /*MC 291 SNESTR - Newton based nonlinear solver that uses a trust region 292 293 Options Database: 294 + -snes_trtol <tol> Trust region tolerance 295 . -snes_tr_mu <mu> 296 . -snes_tr_eta <eta> 297 . -snes_tr_sigma <sigma> 298 . -snes_tr_delta0 <delta0> 299 . -snes_tr_delta1 <delta1> 300 . -snes_tr_delta2 <delta2> 301 - -snes_tr_delta3 <delta3> 302 303 The basic algorithm is taken from "The Minpack Project", by More', 304 Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development 305 of Mathematical Software", Wayne Cowell, editor. 306 307 This is intended as a model implementation, since it does not 308 necessarily have many of the bells and whistles of other 309 implementations. 310 311 Level: intermediate 312 313 .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESSetTrustRegionTolerance() 314 315 M*/ 316 EXTERN_C_BEGIN 317 #undef __FUNCT__ 318 #define __FUNCT__ "SNESCreate_TR" 319 PetscErrorCode PETSCSNES_DLLEXPORT SNESCreate_TR(SNES snes) 320 { 321 SNES_TR *neP; 322 PetscErrorCode ierr; 323 324 PetscFunctionBegin; 325 snes->ops->setup = SNESSetUp_TR; 326 snes->ops->solve = SNESSolve_TR; 327 snes->ops->destroy = SNESDestroy_TR; 328 snes->ops->setfromoptions = SNESSetFromOptions_TR; 329 snes->ops->view = SNESView_TR; 330 331 ierr = PetscNewLog(snes,SNES_TR,&neP);CHKERRQ(ierr); 332 snes->data = (void*)neP; 333 neP->mu = 0.25; 334 neP->eta = 0.75; 335 neP->delta = 0.0; 336 neP->delta0 = 0.2; 337 neP->delta1 = 0.3; 338 neP->delta2 = 0.75; 339 neP->delta3 = 2.0; 340 neP->sigma = 0.0001; 341 neP->itflag = PETSC_FALSE; 342 neP->rnorm0 = 0; 343 neP->ttol = 0; 344 PetscFunctionReturn(0); 345 } 346 EXTERN_C_END 347 348