Actual source code: anderson.c
1: #include <../src/snes/impls/ngmres/snesngmres.h>
3: static PetscErrorCode SNESSetFromOptions_Anderson(SNES snes, PetscOptionItems *PetscOptionsObject)
4: {
5: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
6: PetscBool monitor = PETSC_FALSE;
8: PetscOptionsHeadBegin(PetscOptionsObject, "SNES NGMRES options");
9: PetscOptionsInt("-snes_anderson_m", "Number of directions", "SNES", ngmres->msize, &ngmres->msize, NULL);
10: PetscOptionsReal("-snes_anderson_beta", "Mixing parameter", "SNES", ngmres->andersonBeta, &ngmres->andersonBeta, NULL);
11: PetscOptionsInt("-snes_anderson_restart", "Iterations before forced restart", "SNES", ngmres->restart_periodic, &ngmres->restart_periodic, NULL);
12: PetscOptionsInt("-snes_anderson_restart_it", "Tolerance iterations before restart", "SNES", ngmres->restart_it, &ngmres->restart_it, NULL);
13: PetscOptionsEnum("-snes_anderson_restart_type", "Restart type", "SNESNGMRESSetRestartType", SNESNGMRESRestartTypes, (PetscEnum)ngmres->restart_type, (PetscEnum *)&ngmres->restart_type, NULL);
14: PetscOptionsBool("-snes_anderson_monitor", "Monitor steps of Anderson Mixing", "SNES", ngmres->monitor ? PETSC_TRUE : PETSC_FALSE, &monitor, NULL);
15: if (monitor) ngmres->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
16: PetscOptionsHeadEnd();
17: return 0;
18: }
20: static PetscErrorCode SNESSolve_Anderson(SNES snes)
21: {
22: SNES_NGMRES *ngmres = (SNES_NGMRES *)snes->data;
23: /* present solution, residual, and preconditioned residual */
24: Vec X, F, B, D;
25: /* candidate linear combination answers */
26: Vec XA, FA, XM, FM;
28: /* coefficients and RHS to the minimization problem */
29: PetscReal fnorm, fMnorm, fAnorm;
30: PetscReal xnorm, ynorm;
31: PetscReal dnorm, dminnorm = 0.0, fminnorm;
32: PetscInt restart_count = 0;
33: PetscInt k, k_restart, l, ivec;
34: PetscBool selectRestart;
35: SNESConvergedReason reason;
39: PetscCitationsRegister(SNESCitation, &SNEScite);
40: /* variable initialization */
41: snes->reason = SNES_CONVERGED_ITERATING;
42: X = snes->vec_sol;
43: F = snes->vec_func;
44: B = snes->vec_rhs;
45: XA = snes->vec_sol_update;
46: FA = snes->work[0];
47: D = snes->work[1];
49: /* work for the line search */
50: XM = snes->work[3];
51: FM = snes->work[4];
53: PetscObjectSAWsTakeAccess((PetscObject)snes);
54: snes->iter = 0;
55: snes->norm = 0.;
56: PetscObjectSAWsGrantAccess((PetscObject)snes);
58: /* initialization */
60: /* r = F(x) */
62: if (snes->npc && snes->npcside == PC_LEFT) {
63: SNESApplyNPC(snes, X, NULL, F);
64: SNESGetConvergedReason(snes->npc, &reason);
65: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
66: snes->reason = SNES_DIVERGED_INNER;
67: return 0;
68: }
69: VecNorm(F, NORM_2, &fnorm);
70: } else {
71: if (!snes->vec_func_init_set) {
72: SNESComputeFunction(snes, X, F);
73: } else snes->vec_func_init_set = PETSC_FALSE;
75: VecNorm(F, NORM_2, &fnorm);
76: SNESCheckFunctionNorm(snes, fnorm);
77: }
78: fminnorm = fnorm;
80: PetscObjectSAWsTakeAccess((PetscObject)snes);
81: snes->norm = fnorm;
82: PetscObjectSAWsGrantAccess((PetscObject)snes);
83: SNESLogConvergenceHistory(snes, fnorm, 0);
84: SNESMonitor(snes, 0, fnorm);
85: PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
86: if (snes->reason) return 0;
88: k_restart = 0;
89: l = 0;
90: ivec = 0;
91: for (k = 1; k < snes->max_its + 1; k++) {
92: /* select which vector of the stored subspace will be updated */
93: if (snes->npc && snes->npcside == PC_RIGHT) {
94: VecCopy(X, XM);
95: SNESSetInitialFunction(snes->npc, F);
97: PetscLogEventBegin(SNES_NPCSolve, snes->npc, XM, B, 0);
98: SNESSolve(snes->npc, B, XM);
99: PetscLogEventEnd(SNES_NPCSolve, snes->npc, XM, B, 0);
101: SNESGetConvergedReason(snes->npc, &reason);
102: if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
103: snes->reason = SNES_DIVERGED_INNER;
104: return 0;
105: }
106: SNESGetNPCFunction(snes, FM, &fMnorm);
107: if (ngmres->andersonBeta != 1.0) VecAXPBY(XM, (1.0 - ngmres->andersonBeta), ngmres->andersonBeta, X);
108: } else {
109: VecCopy(F, FM);
110: VecCopy(X, XM);
111: VecAXPY(XM, -ngmres->andersonBeta, FM);
112: fMnorm = fnorm;
113: }
115: SNESNGMRESFormCombinedSolution_Private(snes, ivec, l, XM, FM, fMnorm, X, XA, FA);
116: ivec = k_restart % ngmres->msize;
117: if (ngmres->restart_type == SNES_NGMRES_RESTART_DIFFERENCE) {
118: SNESNGMRESNorms_Private(snes, l, X, F, XM, FM, XA, FA, D, &dnorm, &dminnorm, NULL, NULL, NULL, &xnorm, &fAnorm, &ynorm);
119: SNESNGMRESSelectRestart_Private(snes, l, fMnorm, fnorm, dnorm, fminnorm, dminnorm, &selectRestart);
120: /* if the restart conditions persist for more than restart_it iterations, restart. */
121: if (selectRestart) restart_count++;
122: else restart_count = 0;
123: } else if (ngmres->restart_type == SNES_NGMRES_RESTART_PERIODIC) {
124: SNESNGMRESNorms_Private(snes, l, X, F, XM, FM, XA, FA, D, NULL, NULL, NULL, NULL, NULL, &xnorm, &fAnorm, &ynorm);
125: if (k_restart > ngmres->restart_periodic) {
126: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "periodic restart after %" PetscInt_FMT " iterations\n", k_restart);
127: restart_count = ngmres->restart_it;
128: }
129: } else {
130: SNESNGMRESNorms_Private(snes, l, X, F, XM, FM, XA, FA, D, NULL, NULL, NULL, NULL, NULL, &xnorm, &fAnorm, &ynorm);
131: }
132: /* restart after restart conditions have persisted for a fixed number of iterations */
133: if (restart_count >= ngmres->restart_it) {
134: if (ngmres->monitor) PetscViewerASCIIPrintf(ngmres->monitor, "Restarted at iteration %" PetscInt_FMT "\n", k_restart);
135: restart_count = 0;
136: k_restart = 0;
137: l = 0;
138: ivec = 0;
139: } else {
140: if (l < ngmres->msize) l++;
141: k_restart++;
142: SNESNGMRESUpdateSubspace_Private(snes, ivec, l, FM, fnorm, XM);
143: }
145: fnorm = fAnorm;
146: if (fminnorm > fnorm) fminnorm = fnorm;
148: VecCopy(XA, X);
149: VecCopy(FA, F);
151: PetscObjectSAWsTakeAccess((PetscObject)snes);
152: snes->iter = k;
153: snes->norm = fnorm;
154: snes->xnorm = xnorm;
155: snes->ynorm = ynorm;
156: PetscObjectSAWsGrantAccess((PetscObject)snes);
157: SNESLogConvergenceHistory(snes, snes->norm, snes->iter);
158: SNESMonitor(snes, snes->iter, snes->norm);
159: PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
160: if (snes->reason) return 0;
161: }
162: snes->reason = SNES_DIVERGED_MAX_IT;
163: return 0;
164: }
166: /*MC
167: SNESANDERSON - Anderson Mixing nonlinear solver
169: Level: beginner
171: Options Database Keys:
172: + -snes_anderson_m - Number of stored previous solutions and residuals
173: . -snes_anderson_beta - Anderson mixing parameter
174: . -snes_anderson_restart_type - Type of restart (see SNESNGMRES)
175: . -snes_anderson_restart_it - Number of iterations of restart conditions before restart
176: . -snes_anderson_restart - Number of iterations before periodic restart
177: - -snes_anderson_monitor - Prints relevant information about the ngmres iteration
179: Notes:
180: The Anderson Mixing method combines m previous solutions into a minimum-residual solution by solving a small linearized
181: optimization problem at each iteration.
183: Very similar to the `SNESNGMRES` algorithm.
185: References:
186: + * - D. G. Anderson. Iterative procedures for nonlinear integral equations.
187: J. Assoc. Comput. Mach., 12, 1965."
188: - * - Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu,"Composing Scalable Nonlinear Algebraic Solvers",
189: SIAM Review, 57(4), 2015
191: .seealso: `SNESNGMRES`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESType`
192: M*/
194: PETSC_EXTERN PetscErrorCode SNESCreate_Anderson(SNES snes)
195: {
196: SNES_NGMRES *ngmres;
197: SNESLineSearch linesearch;
199: snes->ops->destroy = SNESDestroy_NGMRES;
200: snes->ops->setup = SNESSetUp_NGMRES;
201: snes->ops->setfromoptions = SNESSetFromOptions_Anderson;
202: snes->ops->view = SNESView_NGMRES;
203: snes->ops->solve = SNESSolve_Anderson;
204: snes->ops->reset = SNESReset_NGMRES;
206: snes->usesnpc = PETSC_TRUE;
207: snes->usesksp = PETSC_FALSE;
208: snes->npcside = PC_RIGHT;
210: snes->alwayscomputesfinalresidual = PETSC_TRUE;
212: PetscNew(&ngmres);
213: snes->data = (void *)ngmres;
214: ngmres->msize = 30;
216: if (!snes->tolerancesset) {
217: snes->max_funcs = 30000;
218: snes->max_its = 10000;
219: }
221: SNESGetLineSearch(snes, &linesearch);
222: if (!((PetscObject)linesearch)->type_name) SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);
224: ngmres->additive_linesearch = NULL;
225: ngmres->approxfunc = PETSC_FALSE;
226: ngmres->restart_type = SNES_NGMRES_RESTART_NONE;
227: ngmres->restart_it = 2;
228: ngmres->restart_periodic = 30;
229: ngmres->gammaA = 2.0;
230: ngmres->gammaC = 2.0;
231: ngmres->deltaB = 0.9;
232: ngmres->epsilonB = 0.1;
234: ngmres->andersonBeta = 1.0;
235: return 0;
236: }