1 #include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
2 #include <petscdm.h>
3
4 const char *const SNESQNScaleTypes[] = {"DEFAULT", "NONE", "SCALAR", "DIAGONAL", "JACOBIAN", "SNESQNScaleType", "SNES_QN_SCALING_", NULL};
5 const char *const SNESQNRestartTypes[] = {"DEFAULT", "NONE", "POWELL", "PERIODIC", "SNESQNRestartType", "SNES_QN_RESTART_", NULL};
6 const char *const SNESQNTypes[] = {"LBFGS", "BROYDEN", "BADBROYDEN", "SNESQNType", "SNES_QN_", NULL};
7
8 typedef struct {
9 Mat B; /* Quasi-Newton approximation Matrix (MATLMVM) */
10 PetscInt m; /* The number of kept previous steps */
11 PetscReal *lambda; /* The line search history of the method */
12 PetscBool monflg;
13 PetscViewer monitor;
14 PetscReal powell_gamma; /* Powell angle restart condition */
15 PetscReal scaling; /* scaling of H0 */
16 SNESQNType type; /* the type of quasi-newton method used */
17 SNESQNScaleType scale_type; /* the type of scaling used */
18 SNESQNRestartType restart_type; /* determine the frequency and type of restart conditions */
19 } SNES_QN;
20
SNESQNGetMatrix_Private(SNES snes,Mat * B)21 static PetscErrorCode SNESQNGetMatrix_Private(SNES snes, Mat *B)
22 {
23 SNES_QN *qn = (SNES_QN *)snes->data;
24 const char *optionsprefix;
25
26 PetscFunctionBegin;
27 if (!qn->B) {
28 PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &qn->B));
29 PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
30 PetscCall(MatSetOptionsPrefix(qn->B, optionsprefix));
31 PetscCall(MatAppendOptionsPrefix(qn->B, "qn_"));
32 switch (qn->type) {
33 case SNES_QN_BROYDEN:
34 PetscCall(MatSetType(qn->B, MATLMVMBROYDEN));
35 qn->scale_type = SNES_QN_SCALE_NONE;
36 break;
37 case SNES_QN_BADBROYDEN:
38 PetscCall(MatSetType(qn->B, MATLMVMBADBROYDEN));
39 qn->scale_type = SNES_QN_SCALE_NONE;
40 break;
41 default:
42 PetscCall(MatSetType(qn->B, MATLMVMBFGS));
43 switch (qn->scale_type) {
44 case SNES_QN_SCALE_NONE:
45 case SNES_QN_SCALE_JACOBIAN:
46 PetscCall(MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_NONE));
47 break;
48 case SNES_QN_SCALE_SCALAR:
49 case SNES_QN_SCALE_DEFAULT:
50 PetscCall(MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_SCALAR));
51 break;
52 default:
53 break;
54 }
55 break;
56 }
57 }
58 *B = qn->B;
59 PetscFunctionReturn(PETSC_SUCCESS);
60 }
61
SNESSolve_QN(SNES snes)62 static PetscErrorCode SNESSolve_QN(SNES snes)
63 {
64 SNES_QN *qn = (SNES_QN *)snes->data;
65 Vec X, F, W;
66 Vec Y, D, Dold;
67 PetscInt i, i_r;
68 PetscReal fnorm, xnorm, ynorm;
69 SNESLineSearchReason lsreason;
70 PetscBool powell, periodic, restart;
71 PetscScalar DolddotD, DolddotDold;
72 SNESConvergedReason reason;
73 #if defined(PETSC_USE_INFO)
74 PetscReal gnorm;
75 #endif
76
77 PetscFunctionBegin;
78 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);
79
80 PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
81 F = snes->vec_func; /* residual vector */
82 Y = snes->vec_sol_update; /* search direction generated by J^-1D*/
83 X = snes->vec_sol; /* solution vector */
84
85 /* directions */
86 W = snes->work[0];
87 D = snes->work[1];
88 Dold = snes->work[2];
89
90 snes->reason = SNES_CONVERGED_ITERATING;
91
92 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
93 snes->iter = 0;
94 snes->norm = 0.;
95 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
96
97 if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
98 /* we need to sample the preconditioned function only once here,
99 since linesearch will always use the preconditioned function */
100 PetscCall(SNESApplyNPC(snes, X, NULL, F));
101 PetscCall(SNESGetConvergedReason(snes->npc, &reason));
102 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
103 snes->reason = SNES_DIVERGED_INNER;
104 PetscFunctionReturn(PETSC_SUCCESS);
105 }
106 } else {
107 if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
108 else snes->vec_func_init_set = PETSC_FALSE;
109 }
110 PetscCall(VecNorm(F, NORM_2, &fnorm));
111 SNESCheckFunctionDomainError(snes, fnorm);
112
113 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
114 snes->norm = fnorm;
115 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
116 PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
117 PetscCall(SNESMonitor(snes, 0, fnorm));
118
119 /* test convergence */
120 PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
121 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
122
123 restart = PETSC_TRUE;
124 for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
125 PetscScalar gdx;
126
127 /* general purpose update */
128 PetscTryTypeMethod(snes, update, snes->iter);
129
130 if (snes->npc) {
131 if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
132 PetscCall(SNESApplyNPC(snes, X, F, D));
133 PetscCall(SNESGetConvergedReason(snes->npc, &reason));
134 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
135 snes->reason = SNES_DIVERGED_INNER;
136 PetscFunctionReturn(PETSC_SUCCESS);
137 }
138 } else if (snes->npcside == PC_RIGHT) {
139 PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, 0, 0));
140 PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
141 PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, 0, 0));
142 PetscCall(SNESGetConvergedReason(snes->npc, &reason));
143 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
144 snes->reason = SNES_DIVERGED_INNER;
145 PetscFunctionReturn(PETSC_SUCCESS);
146 }
147 PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
148 PetscCall(VecCopy(F, D));
149 } else {
150 PetscCall(VecCopy(F, D));
151 }
152 } else {
153 PetscCall(VecCopy(F, D));
154 }
155
156 /* scale the initial update */
157 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN && restart) {
158 PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
159 SNESCheckJacobianDomainError(snes);
160 PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
161 PetscCall(MatLMVMSetJ0KSP(qn->B, snes->ksp));
162 }
163
164 /* update QN approx and calculate step */
165 PetscCall(MatLMVMUpdate(qn->B, X, D));
166 PetscCall(MatSolve(qn->B, D, Y));
167 PetscCall(VecDot(F, Y, &gdx));
168 if (PetscAbsScalar(gdx) <= 0) {
169 /* Step is not descent or solve was not successful
170 Use steepest descent direction */
171 PetscCall(VecCopy(F, Y));
172 PetscCall(MatLMVMReset(qn->B, PETSC_FALSE));
173 }
174
175 /* line search for lambda */
176 ynorm = 1;
177 #if defined(PETSC_USE_INFO)
178 gnorm = fnorm;
179 #endif
180 PetscCall(VecCopy(D, Dold));
181 PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y));
182 if (snes->reason) break;
183 SNESCheckLineSearchFailure(snes);
184 PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsreason));
185
186 PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm));
187 /* convergence monitoring */
188 PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lsreason));
189
190 PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
191 snes->iter = i + 1;
192 snes->norm = fnorm;
193 snes->xnorm = xnorm;
194 snes->ynorm = ynorm;
195 PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
196
197 PetscCall(SNESLogConvergenceHistory(snes, snes->norm, snes->iter));
198
199 /* set parameter for default relative tolerance convergence test */
200 PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
201 PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
202 if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
203
204 if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
205 PetscCall(SNESApplyNPC(snes, X, F, D));
206 PetscCall(SNESGetConvergedReason(snes->npc, &reason));
207 if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
208 snes->reason = SNES_DIVERGED_INNER;
209 PetscFunctionReturn(PETSC_SUCCESS);
210 }
211 } else {
212 PetscCall(VecCopy(F, D));
213 }
214
215 /* restart conditions */
216 powell = PETSC_FALSE;
217 if (qn->restart_type == SNES_QN_RESTART_POWELL && i_r > 1) {
218 /* check restart by Powell's Criterion: |F^T H_0 Fold| > powell_gamma * |Fold^T H_0 Fold| */
219 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
220 PetscCall(MatMult(snes->jacobian, Dold, W));
221 } else {
222 PetscCall(VecCopy(Dold, W));
223 }
224 PetscCall(VecDotBegin(W, Dold, &DolddotDold));
225 PetscCall(VecDotBegin(W, D, &DolddotD));
226 PetscCall(VecDotEnd(W, Dold, &DolddotDold));
227 PetscCall(VecDotEnd(W, D, &DolddotD));
228 if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma * PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE;
229 }
230 periodic = PETSC_FALSE;
231 if (qn->restart_type == SNES_QN_RESTART_PERIODIC) {
232 if (i_r > qn->m - 1) periodic = PETSC_TRUE;
233 }
234
235 /* restart if either powell or periodic restart is satisfied. */
236 restart = PETSC_FALSE;
237 if (lsreason || powell || periodic) {
238 restart = PETSC_TRUE;
239 if (qn->monflg) {
240 PetscCall(PetscViewerASCIIAddTab(qn->monitor, ((PetscObject)snes)->tablevel + 2));
241 if (powell) {
242 PetscCall(PetscViewerASCIIPrintf(qn->monitor, "Powell restart! |%14.12e| > %6.4f*|%14.12e| i_r = %" PetscInt_FMT "\n", (double)PetscRealPart(DolddotD), (double)qn->powell_gamma, (double)PetscRealPart(DolddotDold), i_r));
243 } else {
244 PetscCall(PetscViewerASCIIPrintf(qn->monitor, "Periodic restart! i_r = %" PetscInt_FMT "\n", i_r));
245 }
246 PetscCall(PetscViewerASCIISubtractTab(qn->monitor, ((PetscObject)snes)->tablevel + 2));
247 }
248 i_r = -1;
249 PetscCall(MatLMVMReset(qn->B, PETSC_FALSE));
250 }
251 }
252 PetscFunctionReturn(PETSC_SUCCESS);
253 }
254
SNESSetUp_QN(SNES snes)255 static PetscErrorCode SNESSetUp_QN(SNES snes)
256 {
257 SNES_QN *qn = (SNES_QN *)snes->data;
258 DM dm;
259 PetscInt n, N;
260
261 PetscFunctionBegin;
262 if (!snes->vec_sol) {
263 PetscCall(SNESGetDM(snes, &dm));
264 PetscCall(DMCreateGlobalVector(dm, &snes->vec_sol));
265 }
266 PetscCall(SNESSetWorkVecs(snes, 3));
267 PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));
268
269 if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) PetscCall(SNESSetUpMatrices(snes));
270 if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;
271
272 /* set method defaults */
273 if (qn->scale_type == SNES_QN_SCALE_DEFAULT) {
274 if (qn->type == SNES_QN_BADBROYDEN) {
275 qn->scale_type = SNES_QN_SCALE_NONE;
276 } else {
277 qn->scale_type = SNES_QN_SCALE_SCALAR;
278 }
279 }
280 if (qn->restart_type == SNES_QN_RESTART_DEFAULT) {
281 if (qn->type == SNES_QN_LBFGS) {
282 qn->restart_type = SNES_QN_RESTART_POWELL;
283 } else {
284 qn->restart_type = SNES_QN_RESTART_PERIODIC;
285 }
286 }
287 /* Set up the LMVM matrix */
288 PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));
289 PetscCall(VecGetLocalSize(snes->vec_sol, &n));
290 PetscCall(VecGetSize(snes->vec_sol, &N));
291 PetscCall(MatSetSizes(qn->B, n, n, N, N));
292 PetscCall(MatSetUp(qn->B));
293 PetscCall(MatLMVMReset(qn->B, PETSC_TRUE));
294 PetscCall(MatLMVMSetHistorySize(qn->B, qn->m));
295 PetscCall(MatLMVMAllocate(qn->B, snes->vec_sol, snes->vec_func));
296 PetscFunctionReturn(PETSC_SUCCESS);
297 }
298
SNESReset_QN(SNES snes)299 static PetscErrorCode SNESReset_QN(SNES snes)
300 {
301 SNES_QN *qn = (SNES_QN *)snes->data;
302
303 PetscFunctionBegin;
304 PetscCall(MatDestroy(&qn->B));
305 PetscFunctionReturn(PETSC_SUCCESS);
306 }
307
SNESDestroy_QN(SNES snes)308 static PetscErrorCode SNESDestroy_QN(SNES snes)
309 {
310 PetscFunctionBegin;
311 PetscCall(SNESReset_QN(snes));
312 PetscCall(PetscFree(snes->data));
313 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", NULL));
314 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", NULL));
315 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", NULL));
316 PetscFunctionReturn(PETSC_SUCCESS);
317 }
318
SNESSetFromOptions_QN(SNES snes,PetscOptionItems PetscOptionsObject)319 static PetscErrorCode SNESSetFromOptions_QN(SNES snes, PetscOptionItems PetscOptionsObject)
320 {
321 SNES_QN *qn = (SNES_QN *)snes->data;
322 PetscBool flg;
323 SNESLineSearch linesearch;
324 SNESQNRestartType rtype = qn->restart_type;
325 SNESQNScaleType stype = qn->scale_type;
326 SNESQNType qtype = qn->type;
327
328 PetscFunctionBegin;
329 PetscOptionsHeadBegin(PetscOptionsObject, "SNES QN options");
330 PetscCall(PetscOptionsInt("-snes_qn_m", "Number of past states saved for L-BFGS methods", "SNESQN", qn->m, &qn->m, NULL));
331 PetscCall(PetscOptionsReal("-snes_qn_powell_gamma", "Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL));
332 PetscCall(PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", qn->monflg, &qn->monflg, NULL));
333 PetscCall(PetscOptionsEnum("-snes_qn_scale_type", "Scaling type", "SNESQNSetScaleType", SNESQNScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg));
334 if (flg) PetscCall(SNESQNSetScaleType(snes, stype));
335
336 PetscCall(PetscOptionsEnum("-snes_qn_restart_type", "Restart type", "SNESQNSetRestartType", SNESQNRestartTypes, (PetscEnum)rtype, (PetscEnum *)&rtype, &flg));
337 if (flg) PetscCall(SNESQNSetRestartType(snes, rtype));
338
339 PetscCall(PetscOptionsEnum("-snes_qn_type", "Quasi-Newton update type", "", SNESQNTypes, (PetscEnum)qtype, (PetscEnum *)&qtype, &flg));
340 if (flg) PetscCall(SNESQNSetType(snes, qtype));
341 PetscOptionsHeadEnd();
342 PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));
343 PetscCall(MatSetFromOptions(qn->B));
344 if (!snes->linesearch) {
345 PetscCall(SNESGetLineSearch(snes, &linesearch));
346 if (!((PetscObject)linesearch)->type_name) {
347 if (qn->type == SNES_QN_LBFGS) {
348 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHCP));
349 } else if (qn->type == SNES_QN_BROYDEN) {
350 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC));
351 } else {
352 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHSECANT));
353 }
354 }
355 }
356 if (qn->monflg) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &qn->monitor));
357 PetscFunctionReturn(PETSC_SUCCESS);
358 }
359
SNESView_QN(SNES snes,PetscViewer viewer)360 static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
361 {
362 SNES_QN *qn = (SNES_QN *)snes->data;
363 PetscBool isascii;
364
365 PetscFunctionBegin;
366 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
367 if (isascii) {
368 PetscCall(PetscViewerASCIIPrintf(viewer, " type is %s, restart type is %s, scale type is %s\n", SNESQNTypes[qn->type], SNESQNRestartTypes[qn->restart_type], SNESQNScaleTypes[qn->scale_type]));
369 PetscCall(PetscViewerASCIIPrintf(viewer, " Stored subspace size: %" PetscInt_FMT "\n", qn->m));
370 PetscCall(MatView(qn->B, viewer));
371 }
372 PetscFunctionReturn(PETSC_SUCCESS);
373 }
374
375 /*@
376 SNESQNSetRestartType - Sets the restart type for `SNESQN`.
377
378 Logically Collective
379
380 Input Parameters:
381 + snes - the iterative context
382 - rtype - restart type, see `SNESQNRestartType`
383
384 Options Database Keys:
385 + -snes_qn_restart_type <powell,periodic,none> - set the restart type
386 - -snes_qn_m <m> - sets the number of stored updates and the restart period for periodic
387
388 Level: intermediate
389
390 .seealso: [](ch_snes), `SNES`, `SNESQN`, `SNESQNRestartType`, `SNES_QN_RESTART_NONE`, `SNES_QN_RESTART_POWELL`, `SNES_QN_RESTART_PERIODIC`,
391 `SNESQNType`, `SNESQNScaleType`
392 @*/
SNESQNSetRestartType(SNES snes,SNESQNRestartType rtype)393 PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
394 {
395 PetscFunctionBegin;
396 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
397 PetscTryMethod(snes, "SNESQNSetRestartType_C", (SNES, SNESQNRestartType), (snes, rtype));
398 PetscFunctionReturn(PETSC_SUCCESS);
399 }
400
401 /*@
402 SNESQNSetScaleType - Sets the scaling type for the inner inverse Jacobian in `SNESQN`.
403
404 Logically Collective
405
406 Input Parameters:
407 + snes - the nonlinear solver context
408 - stype - scale type, see `SNESQNScaleType`
409
410 Options Database Key:
411 . -snes_qn_scale_type <diagonal,none,scalar,jacobian> - Scaling type
412
413 Level: intermediate
414
415 .seealso: [](ch_snes), `SNES`, `SNESQN`, `SNESLineSearch`, `SNESQNScaleType`, `SNESSetJacobian()`, `SNESQNType`, `SNESQNRestartType`
416 @*/
SNESQNSetScaleType(SNES snes,SNESQNScaleType stype)417 PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
418 {
419 PetscFunctionBegin;
420 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
421 PetscTryMethod(snes, "SNESQNSetScaleType_C", (SNES, SNESQNScaleType), (snes, stype));
422 PetscFunctionReturn(PETSC_SUCCESS);
423 }
424
SNESQNSetScaleType_QN(SNES snes,SNESQNScaleType stype)425 static PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
426 {
427 SNES_QN *qn = (SNES_QN *)snes->data;
428
429 PetscFunctionBegin;
430 qn->scale_type = stype;
431 if (stype == SNES_QN_SCALE_JACOBIAN) snes->usesksp = PETSC_TRUE;
432 PetscFunctionReturn(PETSC_SUCCESS);
433 }
434
SNESQNSetRestartType_QN(SNES snes,SNESQNRestartType rtype)435 static PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
436 {
437 SNES_QN *qn = (SNES_QN *)snes->data;
438
439 PetscFunctionBegin;
440 qn->restart_type = rtype;
441 PetscFunctionReturn(PETSC_SUCCESS);
442 }
443
444 /*@
445 SNESQNSetType - Sets the quasi-Newton variant to be used in `SNESQN`.
446
447 Logically Collective
448
449 Input Parameters:
450 + snes - the iterative context
451 - qtype - variant type, see `SNESQNType`
452
453 Options Database Key:
454 . -snes_qn_type <lbfgs,broyden,badbroyden> - quasi-Newton type
455
456 Level: intermediate
457
458 .seealso: [](ch_snes), `SNESQN`, `SNES_QN_LBFGS`, `SNES_QN_BROYDEN`, `SNES_QN_BADBROYDEN`, `SNESQNType`, `SNESQNScaleType`, `TAOLMVM`, `TAOBLMVM`
459 @*/
SNESQNSetType(SNES snes,SNESQNType qtype)460 PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
461 {
462 PetscFunctionBegin;
463 PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
464 PetscTryMethod(snes, "SNESQNSetType_C", (SNES, SNESQNType), (snes, qtype));
465 PetscFunctionReturn(PETSC_SUCCESS);
466 }
467
SNESQNSetType_QN(SNES snes,SNESQNType qtype)468 static PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
469 {
470 SNES_QN *qn = (SNES_QN *)snes->data;
471
472 PetscFunctionBegin;
473 qn->type = qtype;
474 PetscFunctionReturn(PETSC_SUCCESS);
475 }
476
477 /*MC
478 SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.
479
480 Options Database Keys:
481 + -snes_qn_m <m> - Number of past states saved for the L-Broyden methods.
482 . -snes_qn_restart_type <powell,periodic,none> - set the restart type
483 . -snes_qn_powell_gamma - Angle condition for restart.
484 . -snes_qn_powell_descent - Descent condition for restart.
485 . -snes_qn_type <lbfgs,broyden,badbroyden> - QN type
486 . -snes_qn_scale_type <diagonal,none,scalar,jacobian> - scaling performed on inner Jacobian
487 . -snes_linesearch_type <cp, l2, basic> - Type of line search.
488 - -snes_qn_monitor - Monitors the quasi-newton Jacobian.
489
490 Level: beginner
491
492 Notes:
493 This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of $F(x) = b$ using
494 previous change in $F(x)$ and $x$ to form the approximate inverse Jacobian using a series of multiplicative rank-one
495 updates.
496
497 When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied. The first of
498 these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this
499 iteration as the current iteration's values when constructing the approximate Jacobian. The second, composed,
500 perturbs the problem the Jacobian represents to be $P(x, b) - x = 0 $, where $P(x, b)$ is the preconditioner.
501
502 Uses left nonlinear preconditioning by default.
503
504 See {cite}`byrd1994representations`, {cite}`kelley95`, {cite}`brown1985experiments`, {cite}`griewank2012broyden`, {cite}`gilbert1989some`,
505 {cite}`dener2019accelerating`, and {cite}`bruneknepleysmithtu15`
506
507 .seealso: [](ch_snes), `SNESQNRestartType`, `SNESQNSetRestartType()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`,
508 `SNESQNScaleType`, `SNESQNSetScaleType()`, `SNESQNType`, `SNESQNSetType()`
509 M*/
SNESCreate_QN(SNES snes)510 PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
511 {
512 SNES_QN *qn;
513
514 PetscFunctionBegin;
515 snes->ops->setup = SNESSetUp_QN;
516 snes->ops->solve = SNESSolve_QN;
517 snes->ops->destroy = SNESDestroy_QN;
518 snes->ops->setfromoptions = SNESSetFromOptions_QN;
519 snes->ops->view = SNESView_QN;
520 snes->ops->reset = SNESReset_QN;
521
522 snes->npcside = PC_LEFT;
523
524 snes->usesnpc = PETSC_TRUE;
525 snes->usesksp = PETSC_FALSE;
526
527 snes->alwayscomputesfinalresidual = PETSC_TRUE;
528
529 PetscCall(SNESParametersInitialize(snes));
530 PetscObjectParameterSetDefault(snes, max_funcs, 30000);
531 PetscObjectParameterSetDefault(snes, max_its, 10000);
532
533 PetscCall(PetscNew(&qn));
534 snes->data = (void *)qn;
535 qn->m = 10;
536 qn->scaling = 1.0;
537 qn->monitor = NULL;
538 qn->monflg = PETSC_FALSE;
539 qn->powell_gamma = 0.9999;
540 qn->scale_type = SNES_QN_SCALE_DEFAULT;
541 qn->restart_type = SNES_QN_RESTART_DEFAULT;
542 qn->type = SNES_QN_LBFGS;
543
544 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", SNESQNSetScaleType_QN));
545 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", SNESQNSetRestartType_QN));
546 PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", SNESQNSetType_QN));
547 PetscFunctionReturn(PETSC_SUCCESS);
548 }
549