Actual source code: snesncg.c
petsc-3.3-p7 2013-05-11
1: #include <../src/snes/impls/ncg/snesncgimpl.h> /*I "petscsnes.h" I*/
2: const char *SNESNCGTypes[] = {"FR","PRP","HS","DY","CD","SNESNCGType","SNES_NCG_",0};
6: PetscErrorCode SNESReset_NCG(SNES snes)
7: {
10: return(0);
11: }
13: #define SNESLINESEARCHNCGLINEAR "linear"
15: /*
16: SNESDestroy_NCG - Destroys the private SNES_NCG context that was created with SNESCreate_NCG().
18: Input Parameter:
19: . snes - the SNES context
21: Application Interface Routine: SNESDestroy()
22: */
25: PetscErrorCode SNESDestroy_NCG(SNES snes)
26: {
27: PetscErrorCode ierr;
30: PetscFree(snes->data);
31: return(0);
32: }
34: /*
35: SNESSetUp_NCG - Sets up the internal data structures for the later use
36: of the SNESNCG nonlinear solver.
38: Input Parameters:
39: + snes - the SNES context
40: - x - the solution vector
42: Application Interface Routine: SNESSetUp()
43: */
45: EXTERN_C_BEGIN
46: extern PetscErrorCode SNESLineSearchCreate_NCGLinear(SNESLineSearch);
47: EXTERN_C_END
51: PetscErrorCode SNESSetUp_NCG(SNES snes)
52: {
56: SNESDefaultGetWork(snes,2);
57: SNESSetUpMatrices(snes);
58: SNESLineSearchRegisterDynamic(SNESLINESEARCHNCGLINEAR, PETSC_NULL,"SNESLineSearchCreate_NCGLinear", SNESLineSearchCreate_NCGLinear);
60: return(0);
61: }
62: /*
63: SNESSetFromOptions_NCG - Sets various parameters for the SNESNCG method.
65: Input Parameter:
66: . snes - the SNES context
68: Application Interface Routine: SNESSetFromOptions()
69: */
72: static PetscErrorCode SNESSetFromOptions_NCG(SNES snes)
73: {
74: SNES_NCG *ncg = (SNES_NCG *)snes->data;
75: PetscErrorCode ierr;
76: PetscBool debug;
77: SNESLineSearch linesearch;
78: SNESNCGType ncgtype;
81: PetscOptionsHead("SNES NCG options");
82: PetscOptionsBool("-snes_ncg_monitor","Monitor NCG iterations","SNES",ncg->monitor ? PETSC_TRUE: PETSC_FALSE, &debug, PETSC_NULL);
83: PetscOptionsEnum("-snes_ncg_type","NCG Beta type used","SNESNCGSetType",SNESNCGTypes,(PetscEnum)ncg->type,(PetscEnum*)&ncgtype,PETSC_NULL);
84: SNESNCGSetType(snes, ncgtype);
85: if (debug) {
86: ncg->monitor = PETSC_VIEWER_STDOUT_(((PetscObject)snes)->comm);
87: }
88: PetscOptionsTail();
89: if (!snes->linesearch) {
90: SNESGetSNESLineSearch(snes, &linesearch);
91: if (!snes->pc) {
92: SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);
93: } else {
94: SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
95: }
96: }
97: return(0);
98: }
100: /*
101: SNESView_NCG - Prints info from the SNESNCG data structure.
103: Input Parameters:
104: + SNES - the SNES context
105: - viewer - visualization context
107: Application Interface Routine: SNESView()
108: */
111: static PetscErrorCode SNESView_NCG(SNES snes, PetscViewer viewer)
112: {
113: PetscBool iascii;
114: PetscErrorCode ierr;
117: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
118: if (iascii) {
119: }
120: return(0);
121: }
125: PetscErrorCode SNESLineSearchApply_NCGLinear(SNESLineSearch linesearch)
126: {
127: PetscScalar alpha, ptAp;
128: Vec X, Y, F, W;
129: SNES snes;
130: PetscErrorCode ierr;
131: PetscReal *fnorm, *xnorm, *ynorm;
132: MatStructure flg = DIFFERENT_NONZERO_PATTERN;
136: SNESLineSearchGetSNES(linesearch, &snes);
137: X = linesearch->vec_sol;
138: W = linesearch->vec_sol_new;
139: F = linesearch->vec_func;
140: Y = linesearch->vec_update;
141: fnorm = &linesearch->fnorm;
142: xnorm = &linesearch->xnorm;
143: ynorm = &linesearch->ynorm;
145: /*
147: The exact step size for unpreconditioned linear CG is just:
148: alpha = (r, r) / (p, Ap) = (f, f) / (y, Jy)
149: */
150: SNESComputeJacobian(snes, X, &snes->jacobian, &snes->jacobian_pre, &flg);
151: VecDot(F, F, &alpha);
152: MatMult(snes->jacobian, Y, W);
153: VecDot(Y, W, &ptAp);
154: alpha = alpha / ptAp;
155: PetscPrintf(((PetscObject)snes)->comm, "alpha: %G\n", PetscRealPart(alpha));
156: VecAXPY(X, alpha, Y);
157: SNESComputeFunction(snes, X, F);
159: VecNorm(F, NORM_2, fnorm);
160: VecNorm(X, NORM_2, xnorm);
161: VecNorm(Y, NORM_2, ynorm);
163: return(0);
164: }
166: EXTERN_C_BEGIN
170: PetscErrorCode SNESLineSearchCreate_NCGLinear(SNESLineSearch linesearch)
171: {
173: linesearch->ops->apply = SNESLineSearchApply_NCGLinear;
174: linesearch->ops->destroy = PETSC_NULL;
175: linesearch->ops->setfromoptions = PETSC_NULL;
176: linesearch->ops->reset = PETSC_NULL;
177: linesearch->ops->view = PETSC_NULL;
178: linesearch->ops->setup = PETSC_NULL;
179: return(0);
180: }
181: EXTERN_C_END
185: /*
187: Assuming F = SNESComputeFunction(X) compute Y^tJ^tF using a simple secant approximation of the jacobian.
189: */
190: PetscErrorCode SNESNCGComputeYtJtF_Private(SNES snes, Vec X, Vec F, Vec Y, Vec W, Vec G, PetscScalar * ytJtf) {
192: PetscScalar ftf, ftg, fty, h;
194: VecDot(F, F, &ftf);
195: VecDot(F, Y, &fty);
196: h = 1e-5*fty / fty;
197: VecCopy(X, W);
198: VecAXPY(W, -h, Y); /* this is arbitrary */
199: SNESComputeFunction(snes, W, G);
200: VecDot(G, F, &ftg);
201: *ytJtf = (ftg - ftf) / h;
202: return(0);
203: }
207: /*@
208: SNESNCGSetType - Sets the conjugate update type for SNESNCG.
210: Logically Collective on SNES
212: Input Parameters:
213: + snes - the iterative context
214: - btype - update type
216: Options Database:
217: . -snes_ncg_type<prp,fr,hs,dy,cd>
219: Level: intermediate
221: SNESNCGSelectTypes:
222: + SNES_NCG_FR - Fletcher-Reeves update
223: . SNES_NCG_PRP - Polak-Ribiere-Polyak update
224: . SNES_NCG_HS - Hestenes-Steifel update
225: . SNES_NCG_DY - Dai-Yuan update
226: - SNES_NCG_CD - Conjugate Descent update
228: Notes:
229: PRP is the default, and the only one that tolerates generalized search directions.
231: .keywords: SNES, SNESNCG, selection, type, set
232: @*/
233: PetscErrorCode SNESNCGSetType(SNES snes, SNESNCGType btype) {
237: PetscTryMethod(snes,"SNESNCGSetType_C",(SNES,SNESNCGType),(snes,btype));
238: return(0);
239: }
242: EXTERN_C_BEGIN
245: PetscErrorCode SNESNCGSetType_NCG(SNES snes, SNESNCGType btype) {
246: SNES_NCG *ncg = (SNES_NCG *)snes->data;
248: ncg->type = btype;
249: return(0);
250: }
251: EXTERN_C_END
253: /*
254: SNESSolve_NCG - Solves a nonlinear system with the Nonlinear Conjugate Gradient method.
256: Input Parameters:
257: . snes - the SNES context
259: Output Parameter:
260: . outits - number of iterations until termination
262: Application Interface Routine: SNESSolve()
263: */
266: PetscErrorCode SNESSolve_NCG(SNES snes)
267: {
268: SNES_NCG *ncg = (SNES_NCG *)snes->data;
269: Vec X, dX, lX, F, B, Fold;
270: PetscReal fnorm, ynorm, xnorm, beta = 0.0;
271: PetscScalar dXdotF, dXolddotFold, dXdotFold, lXdotF, lXdotFold;
272: PetscInt maxits, i;
273: PetscErrorCode ierr;
274: SNESConvergedReason reason;
275: PetscBool lsSuccess = PETSC_TRUE;
276: SNESLineSearch linesearch;
279: snes->reason = SNES_CONVERGED_ITERATING;
281: maxits = snes->max_its; /* maximum number of iterations */
282: X = snes->vec_sol; /* X^n */
283: Fold = snes->work[0]; /* The previous iterate of X */
284: dX = snes->work[1]; /* the preconditioned direction */
285: lX = snes->vec_sol_update; /* the conjugate direction */
286: F = snes->vec_func; /* residual vector */
287: B = snes->vec_rhs; /* the right hand side */
289: SNESGetSNESLineSearch(snes, &linesearch);
291: PetscObjectTakeAccess(snes);
292: snes->iter = 0;
293: snes->norm = 0.;
294: PetscObjectGrantAccess(snes);
296: /* compute the initial function and preconditioned update dX */
297: if (!snes->vec_func_init_set) {
298: SNESComputeFunction(snes,X,F);
299: if (snes->domainerror) {
300: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
301: return(0);
302: }
303: } else {
304: snes->vec_func_init_set = PETSC_FALSE;
305: }
306: if (!snes->norm_init_set) {
307: /* convergence test */
308: VecNorm(F, NORM_2, &fnorm); /* fnorm <- ||F|| */
309: if (PetscIsInfOrNanReal(fnorm)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FP,"Infinite or not-a-number generated in norm");
310: } else {
311: fnorm = snes->norm_init;
312: snes->norm_init_set = PETSC_FALSE;
313: }
314: PetscObjectTakeAccess(snes);
315: snes->norm = fnorm;
316: PetscObjectGrantAccess(snes);
317: SNESLogConvHistory(snes,fnorm,0);
318: SNESMonitor(snes,0,fnorm);
320: /* set parameter for default relative tolerance convergence test */
321: snes->ttol = fnorm*snes->rtol;
322: /* test convergence */
323: (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
324: if (snes->reason) return(0);
326: /* Call general purpose update function */
327: if (snes->ops->update) {
328: (*snes->ops->update)(snes, snes->iter);
329: }
331: /* first update -- just use the (preconditioned) residual direction for the initial conjugate direction */
333: if (snes->pc) {
334: VecCopy(X, dX);
335: SNESSetInitialFunction(snes->pc, F);
336: SNESSetInitialFunctionNorm(snes->pc, fnorm);
337: SNESSolve(snes->pc, B, dX);
338: SNESGetConvergedReason(snes->pc,&reason);
339: if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) {
340: snes->reason = SNES_DIVERGED_INNER;
341: return(0);
342: }
343: VecAYPX(dX,-1.0,X);
344: } else {
345: VecCopy(F, dX);
346: }
347: VecCopy(dX, lX);
348: VecDot(F, dX, &dXdotF);
349: /*
350: } else {
351: SNESNCGComputeYtJtF_Private(snes, X, F, dX, W, G, &dXdotF);
352: }
353: */
354: for(i = 1; i < maxits + 1; i++) {
355: lsSuccess = PETSC_TRUE;
356: /* some update types require the old update direction or conjugate direction */
357: if (ncg->type != SNES_NCG_FR) {
358: VecCopy(F, Fold);
359: }
360: SNESLineSearchApply(linesearch, X, F, &fnorm, lX);
361: SNESLineSearchGetSuccess(linesearch, &lsSuccess);
362: if (!lsSuccess) {
363: if (++snes->numFailures >= snes->maxFailures) {
364: snes->reason = SNES_DIVERGED_LINE_SEARCH;
365: return(0);
366: }
367: }
368: if (snes->nfuncs >= snes->max_funcs) {
369: snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
370: return(0);
371: }
372: if (snes->domainerror) {
373: snes->reason = SNES_DIVERGED_FUNCTION_DOMAIN;
374: return(0);
375: }
376: SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm);
377: /* Monitor convergence */
378: PetscObjectTakeAccess(snes);
379: snes->iter = i;
380: snes->norm = fnorm;
381: PetscObjectGrantAccess(snes);
382: SNESLogConvHistory(snes,snes->norm,0);
383: SNESMonitor(snes,snes->iter,snes->norm);
385: /* Test for convergence */
386: (*snes->ops->converged)(snes,snes->iter,xnorm,ynorm,fnorm,&snes->reason,snes->cnvP);
387: if (snes->reason) return(0);
389: /* Call general purpose update function */
390: if (snes->ops->update) {
391: (*snes->ops->update)(snes, snes->iter);
392: }
393: if (snes->pc) {
394: VecCopy(X,dX);
395: SNESSetInitialFunction(snes->pc, F);
396: SNESSetInitialFunctionNorm(snes->pc, fnorm);
397: SNESSolve(snes->pc, B, dX);
398: SNESGetConvergedReason(snes->pc,&reason);
399: if (reason < 0 && (reason != SNES_DIVERGED_MAX_IT)) {
400: snes->reason = SNES_DIVERGED_INNER;
401: return(0);
402: }
403: VecAYPX(dX,-1.0,X);
404: } else {
405: VecCopy(F, dX);
406: }
408: /* compute the conjugate direction lX = dX + beta*lX with beta = ((dX, dX) / (dX_old, dX_old) (Fletcher-Reeves update)*/
409: switch(ncg->type) {
410: case SNES_NCG_FR: /* Fletcher-Reeves */
411: dXolddotFold = dXdotF;
412: VecDot(dX, dX, &dXdotF);
413: beta = PetscRealPart(dXdotF / dXolddotFold);
414: break;
415: case SNES_NCG_PRP: /* Polak-Ribiere-Poylak */
416: dXolddotFold = dXdotF;
417: VecDotBegin(F, dX, &dXdotF);
418: VecDotBegin(Fold, dX, &dXdotFold);
419: VecDotEnd(F, dX, &dXdotF);
420: VecDotEnd(Fold, dX, &dXdotFold);
421: beta = PetscRealPart(((dXdotF - dXdotFold) / dXolddotFold));
422: if (beta < 0.0) beta = 0.0; /* restart */
423: break;
424: case SNES_NCG_HS: /* Hestenes-Stiefel */
425: VecDotBegin(dX, F, &dXdotF);
426: VecDotBegin(dX, Fold, &dXdotFold);
427: VecDotBegin(lX, F, &lXdotF);
428: VecDotBegin(lX, Fold, &lXdotFold);
429: VecDotEnd(dX, F, &dXdotF);
430: VecDotEnd(dX, Fold, &dXdotFold);
431: VecDotEnd(lX, F, &lXdotF);
432: VecDotEnd(lX, Fold, &lXdotFold);
433: beta = PetscRealPart((dXdotF - dXdotFold) / (lXdotF - lXdotFold));
434: break;
435: case SNES_NCG_DY: /* Dai-Yuan */
436: VecDotBegin(dX, F, &dXdotF);
437: VecDotBegin(lX, F, &lXdotF);
438: VecDotBegin(lX, Fold, &lXdotFold);
439: VecDotEnd(dX, F, &dXdotF);
440: VecDotEnd(lX, F, &lXdotF);
441: VecDotEnd(lX, Fold, &lXdotFold);
442: beta = PetscRealPart(dXdotF / (lXdotFold - lXdotF));
443: break;
444: case SNES_NCG_CD: /* Conjugate Descent */
445: VecDotBegin(dX, F, &dXdotF);
446: VecDotBegin(lX, Fold, &lXdotFold);
447: VecDotEnd(dX, F, &dXdotF);
448: VecDotEnd(lX, Fold, &lXdotFold);
449: beta = PetscRealPart(dXdotF / lXdotFold);
450: break;
451: }
452: if (ncg->monitor) {
453: PetscViewerASCIIPrintf(ncg->monitor, "beta = %e\n", beta);
454: }
455: VecAYPX(lX, beta, dX);
456: }
457: PetscInfo1(snes, "Maximum number of iterations has been reached: %D\n", maxits);
458: if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
459: return(0);
460: }
464: /*MC
465: SNESNCG - Nonlinear Conjugate-Gradient method for the solution of general nonlinear systems.
467: Level: beginner
469: Options Database:
470: + -snes_ncg_type <fr, prp, dy, hs, cd> - Choice of conjugate-gradient update parameter.
471: . -snes_linesearch_type <cp,l2,basic> - Line search type.
472: - -snes_ncg_monitor - Print relevant information about the ncg iteration.
474: Notes: This solves the nonlinear system of equations F(x) = 0 using the nonlinear generalization of the conjugate
475: gradient method. This may be used with a nonlinear preconditioner used to pick the new search directions, but otherwise
476: chooses the initial search direction as F(x) for the initial guess x.
479: .seealso: SNESCreate(), SNES, SNESSetType(), SNESLS, SNESTR, SNESNGMRES, SNESNQN
480: M*/
481: EXTERN_C_BEGIN
484: PetscErrorCode SNESCreate_NCG(SNES snes)
485: {
486: PetscErrorCode ierr;
487: SNES_NCG * neP;
490: snes->ops->destroy = SNESDestroy_NCG;
491: snes->ops->setup = SNESSetUp_NCG;
492: snes->ops->setfromoptions = SNESSetFromOptions_NCG;
493: snes->ops->view = SNESView_NCG;
494: snes->ops->solve = SNESSolve_NCG;
495: snes->ops->reset = SNESReset_NCG;
497: snes->usesksp = PETSC_FALSE;
498: snes->usespc = PETSC_TRUE;
500: if (!snes->tolerancesset) {
501: snes->max_funcs = 30000;
502: snes->max_its = 10000;
503: snes->stol = 1e-20;
504: }
506: PetscNewLog(snes, SNES_NCG, &neP);
507: snes->data = (void*) neP;
508: neP->monitor = PETSC_NULL;
509: neP->type = SNES_NCG_PRP;
510: PetscObjectComposeFunctionDynamic((PetscObject)snes,"SNESNCGSetType_C","SNESNCGSetType_NCG", SNESNCGSetType_NCG);
511: return(0);
512: }
513: EXTERN_C_END