Actual source code: snesncg.c

  1: #include <../src/snes/impls/ncg/snesncgimpl.h>
  2: const char *const SNESNCGTypes[] = {"FR", "PRP", "HS", "DY", "CD", "SNESNCGType", "SNES_NCG_", NULL};

  4: static PetscErrorCode SNESReset_NCG(SNES snes)
  5: {
  6:   return 0;
  7: }

  9: /*
 10:   SNESDestroy_NCG - Destroys the private SNES_NCG context that was created with SNESCreate_NCG().

 12:   Input Parameter:
 13: . snes - the SNES context

 15:   Application Interface Routine: SNESDestroy()
 16: */
 17: static PetscErrorCode SNESDestroy_NCG(SNES snes)
 18: {
 19:   PetscObjectComposeFunction((PetscObject)snes, "SNESNCGSetType_C", NULL);
 20:   PetscFree(snes->data);
 21:   return 0;
 22: }

 24: /*
 25:    SNESSetUp_NCG - Sets up the internal data structures for the later use
 26:    of the SNESNCG nonlinear solver.

 28:    Input Parameters:
 29: +  snes - the SNES context
 30: -  x - the solution vector

 32:    Application Interface Routine: SNESSetUp()
 33:  */

 35: static PetscErrorCode SNESSetUp_NCG(SNES snes)
 36: {
 37:   SNESSetWorkVecs(snes, 2);
 39:   if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED;
 40:   return 0;
 41: }

 43: static PetscErrorCode SNESLineSearchApply_NCGLinear(SNESLineSearch linesearch)
 44: {
 45:   PetscScalar alpha, ptAp;
 46:   Vec         X, Y, F, W;
 47:   SNES        snes;
 48:   PetscReal  *fnorm, *xnorm, *ynorm;

 50:   SNESLineSearchGetSNES(linesearch, &snes);
 51:   X     = linesearch->vec_sol;
 52:   W     = linesearch->vec_sol_new;
 53:   F     = linesearch->vec_func;
 54:   Y     = linesearch->vec_update;
 55:   fnorm = &linesearch->fnorm;
 56:   xnorm = &linesearch->xnorm;
 57:   ynorm = &linesearch->ynorm;

 59:   if (!snes->jacobian) SNESSetUpMatrices(snes);

 61:   /*

 63:    The exact step size for unpreconditioned linear CG is just:
 64:    alpha = (r, r) / (p, Ap) = (f, f) / (y, Jy)
 65:    */
 66:   SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
 67:   VecDot(F, F, &alpha);
 68:   MatMult(snes->jacobian, Y, W);
 69:   VecDot(Y, W, &ptAp);
 70:   alpha = alpha / ptAp;
 71:   VecAXPY(X, -alpha, Y);
 72:   SNESComputeFunction(snes, X, F);

 74:   VecNorm(F, NORM_2, fnorm);
 75:   VecNorm(X, NORM_2, xnorm);
 76:   VecNorm(Y, NORM_2, ynorm);
 77:   return 0;
 78: }

 80: /*MC
 81:    SNESLINESEARCHNCGLINEAR - Special line search only for the nonlinear CG solver `SNESNCG`

 83:    This line search uses the length "as if" the problem is linear (that is what is computed by the linear CG method) using the Jacobian of the function.
 84:    alpha = (r, r) / (p, Ap) = (f, f) / (y, Jy) where r (f) is the current residual (function value), p (y) is the current search direction.

 86:    Notes:
 87:    This requires a Jacobian-vector product but does not require the solution of a linear system with the Jacobian

 89:    This is a "odd-ball" line search, we don't know if it is in the literature or used in practice by anyone.

 91:    Level: advanced

 93: .seealso: `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
 94: M*/

 96: PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NCGLinear(SNESLineSearch linesearch)
 97: {
 98:   linesearch->ops->apply          = SNESLineSearchApply_NCGLinear;
 99:   linesearch->ops->destroy        = NULL;
100:   linesearch->ops->setfromoptions = NULL;
101:   linesearch->ops->reset          = NULL;
102:   linesearch->ops->view           = NULL;
103:   linesearch->ops->setup          = NULL;
104:   return 0;
105: }

107: /*
108:   SNESSetFromOptions_NCG - Sets various parameters for the SNESNCG method.

110:   Input Parameter:
111: . snes - the SNES context

113:   Application Interface Routine: SNESSetFromOptions()
114: */
115: static PetscErrorCode SNESSetFromOptions_NCG(SNES snes, PetscOptionItems *PetscOptionsObject)
116: {
117:   SNES_NCG      *ncg     = (SNES_NCG *)snes->data;
118:   PetscBool      debug   = PETSC_FALSE;
119:   SNESNCGType    ncgtype = ncg->type;
120:   SNESLineSearch linesearch;

122:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES NCG options");
123:   PetscOptionsBool("-snes_ncg_monitor", "Monitor the beta values used in the NCG iterations", "SNES", ncg->monitor ? PETSC_TRUE : PETSC_FALSE, &debug, NULL);
124:   if (debug) ncg->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
125:   PetscOptionsEnum("-snes_ncg_type", "NCG Beta type used", "SNESNCGSetType", SNESNCGTypes, (PetscEnum)ncg->type, (PetscEnum *)&ncgtype, NULL);
126:   SNESNCGSetType(snes, ncgtype);
127:   PetscOptionsHeadEnd();
128:   if (!snes->linesearch) {
129:     SNESGetLineSearch(snes, &linesearch);
130:     if (!((PetscObject)linesearch)->type_name) {
131:       if (!snes->npc) {
132:         SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);
133:       } else {
134:         SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
135:       }
136:     }
137:   }
138:   return 0;
139: }

141: /*
142:   SNESView_NCG - Prints info from the SNESNCG data structure.

144:   Input Parameters:
145: + SNES - the SNES context
146: - viewer - visualization context

148:   Application Interface Routine: SNESView()
149: */
150: static PetscErrorCode SNESView_NCG(SNES snes, PetscViewer viewer)
151: {
152:   SNES_NCG *ncg = (SNES_NCG *)snes->data;
153:   PetscBool iascii;

155:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
156:   if (iascii) PetscViewerASCIIPrintf(viewer, "  type: %s\n", SNESNCGTypes[ncg->type]);
157:   return 0;
158: }

160: /*

162:  Assuming F = SNESComputeFunction(X) compute Y^tJ^tF using a simple secant approximation of the jacobian.

164:  This routine is not currently used. I don't know what its intended purpose is.

166:  Note it has a hardwired differencing parameter of 1e-5

168:  */
169: PetscErrorCode SNESNCGComputeYtJtF_Private(SNES snes, Vec X, Vec F, Vec Y, Vec W, Vec G, PetscScalar *ytJtf)
170: {
171:   PetscScalar ftf, ftg, fty, h;

173:   VecDot(F, F, &ftf);
174:   VecDot(F, Y, &fty);
175:   h = 1e-5 * fty / fty;
176:   VecCopy(X, W);
177:   VecAXPY(W, -h, Y); /* this is arbitrary */
178:   SNESComputeFunction(snes, W, G);
179:   VecDot(G, F, &ftg);
180:   *ytJtf = (ftg - ftf) / h;
181:   return 0;
182: }

184: /*@
185:     SNESNCGSetType - Sets the conjugate update type for nonlinear CG `SNESNCG`.

187:     Logically Collective

189:     Input Parameters:
190: +   snes - the iterative context
191: -   btype - update type

193:     Options Database Key:
194: .   -snes_ncg_type <prp,fr,hs,dy,cd> -strategy for selecting algorithm for computing beta

196:     Level: intermediate

198:     `SNESNCGType`s:
199: +   `SNES_NCG_FR` - Fletcher-Reeves update
200: .   `SNES_NCG_PRP` - Polak-Ribiere-Polyak update
201: .   `SNES_NCG_HS` - Hestenes-Steifel update
202: .   `SNES_NCG_DY` - Dai-Yuan update
203: -   `SNES_NCG_CD` - Conjugate Descent update

205:    Notes:
206:    `SNES_NCG_PRP` is the default, and the only one that tolerates generalized search directions.

208:    It is not clear what "generalized search directions" means, does it mean use with a nonlinear preconditioner,
209:    that is using -npc_snes_type <type>, `SNESSetNPC()`, or `SNESGetNPC()`?

211:    Developer Note:
212:    There should be a `SNESNCGSetType()`

214: .seealso: `SNESNCGType`, `SNES_NCG_FR`, `SNES_NCG_PRP`, `SNES_NCG_HS`, `SNES_NCG_DY`, `SNES_NCG_CD`
215: @*/
216: PetscErrorCode SNESNCGSetType(SNES snes, SNESNCGType btype)
217: {
219:   PetscTryMethod(snes, "SNESNCGSetType_C", (SNES, SNESNCGType), (snes, btype));
220:   return 0;
221: }

223: static PetscErrorCode SNESNCGSetType_NCG(SNES snes, SNESNCGType btype)
224: {
225:   SNES_NCG *ncg = (SNES_NCG *)snes->data;

227:   ncg->type = btype;
228:   return 0;
229: }

231: /*
232:   SNESSolve_NCG - Solves a nonlinear system with the Nonlinear Conjugate Gradient method.

234:   Input Parameters:
235: . snes - the SNES context

237:   Output Parameter:
238: . outits - number of iterations until termination

240:   Application Interface Routine: SNESSolve()
241: */
242: static PetscErrorCode SNESSolve_NCG(SNES snes)
243: {
244:   SNES_NCG            *ncg = (SNES_NCG *)snes->data;
245:   Vec                  X, dX, lX, F, dXold;
246:   PetscReal            fnorm, ynorm, xnorm, beta = 0.0;
247:   PetscScalar          dXdotdX, dXolddotdXold, dXdotdXold, lXdotdX, lXdotdXold;
248:   PetscInt             maxits, i;
249:   SNESLineSearchReason lsresult = SNES_LINESEARCH_SUCCEEDED;
250:   SNESLineSearch       linesearch;
251:   SNESConvergedReason  reason;


255:   PetscCitationsRegister(SNESCitation, &SNEScite);
256:   snes->reason = SNES_CONVERGED_ITERATING;

258:   maxits = snes->max_its;        /* maximum number of iterations */
259:   X      = snes->vec_sol;        /* X^n */
260:   dXold  = snes->work[0];        /* The previous iterate of X */
261:   dX     = snes->work[1];        /* the preconditioned direction */
262:   lX     = snes->vec_sol_update; /* the conjugate direction */
263:   F      = snes->vec_func;       /* residual vector */

265:   SNESGetLineSearch(snes, &linesearch);

267:   PetscObjectSAWsTakeAccess((PetscObject)snes);
268:   snes->iter = 0;
269:   snes->norm = 0.;
270:   PetscObjectSAWsGrantAccess((PetscObject)snes);

272:   /* compute the initial function and preconditioned update dX */

274:   if (snes->npc && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
275:     SNESApplyNPC(snes, X, NULL, dX);
276:     SNESGetConvergedReason(snes->npc, &reason);
277:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
278:       snes->reason = SNES_DIVERGED_INNER;
279:       return 0;
280:     }
281:     VecCopy(dX, F);
282:     VecNorm(F, NORM_2, &fnorm);
283:   } else {
284:     if (!snes->vec_func_init_set) {
285:       SNESComputeFunction(snes, X, F);
286:     } else snes->vec_func_init_set = PETSC_FALSE;

288:     /* convergence test */
289:     VecNorm(F, NORM_2, &fnorm);
290:     SNESCheckFunctionNorm(snes, fnorm);
291:     VecCopy(F, dX);
292:   }
293:   if (snes->npc) {
294:     if (snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
295:       SNESApplyNPC(snes, X, F, dX);
296:       SNESGetConvergedReason(snes->npc, &reason);
297:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
298:         snes->reason = SNES_DIVERGED_INNER;
299:         return 0;
300:       }
301:     }
302:   }
303:   VecCopy(dX, lX);
304:   VecDot(dX, dX, &dXdotdX);

306:   PetscObjectSAWsTakeAccess((PetscObject)snes);
307:   snes->norm = fnorm;
308:   PetscObjectSAWsGrantAccess((PetscObject)snes);
309:   SNESLogConvergenceHistory(snes, fnorm, 0);
310:   SNESMonitor(snes, 0, fnorm);

312:   /* test convergence */
313:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
314:   if (snes->reason) return 0;

316:   /* Call general purpose update function */
317:   PetscTryTypeMethod(snes, update, snes->iter);

319:   /* first update -- just use the (preconditioned) residual direction for the initial conjugate direction */

321:   for (i = 1; i < maxits + 1; i++) {
322:     /* some update types require the old update direction or conjugate direction */
323:     if (ncg->type != SNES_NCG_FR) VecCopy(dX, dXold);
324:     SNESLineSearchApply(linesearch, X, F, &fnorm, lX);
325:     SNESLineSearchGetReason(linesearch, &lsresult);
326:     SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);
327:     if (lsresult) {
328:       if (++snes->numFailures >= snes->maxFailures) {
329:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
330:         return 0;
331:       }
332:     }
333:     if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
334:       snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
335:       return 0;
336:     }
337:     /* Monitor convergence */
338:     PetscObjectSAWsTakeAccess((PetscObject)snes);
339:     snes->iter  = i;
340:     snes->norm  = fnorm;
341:     snes->xnorm = xnorm;
342:     snes->ynorm = ynorm;
343:     PetscObjectSAWsGrantAccess((PetscObject)snes);
344:     SNESLogConvergenceHistory(snes, snes->norm, 0);
345:     SNESMonitor(snes, snes->iter, snes->norm);

347:     /* Test for convergence */
348:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
349:     if (snes->reason) return 0;

351:     /* Call general purpose update function */
352:     PetscTryTypeMethod(snes, update, snes->iter);
353:     if (snes->npc) {
354:       if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
355:         SNESApplyNPC(snes, X, NULL, dX);
356:         SNESGetConvergedReason(snes->npc, &reason);
357:         if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
358:           snes->reason = SNES_DIVERGED_INNER;
359:           return 0;
360:         }
361:         VecCopy(dX, F);
362:       } else {
363:         SNESApplyNPC(snes, X, F, dX);
364:         SNESGetConvergedReason(snes->npc, &reason);
365:         if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
366:           snes->reason = SNES_DIVERGED_INNER;
367:           return 0;
368:         }
369:       }
370:     } else {
371:       VecCopy(F, dX);
372:     }

374:     /* compute the conjugate direction lX = dX + beta*lX with beta = ((dX, dX) / (dX_old, dX_old) (Fletcher-Reeves update)*/
375:     switch (ncg->type) {
376:     case SNES_NCG_FR: /* Fletcher-Reeves */
377:       dXolddotdXold = dXdotdX;
378:       VecDot(dX, dX, &dXdotdX);
379:       beta = PetscRealPart(dXdotdX / dXolddotdXold);
380:       break;
381:     case SNES_NCG_PRP: /* Polak-Ribiere-Poylak */
382:       dXolddotdXold = dXdotdX;
383:       VecDotBegin(dX, dX, &dXdotdXold);
384:       VecDotBegin(dXold, dX, &dXdotdXold);
385:       VecDotEnd(dX, dX, &dXdotdX);
386:       VecDotEnd(dXold, dX, &dXdotdXold);
387:       beta = PetscRealPart(((dXdotdX - dXdotdXold) / dXolddotdXold));
388:       if (beta < 0.0) beta = 0.0; /* restart */
389:       break;
390:     case SNES_NCG_HS: /* Hestenes-Stiefel */
391:       VecDotBegin(dX, dX, &dXdotdX);
392:       VecDotBegin(dX, dXold, &dXdotdXold);
393:       VecDotBegin(lX, dX, &lXdotdX);
394:       VecDotBegin(lX, dXold, &lXdotdXold);
395:       VecDotEnd(dX, dX, &dXdotdX);
396:       VecDotEnd(dX, dXold, &dXdotdXold);
397:       VecDotEnd(lX, dX, &lXdotdX);
398:       VecDotEnd(lX, dXold, &lXdotdXold);
399:       beta = PetscRealPart((dXdotdX - dXdotdXold) / (lXdotdX - lXdotdXold));
400:       break;
401:     case SNES_NCG_DY: /* Dai-Yuan */
402:       VecDotBegin(dX, dX, &dXdotdX);
403:       VecDotBegin(lX, dX, &lXdotdX);
404:       VecDotBegin(lX, dXold, &lXdotdXold);
405:       VecDotEnd(dX, dX, &dXdotdX);
406:       VecDotEnd(lX, dX, &lXdotdX);
407:       VecDotEnd(lX, dXold, &lXdotdXold);
408:       beta = PetscRealPart(dXdotdX / (lXdotdXold - lXdotdX));
409:       break;
410:     case SNES_NCG_CD: /* Conjugate Descent */
411:       VecDotBegin(dX, dX, &dXdotdX);
412:       VecDotBegin(lX, dXold, &lXdotdXold);
413:       VecDotEnd(dX, dX, &dXdotdX);
414:       VecDotEnd(lX, dXold, &lXdotdXold);
415:       beta = PetscRealPart(dXdotdX / lXdotdXold);
416:       break;
417:     }
418:     if (ncg->monitor) PetscViewerASCIIPrintf(ncg->monitor, "beta = %e\n", (double)beta);
419:     VecAYPX(lX, beta, dX);
420:   }
421:   PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", maxits);
422:   if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
423:   return 0;
424: }

426: /*MC
427:   SNESNCG - Nonlinear Conjugate-Gradient method for the solution of nonlinear systems.

429:   Level: beginner

431:   Options Database Keys:
432: +   -snes_ncg_type <fr, prp, dy, hs, cd> - Choice of conjugate-gradient update parameter, default is prp.
433: .   -snes_linesearch_type <cp,l2,basic> - Line search type.
434: -   -snes_ncg_monitor - Print the beta values used in the ncg iteration, .

436:    Notes:
437:    This solves the nonlinear system of equations F(x) = 0 using the nonlinear generalization of the conjugate
438:    gradient method.  This may be used with a nonlinear preconditioner used to pick the new search directions, but otherwise
439:    chooses the initial search direction as F(x) for the initial guess x.

441:    Only supports left non-linear preconditioning.

443:    Default line search is `SNESLINESEARCHCP`, unless a nonlinear preconditioner is used with -npc_snes_type <type>, `SNESSetNPC()`, or `SNESGetNPC()` then
444:    `SNESLINESEARCHL2` is used. Also supports the special purpose line search `SNESLINESEARCHNCGLINEAR`

446:    References:
447: .  * -  Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu,"Composing Scalable Nonlinear Algebraic Solvers",
448:    SIAM Review, 57(4), 2015

450: .seealso: `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNGMRES`, `SNESQN`, `SNESLINESEARCHNCGLINEAR`, `SNESNCGSetType()`, `SNESLineSearchSetType()`
451: M*/
452: PETSC_EXTERN PetscErrorCode SNESCreate_NCG(SNES snes)
453: {
454:   SNES_NCG *neP;

456:   snes->ops->destroy        = SNESDestroy_NCG;
457:   snes->ops->setup          = SNESSetUp_NCG;
458:   snes->ops->setfromoptions = SNESSetFromOptions_NCG;
459:   snes->ops->view           = SNESView_NCG;
460:   snes->ops->solve          = SNESSolve_NCG;
461:   snes->ops->reset          = SNESReset_NCG;

463:   snes->usesksp = PETSC_FALSE;
464:   snes->usesnpc = PETSC_TRUE;
465:   snes->npcside = PC_LEFT;

467:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

469:   if (!snes->tolerancesset) {
470:     snes->max_funcs = 30000;
471:     snes->max_its   = 10000;
472:     snes->stol      = 1e-20;
473:   }

475:   PetscNew(&neP);
476:   snes->data   = (void *)neP;
477:   neP->monitor = NULL;
478:   neP->type    = SNES_NCG_PRP;
479:   PetscObjectComposeFunction((PetscObject)snes, "SNESNCGSetType_C", SNESNCGSetType_NCG);
480:   return 0;
481: }