Actual source code: tr.c
petsc-3.8.4 2018-03-24
2: #include <../src/snes/impls/tr/trimpl.h>
4: typedef struct {
5: void *ctx;
6: SNES snes;
7: } SNES_TR_KSPConverged_Ctx;
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: PetscErrorCode SNES_TR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx)
14: {
15: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
16: SNES snes = ctx->snes;
17: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
18: Vec x;
19: PetscReal nrm;
20: PetscErrorCode ierr;
23: KSPConvergedDefault(ksp,n,rnorm,reason,ctx->ctx);
24: if (*reason) {
25: PetscInfo2(snes,"default convergence test KSP iterations=%D, rnorm=%g\n",n,(double)rnorm);
26: }
27: /* Determine norm of solution */
28: KSPBuildSolution(ksp,0,&x);
29: VecNorm(x,NORM_2,&nrm);
30: if (nrm >= neP->delta) {
31: PetscInfo2(snes,"Ending linear iteration early, delta=%g, length=%g\n",(double)neP->delta,(double)nrm);
32: *reason = KSP_CONVERGED_STEP_LENGTH;
33: }
34: return(0);
35: }
37: PetscErrorCode SNES_TR_KSPConverged_Destroy(void *cctx)
38: {
39: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
40: PetscErrorCode ierr;
43: KSPConvergedDefaultDestroy(ctx->ctx);
44: PetscFree(ctx);
45: return(0);
46: }
48: /* ---------------------------------------------------------------- */
49: /*
50: SNES_TR_Converged_Private -test convergence JUST for
51: the trust region tolerance.
53: */
54: static PetscErrorCode SNES_TR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
55: {
56: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
60: *reason = SNES_CONVERGED_ITERATING;
61: if (neP->delta < xnorm * snes->deltatol) {
62: PetscInfo3(snes,"Converged due to trust region param %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol);
63: *reason = SNES_CONVERGED_TR_DELTA;
64: } else if (snes->nfuncs >= snes->max_funcs) {
65: PetscInfo1(snes,"Exceeded maximum number of function evaluations: %D\n",snes->max_funcs);
66: *reason = SNES_DIVERGED_FUNCTION_COUNT;
67: }
68: return(0);
69: }
72: /*
73: SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust
74: region approach for solving systems of nonlinear equations.
77: */
78: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
79: {
80: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
81: Vec X,F,Y,G,Ytmp;
82: PetscErrorCode ierr;
83: PetscInt maxits,i,lits;
84: PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1;
85: PetscScalar cnorm;
86: KSP ksp;
87: SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
88: PetscBool conv = PETSC_FALSE,breakout = PETSC_FALSE;
91: 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);
93: maxits = snes->max_its; /* maximum number of iterations */
94: X = snes->vec_sol; /* solution vector */
95: F = snes->vec_func; /* residual vector */
96: Y = snes->work[0]; /* work vectors */
97: G = snes->work[1];
98: Ytmp = snes->work[2];
100: PetscObjectSAWsTakeAccess((PetscObject)snes);
101: snes->iter = 0;
102: PetscObjectSAWsGrantAccess((PetscObject)snes);
104: if (!snes->vec_func_init_set) {
105: SNESComputeFunction(snes,X,F); /* F(X) */
106: } else snes->vec_func_init_set = PETSC_FALSE;
108: VecNorm(F,NORM_2,&fnorm); /* fnorm <- || F || */
109: SNESCheckFunctionNorm(snes,fnorm);
110: VecNorm(X,NORM_2,&xnorm); /* fnorm <- || F || */
111: PetscObjectSAWsTakeAccess((PetscObject)snes);
112: snes->norm = fnorm;
113: PetscObjectSAWsGrantAccess((PetscObject)snes);
114: delta = xnorm ? neP->delta0*xnorm : neP->delta0;
115: neP->delta = delta;
116: SNESLogConvergenceHistory(snes,fnorm,0);
117: SNESMonitor(snes,0,fnorm);
119: /* test convergence */
120: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
121: if (snes->reason) return(0);
123: /* Set the stopping criteria to use the More' trick. */
124: PetscOptionsGetBool(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_tr_ksp_regular_convergence_test",&conv,NULL);
125: if (!conv) {
126: SNES_TR_KSPConverged_Ctx *ctx;
127: SNESGetKSP(snes,&ksp);
128: PetscNew(&ctx);
129: ctx->snes = snes;
130: KSPConvergedDefaultCreate(&ctx->ctx);
131: KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,ctx,SNES_TR_KSPConverged_Destroy);
132: PetscInfo(snes,"Using Krylov convergence test SNES_TR_KSPConverged_Private\n");
133: }
135: for (i=0; i<maxits; i++) {
137: /* Call general purpose update function */
138: if (snes->ops->update) {
139: (*snes->ops->update)(snes, snes->iter);
140: }
142: /* Solve J Y = F, where J is Jacobian matrix */
143: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
144: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
145: KSPSolve(snes->ksp,F,Ytmp);
146: KSPGetIterationNumber(snes->ksp,&lits);
148: snes->linear_its += lits;
150: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
151: VecNorm(Ytmp,NORM_2,&nrm);
152: norm1 = nrm;
153: while (1) {
154: VecCopy(Ytmp,Y);
155: nrm = norm1;
157: /* Scale Y if need be and predict new value of F norm */
158: if (nrm >= delta) {
159: nrm = delta/nrm;
160: gpnorm = (1.0 - nrm)*fnorm;
161: cnorm = nrm;
162: PetscInfo1(snes,"Scaling direction by %g\n",(double)nrm);
163: VecScale(Y,cnorm);
164: nrm = gpnorm;
165: ynorm = delta;
166: } else {
167: gpnorm = 0.0;
168: PetscInfo(snes,"Direction is in Trust Region\n");
169: ynorm = nrm;
170: }
171: VecAYPX(Y,-1.0,X); /* Y <- X - Y */
172: VecCopy(X,snes->vec_sol_update);
173: SNESComputeFunction(snes,Y,G); /* F(X) */
174: VecNorm(G,NORM_2,&gnorm); /* gnorm <- || g || */
175: if (fnorm == gpnorm) rho = 0.0;
176: else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);
178: /* Update size of trust region */
179: if (rho < neP->mu) delta *= neP->delta1;
180: else if (rho < neP->eta) delta *= neP->delta2;
181: else delta *= neP->delta3;
182: PetscInfo3(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm);
183: PetscInfo3(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta);
185: neP->delta = delta;
186: if (rho > neP->sigma) break;
187: PetscInfo(snes,"Trying again in smaller region\n");
188: /* check to see if progress is hopeless */
189: neP->itflag = PETSC_FALSE;
190: SNES_TR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
191: if (!reason) { (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP); }
192: if (reason) {
193: /* We're not progressing, so return with the current iterate */
194: SNESMonitor(snes,i+1,fnorm);
195: breakout = PETSC_TRUE;
196: break;
197: }
198: snes->numFailures++;
199: }
200: if (!breakout) {
201: /* Update function and solution vectors */
202: fnorm = gnorm;
203: VecCopy(G,F);
204: VecCopy(Y,X);
205: /* Monitor convergence */
206: PetscObjectSAWsTakeAccess((PetscObject)snes);
207: snes->iter = i+1;
208: snes->norm = fnorm;
209: PetscObjectSAWsGrantAccess((PetscObject)snes);
210: SNESLogConvergenceHistory(snes,snes->norm,lits);
211: SNESMonitor(snes,snes->iter,snes->norm);
212: /* Test for convergence, xnorm = || X || */
213: neP->itflag = PETSC_TRUE;
214: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
215: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
216: if (reason) break;
217: } else break;
218: }
219: if (i == maxits) {
220: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
221: if (!reason) reason = SNES_DIVERGED_MAX_IT;
222: }
223: PetscObjectSAWsTakeAccess((PetscObject)snes);
224: snes->reason = reason;
225: PetscObjectSAWsGrantAccess((PetscObject)snes);
226: return(0);
227: }
228: /*------------------------------------------------------------*/
229: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
230: {
234: SNESSetWorkVecs(snes,3);
235: SNESSetUpMatrices(snes);
236: return(0);
237: }
239: PetscErrorCode SNESReset_NEWTONTR(SNES snes)
240: {
243: return(0);
244: }
246: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
247: {
251: SNESReset_NEWTONTR(snes);
252: PetscFree(snes->data);
253: return(0);
254: }
255: /*------------------------------------------------------------*/
257: static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptionItems *PetscOptionsObject,SNES snes)
258: {
259: SNES_NEWTONTR *ctx = (SNES_NEWTONTR*)snes->data;
263: PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");
264: PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);
265: PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL);
266: PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL);
267: PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL);
268: PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);
269: PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL);
270: PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL);
271: PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL);
272: PetscOptionsTail();
273: return(0);
274: }
276: static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer)
277: {
278: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
280: PetscBool iascii;
283: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
284: if (iascii) {
285: PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma);
286: PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3);
287: }
288: return(0);
289: }
290: /* ------------------------------------------------------------ */
291: /*MC
292: SNESNEWTONTR - Newton based nonlinear solver that uses a trust region
294: Options Database:
295: + -snes_trtol <tol> - trust region tolerance
296: . -snes_tr_mu <mu> - trust region parameter
297: . -snes_tr_eta <eta> - trust region parameter
298: . -snes_tr_sigma <sigma> - trust region parameter
299: . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x)
300: . -snes_tr_delta1 <delta1> - trust region parameter
301: . -snes_tr_delta2 <delta2> - trust region parameter
302: - -snes_tr_delta3 <delta3> - trust region parameter
304: The basic algorithm is taken from "The Minpack Project", by More',
305: Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development
306: of Mathematical Software", Wayne Cowell, editor.
308: Level: intermediate
310: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance()
312: M*/
313: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
314: {
315: SNES_NEWTONTR *neP;
319: snes->ops->setup = SNESSetUp_NEWTONTR;
320: snes->ops->solve = SNESSolve_NEWTONTR;
321: snes->ops->destroy = SNESDestroy_NEWTONTR;
322: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
323: snes->ops->view = SNESView_NEWTONTR;
324: snes->ops->reset = SNESReset_NEWTONTR;
326: snes->usesksp = PETSC_TRUE;
327: snes->usesnpc = PETSC_FALSE;
329: snes->alwayscomputesfinalresidual = PETSC_TRUE;
331: PetscNewLog(snes,&neP);
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.0;
343: neP->ttol = 0.0;
344: return(0);
345: }