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