Actual source code: tr.c
1: /*$Id: tr.c,v 1.128 2001/08/07 03:04:11 balay Exp $*/
2:
3: #include src/snes/impls/tr/tr.h
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__
11: int SNES_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_TR *neP = (SNES_TR*)snes->data;
16: Vec x;
17: PetscReal nrm;
18: int ierr;
21: if (snes->ksp_ewconv) {
22: if (!kctx) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Eisenstat-Walker onvergence context not created");
23: if (!n) {SNES_KSP_EW_ComputeRelativeTolerance_Private(snes,ksp);}
24: kctx->lresid_last = rnorm;
25: }
26: KSPDefaultConverged(ksp,n,rnorm,reason,ctx);
27: if (*reason) {
28: PetscLogInfo(snes,"SNES_TR_KSPConverged_Private: regular convergence test KSP iterations=%d, rnorm=%gn",n,rnorm);
29: }
31: /* Determine norm of solution */
32: KSPBuildSolution(ksp,0,&x);
33: VecNorm(x,NORM_2,&nrm);
34: if (nrm >= neP->delta) {
35: PetscLogInfo(snes,"SNES_TR_KSPConverged_Private: KSP iterations=%d, rnorm=%gn",n,rnorm);
36: PetscLogInfo(snes,"SNES_TR_KSPConverged_Private: Ending linear iteration early, delta=%g, length=%gn",neP->delta,nrm);
37: *reason = KSP_CONVERGED_STEP_LENGTH;
38: }
39: return(0);
40: }
42: /*
43: SNESSolve_TR - Implements Newton's Method with a very simple trust
44: region approach for solving systems of nonlinear equations.
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.
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__
56: static int SNESSolve_TR(SNES snes,int *its)
57: {
58: SNES_TR *neP = (SNES_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,nrm,ynorm,norm1;
63: PetscScalar mone = -1.0,cnorm;
64: KSP ksp;
65: SLES sles;
66: SNESConvergedReason reason;
67: PetscTruth conv;
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];
77: PetscObjectTakeAccess(snes);
78: snes->iter = 0;
79: PetscObjectGrantAccess(snes);
80: VecNorm(X,NORM_2,&xnorm); /* xnorm = || X || */
82: SNESComputeFunction(snes,X,F); /* F(X) */
83: VecNorm(F,NORM_2,&fnorm); /* fnorm <- || F || */
84: PetscObjectTakeAccess(snes);
85: snes->norm = fnorm;
86: PetscObjectGrantAccess(snes);
87: delta = neP->delta0*fnorm;
88: neP->delta = delta;
89: SNESLogConvHistory(snes,fnorm,0);
90: SNESMonitor(snes,0,fnorm);
92: if (fnorm < snes->atol) {*its = 0; snes->reason = SNES_CONVERGED_FNORM_ABS; return(0);}
94: /* set parameter for default relative tolerance convergence test */
95: snes->ttol = fnorm*snes->rtol;
97: /* Set the stopping criteria to use the More' trick. */
98: PetscOptionsHasName(PETSC_NULL,"-snes_tr_ksp_regular_convergence_test",&conv);
99: if (!conv) {
100: SNESGetSLES(snes,&sles);
101: SLESGetKSP(sles,&ksp);
102: KSPSetConvergenceTest(ksp,SNES_TR_KSPConverged_Private,(void *)snes);
103: PetscLogInfo(snes,"SNESSolve_TR: Using Krylov convergence test SNES_TR_KSPConverged_Privaten");
104: }
105:
106: for (i=0; i<maxits; i++) {
107: SNESComputeJacobian(snes,X,&snes->jacobian,&snes->jacobian_pre,&flg);
108: SLESSetOperators(snes->sles,snes->jacobian,snes->jacobian_pre,flg);
110: /* Solve J Y = F, where J is Jacobian matrix */
111: SLESSolve(snes->sles,F,Ytmp,&lits);
112: snes->linear_its += lits;
113: PetscLogInfo(snes,"SNESSolve_TR: iter=%d, linear solve iterations=%dn",snes->iter,lits);
114: VecNorm(Ytmp,NORM_2,&nrm);
115: norm1 = nrm;
116: while(1) {
117: VecCopy(Ytmp,Y);
118: nrm = norm1;
120: /* Scale Y if need be and predict new value of F norm */
121: if (nrm >= delta) {
122: nrm = delta/nrm;
123: gpnorm = (1.0 - nrm)*fnorm;
124: cnorm = nrm;
125: PetscLogInfo(snes,"SNESSolve_TR: Scaling direction by %gn",nrm);
126: VecScale(&cnorm,Y);
127: nrm = gpnorm;
128: ynorm = delta;
129: } else {
130: gpnorm = 0.0;
131: PetscLogInfo(snes,"SNESSolve_TR: Direction is in Trust Regionn");
132: ynorm = nrm;
133: }
134: VecAYPX(&mone,X,Y); /* Y <- X - Y */
135: VecCopy(X,snes->vec_sol_update_always);
136: SNESComputeFunction(snes,Y,G); /* F(X) */
137: VecNorm(G,NORM_2,&gnorm); /* gnorm <- || g || */
138: if (fnorm == gpnorm) rho = 0.0;
139: else rho = (fnorm*fnorm - gnorm*gnorm)/(fnorm*fnorm - gpnorm*gpnorm);
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_TR: fnorm=%g, gnorm=%g, ynorm=%gn",fnorm,gnorm,ynorm);
146: PetscLogInfo(snes,"SNESSolve_TR: gpred=%g, rho=%g, delta=%gn",gpnorm,rho,delta);
147: neP->delta = delta;
148: if (rho > neP->sigma) break;
149: PetscLogInfo(snes,"SNESSolve_TR: Trying again in smaller regionn");
150: /* check to see if progress is hopeless */
151: neP->itflag = 0;
152: (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);
153: if (reason) {
154: /* We're not progressing, so return with the current iterate */
155: SNESMonitor(snes,i+1,fnorm);
156: breakout = 1;
157: break;
158: }
159: snes->numFailures++;
160: }
161: if (!breakout) {
162: fnorm = gnorm;
163: PetscObjectTakeAccess(snes);
164: snes->iter = i+1;
165: snes->norm = fnorm;
166: PetscObjectGrantAccess(snes);
167: TMP = F; F = G; snes->vec_func_always = F; G = TMP;
168: TMP = X; X = Y; snes->vec_sol_always = X; Y = TMP;
169: VecNorm(X,NORM_2,&xnorm); /* xnorm = || X || */
170: SNESLogConvHistory(snes,fnorm,lits);
171: SNESMonitor(snes,i+1,fnorm);
173: /* Test for convergence */
174: neP->itflag = 1;
175: (*snes->converged)(snes,xnorm,ynorm,fnorm,&reason,snes->cnvP);
176: if (reason) {
177: break;
178: }
179: } else {
180: break;
181: }
182: }
183: /* Verify solution is in corect location */
184: if (X != snes->vec_sol) {
185: VecCopy(X,snes->vec_sol);
186: }
187: if (F != snes->vec_func) {
188: VecCopy(F,snes->vec_func);
189: }
190: snes->vec_sol_always = snes->vec_sol;
191: snes->vec_func_always = snes->vec_func;
192: if (i == maxits) {
193: PetscLogInfo(snes,"SNESSolve_TR: Maximum number of iterations has been reached: %dn",maxits);
194: i--;
195: reason = SNES_DIVERGED_MAX_IT;
196: }
197: *its = i+1;
198: PetscObjectTakeAccess(snes);
199: snes->reason = reason;
200: PetscObjectGrantAccess(snes);
201: return(0);
202: }
203: /*------------------------------------------------------------*/
204: #undef __FUNCT__
206: static int SNESSetUp_TR(SNES snes)
207: {
211: snes->nwork = 4;
212: VecDuplicateVecs(snes->vec_sol,snes->nwork,&snes->work);
213: PetscLogObjectParents(snes,snes->nwork,snes->work);
214: snes->vec_sol_update_always = snes->work[3];
215: return(0);
216: }
217: /*------------------------------------------------------------*/
218: #undef __FUNCT__
220: static int SNESDestroy_TR(SNES snes)
221: {
222: int ierr;
225: if (snes->nwork) {
226: VecDestroyVecs(snes->work,snes->nwork);
227: }
228: PetscFree(snes->data);
229: return(0);
230: }
231: /*------------------------------------------------------------*/
233: #undef __FUNCT__
235: static int SNESSetFromOptions_TR(SNES snes)
236: {
237: SNES_TR *ctx = (SNES_TR *)snes->data;
238: int ierr;
241: PetscOptionsHead("SNES trust region options for nonlinear equations");
242: PetscOptionsReal("-snes_trtol","Trust region tolerance","SNESSetTrustRegionTolerance",snes->deltatol,&snes->deltatol,0);
243: PetscOptionsReal("-snes_tr_mu","mu","None",ctx->mu,&ctx->mu,0);
244: PetscOptionsReal("-snes_tr_eta","eta","None",ctx->eta,&ctx->eta,0);
245: PetscOptionsReal("-snes_tr_sigma","sigma","None",ctx->sigma,&ctx->sigma,0);
246: PetscOptionsReal("-snes_tr_delta0","delta0","None",ctx->delta0,&ctx->delta0,0);
247: PetscOptionsReal("-snes_tr_delta1","delta1","None",ctx->delta1,&ctx->delta1,0);
248: PetscOptionsReal("-snes_tr_delta2","delta2","None",ctx->delta2,&ctx->delta2,0);
249: PetscOptionsReal("-snes_tr_delta3","delta3","None",ctx->delta3,&ctx->delta3,0);
250: PetscOptionsTail();
251: return(0);
252: }
254: #undef __FUNCT__
256: static int SNESView_TR(SNES snes,PetscViewer viewer)
257: {
258: SNES_TR *tr = (SNES_TR *)snes->data;
259: int ierr;
260: PetscTruth isascii;
263: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
264: if (isascii) {
265: PetscViewerASCIIPrintf(viewer," mu=%g, eta=%g, sigma=%gn",tr->mu,tr->eta,tr->sigma);
266: PetscViewerASCIIPrintf(viewer," delta0=%g, delta1=%g, delta2=%g, delta3=%gn",tr->delta0,tr->delta1,tr->delta2,tr->delta3);
267: } else {
268: SETERRQ1(1,"Viewer type %s not supported for SNES EQ TR",((PetscObject)viewer)->type_name);
269: }
270: return(0);
271: }
273: /* ---------------------------------------------------------------- */
274: #undef __FUNCT__
276: /*@C
277: SNESConverged_TR - Monitors the convergence of the trust region
278: method SNESTR for solving systems of nonlinear equations (default).
280: Collective on SNES
282: Input Parameters:
283: + snes - the SNES context
284: . xnorm - 2-norm of current iterate
285: . pnorm - 2-norm of current step
286: . fnorm - 2-norm of function
287: - dummy - unused context
289: Output Parameter:
290: . reason - one of
291: $ SNES_CONVERGED_FNORM_ABS - (fnorm < atol),
292: $ SNES_CONVERGED_PNORM_RELATIVE - (pnorm < xtol*xnorm),
293: $ SNES_CONVERGED_FNORM_RELATIVE - (fnorm < rtol*fnorm0),
294: $ SNES_DIVERGED_FUNCTION_COUNT - (nfct > maxf),
295: $ SNES_DIVERGED_FNORM_NAN - (fnorm == NaN),
296: $ SNES_CONVERGED_TR_DELTA - (delta < xnorm*deltatol),
297: $ SNES_CONVERGED_ITERATING - (otherwise)
299: where
300: + delta - trust region paramenter
301: . deltatol - trust region size tolerance,
302: set with SNESSetTrustRegionTolerance()
303: . maxf - maximum number of function evaluations,
304: set with SNESSetTolerances()
305: . nfct - number of function evaluations,
306: . atol - absolute function norm tolerance,
307: set with SNESSetTolerances()
308: - xtol - relative function norm tolerance,
309: set with SNESSetTolerances()
311: Level: intermediate
313: .keywords: SNES, nonlinear, default, converged, convergence
315: .seealso: SNESSetConvergenceTest(), SNESEisenstatWalkerConverged()
316: @*/
317: int SNESConverged_TR(SNES snes,PetscReal xnorm,PetscReal pnorm,PetscReal fnorm,SNESConvergedReason *reason,void *dummy)
318: {
319: SNES_TR *neP = (SNES_TR *)snes->data;
320: int ierr;
323: if (fnorm != fnorm) {
324: PetscLogInfo(snes,"SNESConverged_TR:Failed to converged, function norm is NaNn");
325: *reason = SNES_DIVERGED_FNORM_NAN;
326: } else if (neP->delta < xnorm * snes->deltatol) {
327: PetscLogInfo(snes,"SNESConverged_TR: Converged due to trust region param %g<%g*%gn",neP->delta,xnorm,snes->deltatol);
328: *reason = SNES_CONVERGED_TR_DELTA;
329: } else if (neP->itflag) {
330: SNESConverged_LS(snes,xnorm,pnorm,fnorm,reason,dummy);
331: } else if (snes->nfuncs > snes->max_funcs) {
332: PetscLogInfo(snes,"SNESConverged_TR: Exceeded maximum number of function evaluations: %d > %dn",snes->nfuncs,snes->max_funcs);
333: *reason = SNES_DIVERGED_FUNCTION_COUNT;
334: } else {
335: *reason = SNES_CONVERGED_ITERATING;
336: }
337: return(0);
338: }
339: /* ------------------------------------------------------------ */
340: EXTERN_C_BEGIN
341: #undef __FUNCT__
343: int SNESCreate_TR(SNES snes)
344: {
345: SNES_TR *neP;
346: int ierr;
349: snes->setup = SNESSetUp_TR;
350: snes->solve = SNESSolve_TR;
351: snes->destroy = SNESDestroy_TR;
352: snes->converged = SNESConverged_TR;
353: snes->setfromoptions = SNESSetFromOptions_TR;
354: snes->view = SNESView_TR;
355: snes->nwork = 0;
356:
357: ierr = PetscNew(SNES_TR,&neP);
358: PetscLogObjectMemory(snes,sizeof(SNES_TR));
359: snes->data = (void*)neP;
360: neP->mu = 0.25;
361: neP->eta = 0.75;
362: neP->delta = 0.0;
363: neP->delta0 = 0.2;
364: neP->delta1 = 0.3;
365: neP->delta2 = 0.75;
366: neP->delta3 = 2.0;
367: neP->sigma = 0.0001;
368: neP->itflag = 0;
369: neP->rnorm0 = 0;
370: neP->ttol = 0;
371: return(0);
372: }
373: EXTERN_C_END