1 #include <../src/snes/impls/ncg/snesncgimpl.h> /*I "petscsnes.h" I*/
2 const char *const SNESNCGTypes[] = {"FR", "PRP", "HS", "DY", "CD", "SNESNCGType", "SNES_NCG_", NULL};
3
SNESDestroy_NCG(SNES snes)4 static PetscErrorCode SNESDestroy_NCG(SNES snes)
5 {
6 PetscFunctionBegin;
7 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNCGSetType_C", NULL));
8 PetscCall(PetscFree(snes->data));
9 PetscFunctionReturn(PETSC_SUCCESS);
10 }
11
SNESSetUp_NCG(SNES snes)12 static PetscErrorCode SNESSetUp_NCG(SNES snes)
13 {
14 PetscFunctionBegin;
15 PetscCall(SNESSetWorkVecs(snes, 2));
16 PetscCheck(snes->npcside != PC_RIGHT, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNESNCG only supports left preconditioning");
17 if (snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED;
18 PetscFunctionReturn(PETSC_SUCCESS);
19 }
20
SNESLineSearchApply_NCGLinear(SNESLineSearch linesearch)21 static PetscErrorCode SNESLineSearchApply_NCGLinear(SNESLineSearch linesearch)
22 {
23 PetscScalar alpha, ptAp;
24 Vec X, Y, F, W;
25 SNES snes;
26 PetscReal *fnorm, *xnorm, *ynorm;
27
28 PetscFunctionBegin;
29 PetscCall(SNESLineSearchGetSNES(linesearch, &snes));
30 X = linesearch->vec_sol;
31 W = linesearch->vec_sol_new;
32 F = linesearch->vec_func;
33 Y = linesearch->vec_update;
34 fnorm = &linesearch->fnorm;
35 xnorm = &linesearch->xnorm;
36 ynorm = &linesearch->ynorm;
37
38 if (!snes->jacobian) PetscCall(SNESSetUpMatrices(snes));
39
40 /*
41 The exact step size for unpreconditioned linear CG is just:
42 alpha = (r, r) / (p, Ap) = (f, f) / (y, Jy)
43 */
44 PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
45 SNESLineSearchCheckJacobianDomainError(snes, linesearch);
46
47 PetscCall(VecDot(F, F, &alpha));
48 PetscCall(MatMult(snes->jacobian, Y, W));
49 PetscCall(VecDot(Y, W, &ptAp));
50 alpha = alpha / ptAp;
51 PetscCall(VecAXPY(X, -alpha, Y));
52 PetscCall(SNESComputeFunction(snes, X, F));
53
54 PetscCall(VecNorm(F, NORM_2, fnorm));
55 SNESLineSearchCheckFunctionDomainError(snes, linesearch, *fnorm);
56 PetscCall(VecNorm(X, NORM_2, xnorm));
57 PetscCall(VecNorm(Y, NORM_2, ynorm));
58 PetscFunctionReturn(PETSC_SUCCESS);
59 }
60
61 /*MC
62 SNESLINESEARCHNCGLINEAR - Special line search only for the nonlinear CG solver `SNESNCG`
63
64 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.
65 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.
66
67 Level: advanced
68
69 Notes:
70 This requires a Jacobian-vector product but does not require the solution of a linear system with the Jacobian
71
72 This is a "odd-ball" line search, we don't know if it is in the literature or used in practice by anyone.
73
74 .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESLineSearchCreate()`, `SNESLineSearchSetType()`
75 M*/
76
SNESLineSearchCreate_NCGLinear(SNESLineSearch linesearch)77 PETSC_EXTERN PetscErrorCode SNESLineSearchCreate_NCGLinear(SNESLineSearch linesearch)
78 {
79 PetscFunctionBegin;
80 linesearch->ops->apply = SNESLineSearchApply_NCGLinear;
81 linesearch->ops->destroy = NULL;
82 linesearch->ops->setfromoptions = NULL;
83 linesearch->ops->reset = NULL;
84 linesearch->ops->view = NULL;
85 linesearch->ops->setup = NULL;
86 PetscFunctionReturn(PETSC_SUCCESS);
87 }
88
SNESSetFromOptions_NCG(SNES snes,PetscOptionItems PetscOptionsObject)89 static PetscErrorCode SNESSetFromOptions_NCG(SNES snes, PetscOptionItems PetscOptionsObject)
90 {
91 SNES_NCG *ncg = (SNES_NCG *)snes->data;
92 PetscBool debug = PETSC_FALSE;
93 SNESNCGType ncgtype = ncg->type;
94 SNESLineSearch linesearch;
95
96 PetscFunctionBegin;
97 PetscOptionsHeadBegin(PetscOptionsObject, "SNES NCG options");
98 PetscCall(PetscOptionsBool("-snes_ncg_monitor", "Monitor the beta values used in the NCG iterations", "SNES", ncg->monitor ? PETSC_TRUE : PETSC_FALSE, &debug, NULL));
99 if (debug) ncg->monitor = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));
100 PetscCall(PetscOptionsEnum("-snes_ncg_type", "NCG Beta type used", "SNESNCGSetType", SNESNCGTypes, (PetscEnum)ncg->type, (PetscEnum *)&ncgtype, NULL));
101 PetscCall(SNESNCGSetType(snes, ncgtype));
102 PetscOptionsHeadEnd();
103 if (!snes->linesearch) {
104 PetscCall(SNESGetLineSearch(snes, &linesearch));
105 if (!((PetscObject)linesearch)->type_name) {
106 if (!snes->npc) {
107 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHCP));
108 } else {
109 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHSECANT));
110 }
111 }
112 }
113 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
SNESView_NCG(SNES snes,PetscViewer viewer)116 static PetscErrorCode SNESView_NCG(SNES snes, PetscViewer viewer)
117 {
118 SNES_NCG *ncg = (SNES_NCG *)snes->data;
119 PetscBool isascii;
120
121 PetscFunctionBegin;
122 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
123 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, " type: %s\n", SNESNCGTypes[ncg->type]));
124 PetscFunctionReturn(PETSC_SUCCESS);
125 }
126
127 /*@
128 SNESNCGSetType - Sets the conjugate update type for nonlinear CG `SNESNCG`.
129
130 Logically Collective
131
132 Input Parameters:
133 + snes - the iterative context
134 - btype - update type, see `SNESNCGType`
135
136 Options Database Key:
137 . -snes_ncg_type <prp,fr,hs,dy,cd> - strategy for selecting algorithm for computing beta
138
139 Level: intermediate
140
141 Notes:
142 `SNES_NCG_PRP` is the default, and the only one that tolerates generalized search directions.
143
144 It is not clear what "generalized search directions" means, does it mean use with a nonlinear preconditioner,
145 that is using -npc_snes_type <type>, `SNESSetNPC()`, or `SNESGetNPC()`?
146
147 .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESNCGType`, `SNES_NCG_FR`, `SNES_NCG_PRP`, `SNES_NCG_HS`, `SNES_NCG_DY`, `SNES_NCG_CD`
148 @*/
SNESNCGSetType(SNES snes,SNESNCGType btype)149 PetscErrorCode SNESNCGSetType(SNES snes, SNESNCGType btype)
150 {
151 PetscFunctionBegin;
152 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
153 PetscTryMethod(snes, "SNESNCGSetType_C", (SNES, SNESNCGType), (snes, btype));
154 PetscFunctionReturn(PETSC_SUCCESS);
155 }
156
SNESNCGSetType_NCG(SNES snes,SNESNCGType btype)157 static PetscErrorCode SNESNCGSetType_NCG(SNES snes, SNESNCGType btype)
158 {
159 SNES_NCG *ncg = (SNES_NCG *)snes->data;
160
161 PetscFunctionBegin;
162 ncg->type = btype;
163 PetscFunctionReturn(PETSC_SUCCESS);
164 }
165
166 /*
167 SNESSolve_NCG - Solves a nonlinear system with the Nonlinear Conjugate Gradient method.
168
169 Input Parameter:
170 . snes - the `SNES` context
171
172 Output Parameter:
173 . outits - number of iterations until termination
174
175 Application Interface Routine: SNESSolve()
176 */
SNESSolve_NCG(SNES snes)177 static PetscErrorCode SNESSolve_NCG(SNES snes)
178 {
179 SNES_NCG *ncg = (SNES_NCG *)snes->data;
180 Vec X, dX, lX, F, dXold;
181 PetscReal fnorm, ynorm, xnorm, beta = 0.0;
182 PetscScalar dXdotdX, dXolddotdXold, dXdotdXold, lXdotdX, lXdotdXold;
183 PetscInt maxits, i;
184 SNESLineSearch linesearch;
185 SNESConvergedReason reason;
186
187 PetscFunctionBegin;
188 PetscCheck(!snes->xl && !snes->xu && !snes->ops->computevariablebounds, PetscObjectComm((PetscObject)snes), PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);
189
190 PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
191 snes->reason = SNES_CONVERGED_ITERATING;
192
193 maxits = snes->max_its; /* maximum number of iterations */
194 X = snes->vec_sol; /* X^n */
195 dXold = snes->work[0]; /* The previous iterate of X */
196 dX = snes->work[1]; /* the preconditioned direction */
197 lX = snes->vec_sol_update; /* the conjugate direction */
198 F = snes->vec_func; /* residual vector */
199
200 PetscCall(SNESGetLineSearch(snes, &linesearch));
201
202 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
203 snes->iter = 0;
204 snes->norm = 0.;
205 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
206
207 /* compute the initial function and preconditioned update dX */
208
209 if (snes->npc && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
210 PetscCall(SNESApplyNPC(snes, X, NULL, dX));
211 PetscCall(SNESGetConvergedReason(snes->npc, &reason));
212 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
213 snes->reason = SNES_DIVERGED_INNER;
214 PetscFunctionReturn(PETSC_SUCCESS);
215 }
216 PetscCall(VecCopy(dX, F));
217 PetscCall(VecNorm(F, NORM_2, &fnorm));
218 } else {
219 if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
220 else snes->vec_func_init_set = PETSC_FALSE;
221
222 /* convergence test */
223 PetscCall(VecNorm(F, NORM_2, &fnorm));
224 SNESCheckFunctionDomainError(snes, fnorm);
225 PetscCall(VecCopy(F, dX));
226 }
227 if (snes->npc) {
228 if (snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
229 PetscCall(SNESApplyNPC(snes, X, F, dX));
230 PetscCall(SNESGetConvergedReason(snes->npc, &reason));
231 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
232 snes->reason = SNES_DIVERGED_INNER;
233 PetscFunctionReturn(PETSC_SUCCESS);
234 }
235 }
236 }
237 PetscCall(VecCopy(dX, lX));
238 PetscCall(VecDot(dX, dX, &dXdotdX));
239
240 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
241 snes->norm = fnorm;
242 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
243 PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
244
245 /* test convergence */
246 PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
247 PetscCall(SNESMonitor(snes, 0, fnorm));
248 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
249
250 /* Call general purpose update function */
251 PetscTryTypeMethod(snes, update, snes->iter);
252
253 /* first update -- just use the (preconditioned) residual direction for the initial conjugate direction */
254
255 for (i = 1; i < maxits + 1; i++) {
256 /* some update types require the old update direction or conjugate direction */
257 if (ncg->type != SNES_NCG_FR) PetscCall(VecCopy(dX, dXold));
258 PetscCall(SNESLineSearchApply(linesearch, X, F, &fnorm, lX));
259 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
260 PetscCall(SNESLineSearchGetNorms(linesearch, &xnorm, &fnorm, &ynorm));
261 SNESCheckLineSearchFailure(snes);
262 if (snes->nfuncs >= snes->max_funcs && snes->max_funcs >= 0) {
263 snes->reason = SNES_DIVERGED_FUNCTION_COUNT;
264 PetscFunctionReturn(PETSC_SUCCESS);
265 }
266 /* Monitor convergence */
267 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
268 snes->iter = i;
269 snes->norm = fnorm;
270 snes->xnorm = xnorm;
271 snes->ynorm = ynorm;
272 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
273 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, 0));
274
275 /* Test for convergence */
276 PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
277 PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
278 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
279
280 /* Call general purpose update function */
281 PetscTryTypeMethod(snes, update, snes->iter);
282 if (snes->npc) {
283 if (snes->functype == SNES_FUNCTION_PRECONDITIONED) {
284 PetscCall(SNESApplyNPC(snes, X, NULL, dX));
285 PetscCall(SNESGetConvergedReason(snes->npc, &reason));
286 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
287 snes->reason = SNES_DIVERGED_INNER;
288 PetscFunctionReturn(PETSC_SUCCESS);
289 }
290 PetscCall(VecCopy(dX, F));
291 } else {
292 PetscCall(SNESApplyNPC(snes, X, F, dX));
293 PetscCall(SNESGetConvergedReason(snes->npc, &reason));
294 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
295 snes->reason = SNES_DIVERGED_INNER;
296 PetscFunctionReturn(PETSC_SUCCESS);
297 }
298 }
299 } else {
300 PetscCall(VecCopy(F, dX));
301 }
302
303 /* compute the conjugate direction lX = dX + beta*lX with beta = ((dX, dX) / (dX_old, dX_old) (Fletcher-Reeves update)*/
304 switch (ncg->type) {
305 case SNES_NCG_FR: /* Fletcher-Reeves */
306 dXolddotdXold = dXdotdX;
307 PetscCall(VecDot(dX, dX, &dXdotdX));
308 beta = PetscRealPart(dXdotdX / dXolddotdXold);
309 break;
310 case SNES_NCG_PRP: /* Polak-Ribiere-Poylak */
311 dXolddotdXold = dXdotdX;
312 PetscCall(VecDotBegin(dX, dX, &dXdotdX));
313 PetscCall(VecDotBegin(dXold, dX, &dXdotdXold));
314 PetscCall(VecDotEnd(dX, dX, &dXdotdX));
315 PetscCall(VecDotEnd(dXold, dX, &dXdotdXold));
316 beta = PetscRealPart((dXdotdX - dXdotdXold) / dXolddotdXold);
317 if (beta < 0.0) beta = 0.0; /* restart */
318 break;
319 case SNES_NCG_HS: /* Hestenes-Stiefel */
320 PetscCall(VecDotBegin(dX, dX, &dXdotdX));
321 PetscCall(VecDotBegin(dX, dXold, &dXdotdXold));
322 PetscCall(VecDotBegin(lX, dX, &lXdotdX));
323 PetscCall(VecDotBegin(lX, dXold, &lXdotdXold));
324 PetscCall(VecDotEnd(dX, dX, &dXdotdX));
325 PetscCall(VecDotEnd(dX, dXold, &dXdotdXold));
326 PetscCall(VecDotEnd(lX, dX, &lXdotdX));
327 PetscCall(VecDotEnd(lX, dXold, &lXdotdXold));
328 beta = PetscRealPart((dXdotdX - dXdotdXold) / (lXdotdX - lXdotdXold));
329 break;
330 case SNES_NCG_DY: /* Dai-Yuan */
331 PetscCall(VecDotBegin(dX, dX, &dXdotdX));
332 PetscCall(VecDotBegin(lX, dX, &lXdotdX));
333 PetscCall(VecDotBegin(lX, dXold, &lXdotdXold));
334 PetscCall(VecDotEnd(dX, dX, &dXdotdX));
335 PetscCall(VecDotEnd(lX, dX, &lXdotdX));
336 PetscCall(VecDotEnd(lX, dXold, &lXdotdXold));
337 beta = PetscRealPart(dXdotdX / (lXdotdXold - lXdotdX));
338 break;
339 case SNES_NCG_CD: /* Conjugate Descent */
340 PetscCall(VecDotBegin(dX, dX, &dXdotdX));
341 PetscCall(VecDotBegin(lX, dXold, &lXdotdXold));
342 PetscCall(VecDotEnd(dX, dX, &dXdotdX));
343 PetscCall(VecDotEnd(lX, dXold, &lXdotdXold));
344 beta = PetscRealPart(dXdotdX / lXdotdXold);
345 break;
346 }
347 if (ncg->monitor) PetscCall(PetscViewerASCIIPrintf(ncg->monitor, "beta = %e\n", (double)beta));
348 PetscCall(VecAYPX(lX, beta, dX));
349 }
350 PetscFunctionReturn(PETSC_SUCCESS);
351 }
352
353 /*MC
354 SNESNCG - Nonlinear Conjugate-Gradient method for the solution of nonlinear systems {cite}`bruneknepleysmithtu15`.
355
356 Level: beginner
357
358 Options Database Keys:
359 + -snes_ncg_type <fr, prp, dy, hs, cd> - Choice of conjugate-gradient update parameter, default is `prp`.
360 . -snes_linesearch_type <cp,l2,basic> - Line search type.
361 - -snes_ncg_monitor - Print the beta values nonlinear Conjugate-Gradient used in the iteration, .
362
363 Notes:
364 This solves the nonlinear system of equations $ F(x) = 0 $ using the nonlinear generalization of the conjugate
365 gradient method. This may be used with a nonlinear preconditioner used to pick the new search directions, but otherwise
366 chooses the initial search direction as $ F(x) $ for the initial guess $x$.
367
368 Only supports left non-linear preconditioning.
369
370 Default line search is `SNESLINESEARCHCP`, unless a nonlinear preconditioner is used with `-npc_snes_type` <type>, `SNESSetNPC()`, or `SNESGetNPC()` then
371 `SNESLINESEARCHSECANT` is used. Also supports the special-purpose line search `SNESLINESEARCHNCGLINEAR`
372
373 .seealso: [](ch_snes), `SNES`, `SNESNCG`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`, `SNESNGMRES`, `SNESQN`, `SNESLINESEARCHNCGLINEAR`, `SNESNCGSetType()`, `SNESLineSearchSetType()`
374 M*/
SNESCreate_NCG(SNES snes)375 PETSC_EXTERN PetscErrorCode SNESCreate_NCG(SNES snes)
376 {
377 SNES_NCG *neP;
378
379 PetscFunctionBegin;
380 snes->ops->destroy = SNESDestroy_NCG;
381 snes->ops->setup = SNESSetUp_NCG;
382 snes->ops->setfromoptions = SNESSetFromOptions_NCG;
383 snes->ops->view = SNESView_NCG;
384 snes->ops->solve = SNESSolve_NCG;
385
386 snes->usesksp = PETSC_FALSE;
387 snes->usesnpc = PETSC_TRUE;
388 snes->npcside = PC_LEFT;
389
390 snes->alwayscomputesfinalresidual = PETSC_TRUE;
391
392 PetscCall(SNESParametersInitialize(snes));
393 PetscObjectParameterSetDefault(snes, max_funcs, 30000);
394 PetscObjectParameterSetDefault(snes, max_its, 10000);
395 PetscObjectParameterSetDefault(snes, stol, 1e-20);
396
397 PetscCall(PetscNew(&neP));
398 snes->data = (void *)neP;
399 neP->monitor = NULL;
400 neP->type = SNES_NCG_PRP;
401 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESNCGSetType_C", SNESNCGSetType_NCG));
402 PetscFunctionReturn(PETSC_SUCCESS);
403 }
404