Actual source code: tr.c
petsc-3.6.1 2015-08-06
2: #include <../src/snes/impls/tr/trimpl.h> /*I "petscsnes.h" I*/
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: */
15: PetscErrorCode SNES_TR_KSPConverged_Private(KSP ksp,PetscInt n,PetscReal rnorm,KSPConvergedReason *reason,void *cctx)
16: {
17: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
18: SNES snes = ctx->snes;
19: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
20: Vec x;
21: PetscReal nrm;
22: PetscErrorCode ierr;
25: KSPConvergedDefault(ksp,n,rnorm,reason,ctx->ctx);
26: if (*reason) {
27: PetscInfo2(snes,"default convergence test KSP iterations=%D, rnorm=%g\n",n,(double)rnorm);
28: }
29: /* Determine norm of solution */
30: KSPBuildSolution(ksp,0,&x);
31: VecNorm(x,NORM_2,&nrm);
32: if (nrm >= neP->delta) {
33: PetscInfo2(snes,"Ending linear iteration early, delta=%g, length=%g\n",(double)neP->delta,(double)nrm);
34: *reason = KSP_CONVERGED_STEP_LENGTH;
35: }
36: return(0);
37: }
41: PetscErrorCode SNES_TR_KSPConverged_Destroy(void *cctx)
42: {
43: SNES_TR_KSPConverged_Ctx *ctx = (SNES_TR_KSPConverged_Ctx*)cctx;
44: PetscErrorCode ierr;
47: KSPConvergedDefaultDestroy(ctx->ctx);
48: PetscFree(ctx);
49: return(0);
50: }
52: /* ---------------------------------------------------------------- */
55: /*
56: SNES_TR_Converged_Private -test convergence JUST for
57: the trust region tolerance.
59: */
60: static PetscErrorCode SNES_TR_Converged_Private(SNES snes,PetscInt it,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
61: {
62: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
66: *reason = SNES_CONVERGED_ITERATING;
67: if (neP->delta < xnorm * snes->deltatol) {
68: PetscInfo3(snes,"Converged due to trust region param %g<%g*%g\n",(double)neP->delta,(double)xnorm,(double)snes->deltatol);
69: *reason = SNES_CONVERGED_TR_DELTA;
70: } else if (snes->nfuncs >= snes->max_funcs) {
71: PetscInfo1(snes,"Exceeded maximum number of function evaluations: %D\n",snes->max_funcs);
72: *reason = SNES_DIVERGED_FUNCTION_COUNT;
73: }
74: return(0);
75: }
78: /*
79: SNESSolve_NEWTONTR - Implements Newton's Method with a very simple trust
80: region approach for solving systems of nonlinear equations.
83: */
86: static PetscErrorCode SNESSolve_NEWTONTR(SNES snes)
87: {
88: SNES_NEWTONTR *neP = (SNES_NEWTONTR*)snes->data;
89: Vec X,F,Y,G,Ytmp;
90: PetscErrorCode ierr;
91: PetscInt maxits,i,lits;
92: PetscReal rho,fnorm,gnorm,gpnorm,xnorm=0,delta,nrm,ynorm,norm1;
93: PetscScalar cnorm;
94: KSP ksp;
95: SNESConvergedReason reason = SNES_CONVERGED_ITERATING;
96: PetscBool conv = PETSC_FALSE,breakout = PETSC_FALSE;
100: if (snes->xl || snes->xu || snes->ops->computevariablebounds) {
101: SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
102: }
104: maxits = snes->max_its; /* maximum number of iterations */
105: X = snes->vec_sol; /* solution vector */
106: F = snes->vec_func; /* residual vector */
107: Y = snes->work[0]; /* work vectors */
108: G = snes->work[1];
109: Ytmp = snes->work[2];
111: PetscObjectSAWsTakeAccess((PetscObject)snes);
112: snes->iter = 0;
113: PetscObjectSAWsGrantAccess((PetscObject)snes);
115: if (!snes->vec_func_init_set) {
116: SNESComputeFunction(snes,X,F); /* F(X) */
117: } else snes->vec_func_init_set = PETSC_FALSE;
119: VecNorm(F,NORM_2,&fnorm); /* fnorm <- || F || */
120: SNESCheckFunctionNorm(snes,fnorm);
121: VecNorm(X,NORM_2,&xnorm); /* fnorm <- || F || */
122: PetscObjectSAWsTakeAccess((PetscObject)snes);
123: snes->norm = fnorm;
124: PetscObjectSAWsGrantAccess((PetscObject)snes);
125: delta = xnorm ? neP->delta0*xnorm : neP->delta0;
126: neP->delta = delta;
127: SNESLogConvergenceHistory(snes,fnorm,0);
128: SNESMonitor(snes,0,fnorm);
130: /* test convergence */
131: (*snes->ops->converged)(snes,snes->iter,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
132: if (snes->reason) return(0);
134: /* Set the stopping criteria to use the More' trick. */
135: PetscOptionsGetBool(NULL,"-snes_tr_ksp_regular_convergence_test",&conv,NULL);
136: if (!conv) {
137: SNES_TR_KSPConverged_Ctx *ctx;
138: SNESGetKSP(snes,&ksp);
139: PetscNew(&ctx);
140: ctx->snes = snes;
141: KSPConvergedDefaultCreate(&ctx->ctx);
142: KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,ctx,SNES_TR_KSPConverged_Destroy);
143: PetscInfo(snes,"Using Krylov convergence test SNES_TR_KSPConverged_Private\n");
144: }
146: for (i=0; i<maxits; i++) {
148: /* Call general purpose update function */
149: if (snes->ops->update) {
150: (*snes->ops->update)(snes, snes->iter);
151: }
153: /* Solve J Y = F, where J is Jacobian matrix */
154: SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
155: KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
156: KSPSolve(snes->ksp,F,Ytmp);
157: KSPGetIterationNumber(snes->ksp,&lits);
159: snes->linear_its += lits;
161: PetscInfo2(snes,"iter=%D, linear solve iterations=%D\n",snes->iter,lits);
162: VecNorm(Ytmp,NORM_2,&nrm);
163: norm1 = nrm;
164: while (1) {
165: VecCopy(Ytmp,Y);
166: nrm = norm1;
168: /* Scale Y if need be and predict new value of F norm */
169: if (nrm >= delta) {
170: nrm = delta/nrm;
171: gpnorm = (1.0 - nrm)*fnorm;
172: cnorm = nrm;
173: PetscInfo1(snes,"Scaling direction by %g\n",(double)nrm);
174: VecScale(Y,cnorm);
175: nrm = gpnorm;
176: ynorm = delta;
177: } else {
178: gpnorm = 0.0;
179: PetscInfo(snes,"Direction is in Trust Region\n");
180: ynorm = nrm;
181: }
182: VecAYPX(Y,-1.0,X); /* Y <- X - Y */
183: VecCopy(X,snes->vec_sol_update);
184: SNESComputeFunction(snes,Y,G); /* F(X) */
185: VecNorm(G,NORM_2,&gnorm); /* gnorm <- || g || */
186: if (fnorm == gpnorm) rho = 0.0;
187: else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);
189: /* Update size of trust region */
190: if (rho < neP->mu) delta *= neP->delta1;
191: else if (rho < neP->eta) delta *= neP->delta2;
192: else delta *= neP->delta3;
193: PetscInfo3(snes,"fnorm=%g, gnorm=%g, ynorm=%g\n",(double)fnorm,(double)gnorm,(double)ynorm);
194: PetscInfo3(snes,"gpred=%g, rho=%g, delta=%g\n",(double)gpnorm,(double)rho,(double)delta);
196: neP->delta = delta;
197: if (rho > neP->sigma) break;
198: PetscInfo(snes,"Trying again in smaller region\n");
199: /* check to see if progress is hopeless */
200: neP->itflag = PETSC_FALSE;
201: SNES_TR_Converged_Private(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
202: if (!reason) { (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP); }
203: if (reason) {
204: /* We're not progressing, so return with the current iterate */
205: SNESMonitor(snes,i+1,fnorm);
206: breakout = PETSC_TRUE;
207: break;
208: }
209: snes->numFailures++;
210: }
211: if (!breakout) {
212: /* Update function and solution vectors */
213: fnorm = gnorm;
214: VecCopy(G,F);
215: VecCopy(Y,X);
216: /* Monitor convergence */
217: PetscObjectSAWsTakeAccess((PetscObject)snes);
218: snes->iter = i+1;
219: snes->norm = fnorm;
220: PetscObjectSAWsGrantAccess((PetscObject)snes);
221: SNESLogConvergenceHistory(snes,snes->norm,lits);
222: SNESMonitor(snes,snes->iter,snes->norm);
223: /* Test for convergence, xnorm = || X || */
224: neP->itflag = PETSC_TRUE;
225: if (snes->ops->converged != SNESConvergedSkip) { VecNorm(X,NORM_2,&xnorm); }
226: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&reason,snes->cnvP);
227: if (reason) break;
228: } else break;
229: }
230: if (i == maxits) {
231: PetscInfo1(snes,"Maximum number of iterations has been reached: %D\n",maxits);
232: if (!reason) reason = SNES_DIVERGED_MAX_IT;
233: }
234: PetscObjectSAWsTakeAccess((PetscObject)snes);
235: snes->reason = reason;
236: PetscObjectSAWsGrantAccess((PetscObject)snes);
237: return(0);
238: }
239: /*------------------------------------------------------------*/
242: static PetscErrorCode SNESSetUp_NEWTONTR(SNES snes)
243: {
247: SNESSetWorkVecs(snes,3);
248: SNESSetUpMatrices(snes);
249: return(0);
250: }
254: PetscErrorCode SNESReset_NEWTONTR(SNES snes)
255: {
258: return(0);
259: }
263: static PetscErrorCode SNESDestroy_NEWTONTR(SNES snes)
264: {
268: SNESReset_NEWTONTR(snes);
269: PetscFree(snes->data);
270: return(0);
271: }
272: /*------------------------------------------------------------*/
276: static PetscErrorCode SNESSetFromOptions_NEWTONTR(PetscOptions *PetscOptionsObject,SNES snes)
277: {
278: SNES_NEWTONTR *ctx = (SNES_NEWTONTR*)snes->data;
282: PetscOptionsHead(PetscOptionsObject,"SNES trust region options for nonlinear equations");
283: PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,NULL);
284: PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,NULL);
285: PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,NULL);
286: PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,NULL);
287: PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,NULL);
288: PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,NULL);
289: PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,NULL);
290: PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,NULL);
291: PetscOptionsTail();
292: return(0);
293: }
297: static PetscErrorCode SNESView_NEWTONTR(SNES snes,PetscViewer viewer)
298: {
299: SNES_NEWTONTR *tr = (SNES_NEWTONTR*)snes->data;
301: PetscBool iascii;
304: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
305: if (iascii) {
306: PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%g\n",(double)tr->mu,(double)tr->eta,(double)tr->sigma);
307: PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%g\n",(double)tr->delta0,(double)tr->delta1,(double)tr->delta2,(double)tr->delta3);
308: }
309: return(0);
310: }
311: /* ------------------------------------------------------------ */
312: /*MC
313: SNESNEWTONTR - Newton based nonlinear solver that uses a trust region
315: Options Database:
316: + -snes_trtol <tol> Trust region tolerance
317: . -snes_tr_mu <mu>
318: . -snes_tr_eta <eta>
319: . -snes_tr_sigma <sigma>
320: . -snes_tr_delta0 <delta0> - initial size of the trust region is delta0*norm2(x)
321: . -snes_tr_delta1 <delta1>
322: . -snes_tr_delta2 <delta2>
323: - -snes_tr_delta3 <delta3>
325: The basic algorithm is taken from "The Minpack Project", by More',
326: Sorensen, Garbow, Hillstrom, pages 88-111 of "Sources and Development
327: of Mathematical Software", Wayne Cowell, editor.
329: Level: intermediate
331: .seealso: SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESSetTrustRegionTolerance()
333: M*/
336: PETSC_EXTERN PetscErrorCode SNESCreate_NEWTONTR(SNES snes)
337: {
338: SNES_NEWTONTR *neP;
342: snes->ops->setup = SNESSetUp_NEWTONTR;
343: snes->ops->solve = SNESSolve_NEWTONTR;
344: snes->ops->destroy = SNESDestroy_NEWTONTR;
345: snes->ops->setfromoptions = SNESSetFromOptions_NEWTONTR;
346: snes->ops->view = SNESView_NEWTONTR;
347: snes->ops->reset = SNESReset_NEWTONTR;
349: snes->usesksp = PETSC_TRUE;
350: snes->usespc = PETSC_FALSE;
352: PetscNewLog(snes,&neP);
353: snes->data = (void*)neP;
354: neP->mu = 0.25;
355: neP->eta = 0.75;
356: neP->delta = 0.0;
357: neP->delta0 = 0.2;
358: neP->delta1 = 0.3;
359: neP->delta2 = 0.75;
360: neP->delta3 = 2.0;
361: neP->sigma = 0.0001;
362: neP->itflag = PETSC_FALSE;
363: neP->rnorm0 = 0.0;
364: neP->ttol = 0.0;
365: return(0);
366: }