Actual source code: tr.c
petsc-3.11.4 2019-09-28
2: #include <../src/snes/impls/tr/trimpl.h>
4: typedef struct {
5: SNES snes;
6: /* Information on the regular SNES convergence test; which may have been user provided */
7: PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
8: PetscErrorCode (*convdestroy)(void*);
9: void *convctx;
10: } SNES_TR_KSPConverged_Ctx;
12: static PetscErrorCode SNESTR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx)
13: {
14: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
15: SNES snes = ctx->snes;
16: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
17: Vec x;
18: PetscReal nrm;
19: PetscErrorCode ierr;
22: (*ctx->convtest)(ksp,n,rnorm,reason,ctx->convctx);
23: if (*reason) {
24: PetscInfo2(snes,"Default or user provided convergence test KSP iterations=%D, rnorm=%g\n",n,(double)rnorm);
25: }
26: /* Determine norm of solution */
27: KSPBuildSolution(ksp,0,&x);
28: VecNorm(x,NORM_2,&nrm);
29: if (nrm >= neP->delta) {
30: PetscInfo2(snes,"Ending linear iteration early, delta=%g, length=%g\n",(double)neP->delta,(double)nrm);
31: *reason = KSP_CONVERGED_STEP_LENGTH;
32: }
33: return(0);
34: }
36: static PetscErrorCode SNESTR_KSPConverged_Destroy(void *cctx)
37: {
38: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
39: PetscErrorCode ierr;
42: (*ctx->convdestroy)(ctx->convctx);
43: PetscFree(ctx);
44: return(0);
45: }
47: /* ---------------------------------------------------------------- */
48: /*
49: SNESTR_Converged_Private -test convergence JUST for
50: the trust region tolerance.
52: */
53: static PetscErrorCode SNESTR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
54: {
55: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
59: *reason = SNES_CONVERGED_ITERATING;
60: if (neP->delta < xnorm * snes->deltatol) {
61: PetscInfo3(snes,"Converged due to trust region param %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol);
62: *reason = SNES_CONVERGED_TR_DELTA;
63: } else if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
64: PetscInfo1(snes,"Exceeded maximum number of function evaluations: %D\n",snes->max_funcs);
65: *reason = SNES_DIVERGED_FUNCTION_COUNT;
66: }
67: return(0);
68: }
71: /*
72: SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust
73: region approach for solving systems of nonlinear equations.
76: */
77: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
78: {
79: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
80: Vec X,F,Y,G,Ytmp;
81: PetscErrorCode ierr;
82: PetscInt maxits,i,lits;
83: PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1;
84: PetscScalar cnorm;
85: KSP ksp;
86: SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
87: PetscBool breakout = PETSC_FALSE;
88: SNES_TR_KSPConverged_Ctx *ctx;
89: PetscErrorCode (*convtest)(KSP,PetscInt,PetscReal,KSPConvergedReason*,void*);
92: if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
94: maxits = snes->max_its; /* maximum number of iterations */
95: X = snes->vec_sol; /* solution vector */
96: F = snes->vec_func; /* residual vector */
97: Y = snes->work[0]; /* work vectors */
98: G = snes->work[1];
99: Ytmp = snes->work[2];
101: PetscObjectSAWsTakeAccess((PetscObject)snes);
102: snes->iter = 0;
103: PetscObjectSAWsGrantAccess((PetscObject)snes);
105: /* Set the linear stopping criteria to use the More' trick. */
106: SNESGetKSP(snes,&ksp);
107: KSPGetConvergenceTest(ksp,&convtest,NULL,NULL);
108: if (convtest != SNESTR_KSPConverged_Private) {
109: PetscNew(&ctx);
110: ctx->snes = snes;
111: KSPGetAndClearConvergenceTest(ksp,&ctx->convtest,&ctx->convctx,&ctx->convdestroy);
112: KSPSetConvergenceTest(ksp,SNESTR_KSPConverged_Private,ctx,SNESTR_KSPConverged_Destroy);
113: PetscInfo(snes,"Using Krylov convergence test SNESTR_KSPConverged_Private\n");
114: }
116: if (!snes->vec_func_init_set) {
117: SNESComputeFunction(snes,X,F); /* F(X) */
118: } else snes->vec_func_init_set = PETSC_FALSE;
120: VecNorm(F,NORM_2,&fnorm); /* fnorm <- || F || */
121: SNESCheckFunctionNorm(snes,fnorm);
122: VecNorm(X,NORM_2,&xnorm); /* fnorm <- || F || */
123: PetscObjectSAWsTakeAccess((PetscObject)snes);
124: snes->norm = fnorm;
125: PetscObjectSAWsGrantAccess((PetscObject)snes);
126: delta = xnorm ? neP->delta0*xnorm : neP->delta0;
127: neP->delta = delta;
128: SNESLogConvergenceHistory(snes,fnorm,0);
129: SNESMonitor(snes,0,fnorm);
131: /* test convergence */
132: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
133: if (snes->reason) return(0);
136: for (i=0; i<maxits; i++) {
138: /* Call general purpose update function */
139: if (snes->ops->update) {
140: (*snes->ops->update)(snes, snes->iter);
141: }
143: /* Solve J Y = F, where J is Jacobian matrix */
144: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
145: SNESCheckJacobianDomainerror(snes);
146: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
147: KSPSolve(snes->ksp,F,Ytmp);
148: KSPGetIterationNumber(snes->ksp,&lits);
150: snes->linear_its += lits;
152: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
153: VecNorm(Ytmp,NORM_2,&nrm);
154: norm1 = nrm;
155: while (1) {
156: VecCopy(Ytmp,Y);
157: nrm = norm1;
159: /* Scale Y if need be and predict new value of F norm */
160: if (nrm >= delta) {
161: nrm = delta/nrm;
162: gpnorm = (1.0 - nrm)*fnorm;
163: cnorm = nrm;
164: PetscInfo1(snes,"Scaling direction by %g\n",(double)nrm);
165: VecScale(Y,cnorm);
166: nrm = gpnorm;
167: ynorm = delta;
168: } else {
169: gpnorm = 0.0;
170: PetscInfo(snes,"Direction is in Trust Region\n");
171: ynorm = nrm;
172: }
173: VecCopy(Y,snes->vec_sol_update);
174: VecAYPX(Y,-1.0,X); /* Y <- X - Y */
175: SNESComputeFunction(snes,Y,G); /* F(X) */
176: VecNorm(G,NORM_2,&gnorm); /* gnorm <- || g || */
177: SNESCheckFunctionNorm(snes,gnorm);
178: if (fnorm == gpnorm) rho = 0.0;
179: else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);
181: /* Update size of trust region */
182: if (rho < neP->mu) delta *= neP->delta1;
183: else if (rho < neP->eta) delta *= neP->delta2;
184: else delta *= neP->delta3;
185: PetscInfo3(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm);
186: PetscInfo3(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta);
188: neP->delta = delta;
189: if (rho > neP->sigma) break;
190: PetscInfo(snes,"Trying again in smaller region\n");
191: /* check to see if progress is hopeless */
192: neP->itflag = PETSC_FALSE;
193: SNESTR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
194: if (!reason) { (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP); }
195: if (reason) {
196: /* We're not progressing, so return with the current iterate */
197: SNESMonitor(snes,i+1,fnorm);
198: breakout = PETSC_TRUE;
199: break;
200: }
201: snes->numFailures++;
202: }
203: if (!breakout) {
204: /* Update function and solution vectors */
205: fnorm = gnorm;
206: VecCopy(G,F);
207: VecCopy(Y,X);
208: /* Monitor convergence */
209: PetscObjectSAWsTakeAccess((PetscObject)snes);
210: snes->iter = i+1;
211: snes->norm = fnorm;
212: snes->xnorm = xnorm;
213: snes->ynorm = ynorm;
214: PetscObjectSAWsGrantAccess((PetscObject)snes);
215: SNESLogConvergenceHistory(snes,snes->norm,lits);
216: SNESMonitor(snes,snes->iter,snes->norm);
217: /* Test for convergence, xnorm = || X || */
218: neP->itflag = PETSC_TRUE;
219: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
220: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
221: if (reason) break;
222: } else break;
223: }
224: if (i == maxits) {
225: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
226: if (!reason) reason = SNES_DIVERGED_MAX_IT;
227: }
228: PetscObjectSAWsTakeAccess((PetscObject)snes);
229: snes->reason = reason;
230: PetscObjectSAWsGrantAccess((PetscObject)snes);
231: return(0);
232: }
233: /*------------------------------------------------------------*/
234: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
235: {
239: SNESSetWorkVecs(snes,3);
240: SNESSetUpMatrices(snes);
241: return(0);
242: }
244: PetscErrorCode SNESReset_NEWTONTR(SNES snes)
245: {
248: return(0);
249: }
251: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
252: {
256: SNESReset_NEWTONTR(snes);
257: PetscFree(snes->data);
258: return(0);
259: }
260: /*------------------------------------------------------------*/
262: static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes)
263: {
264: SNES_NEWTONTR *ctx = (SNES_NEWTONTR*)snes->data;
268: PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");
269: PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);
270: PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL);
271: PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL);
272: PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL);
273: PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);
274: PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL);
275: PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL);
276: PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL);
277: PetscOptionsTail();
278: return(0);
279: }
281: static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer)
282: {
283: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
285: PetscBool iascii;
288: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
289: if (iascii) {
290: PetscViewerASCIIPrintf(viewer," Trust region tolerance (-snes_trtol)\n",(double)snes->deltatol);
291: PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma);
292: PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3);
293: }
294: return(0);
295: }
296: /* ------------------------------------------------------------ */
297: /*MC
298: SNESNEWTONTR - Newton based nonlinear solver that uses a trust region
300: Options Database:
301: + -snes_trtol <tol> - trust region tolerance
302: . -snes_tr_mu <mu> - trust region parameter
303: . -snes_tr_eta <eta> - trust region parameter
304: . -snes_tr_sigma <sigma> - trust region parameter
305: . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x)
306: . -snes_tr_delta1 <delta1> - trust region parameter
307: . -snes_tr_delta2 <delta2> - trust region parameter
308: - -snes_tr_delta3 <delta3> - trust region parameter
310: The basic algorithm is taken from "The Minpack Project", by More',
311: Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development
312: of Mathematical Software", Wayne Cowell, editor.
314: Level: intermediate
316: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance()
318: M*/
319: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
320: {
321: SNES_NEWTONTR *neP;
325: snes->ops->setup = SNESSetUp_NEWTONTR;
326: snes->ops->solve = SNESSolve_NEWTONTR;
327: snes->ops->destroy = SNESDestroy_NEWTONTR;
328: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
329: snes->ops->view = SNESView_NEWTONTR;
330: snes->ops->reset = SNESReset_NEWTONTR;
332: snes->usesksp = PETSC_TRUE;
333: snes->usesnpc = PETSC_FALSE;
335: snes->alwayscomputesfinalresidual = PETSC_TRUE;
337: PetscNewLog(snes,&neP);
338: snes->data = (void*)neP;
339: neP->mu = 0.25;
340: neP->eta = 0.75;
341: neP->delta = 0.0;
342: neP->delta0 = 0.2;
343: neP->delta1 = 0.3;
344: neP->delta2 = 0.75;
345: neP->delta3 = 2.0;
346: neP->sigma = 0.0001;
347: neP->itflag = PETSC_FALSE;
348: neP->rnorm0 = 0.0;
349: neP->ttol = 0.0;
350: return(0);
351: }