Actual source code: qn.c

  1: #include <petsc/private/snesimpl.h>
  2: #include <petscdm.h>

  4: #define H(i, j) qn->dXdFmat[i * qn->m + j]

  6: const char *const SNESQNScaleTypes[]   = {"DEFAULT", "NONE", "SCALAR", "DIAGONAL", "JACOBIAN", "SNESQNScaleType", "SNES_QN_SCALING_", NULL};
  7: const char *const SNESQNRestartTypes[] = {"DEFAULT", "NONE", "POWELL", "PERIODIC", "SNESQNRestartType", "SNES_QN_RESTART_", NULL};
  8: const char *const SNESQNTypes[]        = {"LBFGS", "BROYDEN", "BADBROYDEN", "SNESQNType", "SNES_QN_", NULL};

 10: typedef struct {
 11:   Mat               B;      /* Quasi-Newton approximation Matrix (MATLMVM) */
 12:   PetscInt          m;      /* The number of kept previous steps */
 13:   PetscReal        *lambda; /* The line search history of the method */
 14:   PetscBool         monflg;
 15:   PetscViewer       monitor;
 16:   PetscReal         powell_gamma; /* Powell angle restart condition */
 17:   PetscReal         scaling;      /* scaling of H0 */
 18:   SNESQNType        type;         /* the type of quasi-newton method used */
 19:   SNESQNScaleType   scale_type;   /* the type of scaling used */
 20:   SNESQNRestartType restart_type; /* determine the frequency and type of restart conditions */
 21: } SNES_QN;

 23: static PetscErrorCode SNESQNGetMatrix_Private(SNES snes, Mat *B)
 24: {
 25:   SNES_QN    *qn = (SNES_QN *)snes->data;
 26:   const char *optionsprefix;

 28:   PetscFunctionBegin;
 29:   if (!qn->B) {
 30:     PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &qn->B));
 31:     PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
 32:     PetscCall(MatSetOptionsPrefix(qn->B, optionsprefix));
 33:     PetscCall(MatAppendOptionsPrefix(qn->B, "qn_"));
 34:     switch (qn->type) {
 35:     case SNES_QN_BROYDEN:
 36:       PetscCall(MatSetType(qn->B, MATLMVMBROYDEN));
 37:       qn->scale_type = SNES_QN_SCALE_NONE;
 38:       break;
 39:     case SNES_QN_BADBROYDEN:
 40:       PetscCall(MatSetType(qn->B, MATLMVMBADBROYDEN));
 41:       qn->scale_type = SNES_QN_SCALE_NONE;
 42:       break;
 43:     default:
 44:       PetscCall(MatSetType(qn->B, MATLMVMBFGS));
 45:       switch (qn->scale_type) {
 46:       case SNES_QN_SCALE_NONE:
 47:         PetscCall(MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_NONE));
 48:         break;
 49:       case SNES_QN_SCALE_SCALAR:
 50:       case SNES_QN_SCALE_DEFAULT:
 51:         PetscCall(MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_SCALAR));
 52:         break;
 53:       case SNES_QN_SCALE_JACOBIAN:
 54:         PetscCall(MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_USER));
 55:         break;
 56:       case SNES_QN_SCALE_DIAGONAL:
 57:       default:
 58:         break;
 59:       }
 60:       break;
 61:     }
 62:   }
 63:   *B = qn->B;
 64:   PetscFunctionReturn(PETSC_SUCCESS);
 65: }

 67: static PetscErrorCode SNESSolve_QN(SNES snes)
 68: {
 69:   SNES_QN             *qn = (SNES_QN *)snes->data;
 70:   Vec                  X, F, W;
 71:   Vec                  Y, D, Dold;
 72:   PetscInt             i, i_r;
 73:   PetscReal            fnorm, xnorm, ynorm;
 74:   SNESLineSearchReason lssucceed;
 75:   PetscBool            badstep, powell, periodic, restart;
 76:   PetscScalar          DolddotD, DolddotDold;
 77:   SNESConvergedReason  reason;
 78: #if defined(PETSC_USE_INFO)
 79:   PetscReal gnorm;
 80: #endif

 82:   PetscFunctionBegin;
 83:   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);

 85:   PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
 86:   F = snes->vec_func;       /* residual vector */
 87:   Y = snes->vec_sol_update; /* search direction generated by J^-1D*/
 88:   X = snes->vec_sol;        /* solution vector */

 90:   /* directions */
 91:   W    = snes->work[0];
 92:   D    = snes->work[1];
 93:   Dold = snes->work[2];

 95:   snes->reason = SNES_CONVERGED_ITERATING;

 97:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
 98:   snes->iter = 0;
 99:   snes->norm = 0.;
100:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

102:   if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
103:     /* we need to sample the preconditioned function only once here,
104:        since linesearch will always use the preconditioned function */
105:     PetscCall(SNESApplyNPC(snes, X, NULL, F));
106:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
107:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
108:       snes->reason = SNES_DIVERGED_INNER;
109:       PetscFunctionReturn(PETSC_SUCCESS);
110:     }
111:   } else {
112:     if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
113:     else snes->vec_func_init_set = PETSC_FALSE;
114:   }
115:   PetscCall(VecNorm(F, NORM_2, &fnorm));
116:   SNESCheckFunctionNorm(snes, fnorm);

118:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
119:   snes->norm = fnorm;
120:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
121:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
122:   PetscCall(SNESMonitor(snes, 0, fnorm));

124:   /* test convergence */
125:   PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
126:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

128:   restart = PETSC_TRUE;
129:   for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
130:     PetscScalar gdx;

132:     /* general purpose update */
133:     PetscTryTypeMethod(snes, update, snes->iter);

135:     if (snes->npc) {
136:       if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
137:         PetscCall(SNESApplyNPC(snes, X, F, D));
138:         PetscCall(SNESGetConvergedReason(snes->npc, &reason));
139:         if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
140:           snes->reason = SNES_DIVERGED_INNER;
141:           PetscFunctionReturn(PETSC_SUCCESS);
142:         }
143:       } else if (snes->npcside == PC_RIGHT) {
144:         PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, 0, 0));
145:         PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
146:         PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, 0, 0));
147:         PetscCall(SNESGetConvergedReason(snes->npc, &reason));
148:         if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
149:           snes->reason = SNES_DIVERGED_INNER;
150:           PetscFunctionReturn(PETSC_SUCCESS);
151:         }
152:         PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
153:         PetscCall(VecCopy(F, D));
154:       } else {
155:         PetscCall(VecCopy(F, D));
156:       }
157:     } else {
158:       PetscCall(VecCopy(F, D));
159:     }

161:     /* scale the initial update */
162:     if (qn->scale_type == SNES_QN_SCALE_JACOBIAN && restart) {
163:       PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
164:       SNESCheckJacobianDomainerror(snes);
165:       PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
166:       PetscCall(MatLMVMSetJ0KSP(qn->B, snes->ksp));
167:     }

169:     /* update QN approx and calculate step */
170:     PetscCall(MatLMVMUpdate(qn->B, X, D));
171:     PetscCall(MatSolve(qn->B, D, Y));
172:     PetscCall(VecDot(F, Y, &gdx));
173:     if (PetscAbsScalar(gdx) <= 0) {
174:       /* Step is not descent or solve was not successful
175:          Use steepest descent direction */
176:       PetscCall(VecCopy(F, Y));
177:       PetscCall(MatLMVMReset(qn->B, PETSC_FALSE));
178:     }

180:     /* line search for lambda */
181:     ynorm = 1;
182: #if defined(PETSC_USE_INFO)
183:     gnorm = fnorm;
184: #endif
185:     PetscCall(VecCopy(D, Dold));
186:     PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y));
187:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
188:     PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed));
189:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm));
190:     badstep = PETSC_FALSE;
191:     if (lssucceed) {
192:       if (++snes->numFailures >= snes->maxFailures) {
193:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
194:         break;
195:       }
196:       badstep = PETSC_TRUE;
197:     }

199:     /* convergence monitoring */
200:     PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lssucceed));

202:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
203:     snes->iter  = i + 1;
204:     snes->norm  = fnorm;
205:     snes->xnorm = xnorm;
206:     snes->ynorm = ynorm;
207:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

209:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, snes->iter));

211:     /* set parameter for default relative tolerance convergence test */
212:     PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
213:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
214:     if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

216:     if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
217:       PetscCall(SNESApplyNPC(snes, X, F, D));
218:       PetscCall(SNESGetConvergedReason(snes->npc, &reason));
219:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
220:         snes->reason = SNES_DIVERGED_INNER;
221:         PetscFunctionReturn(PETSC_SUCCESS);
222:       }
223:     } else {
224:       PetscCall(VecCopy(F, D));
225:     }

227:     /* restart conditions */
228:     powell = PETSC_FALSE;
229:     if (qn->restart_type == SNES_QN_RESTART_POWELL && i_r > 1) {
230:       /* check restart by Powell's Criterion: |F^T H_0 Fold| > powell_gamma * |Fold^T H_0 Fold| */
231:       if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
232:         PetscCall(MatMult(snes->jacobian, Dold, W));
233:       } else {
234:         PetscCall(VecCopy(Dold, W));
235:       }
236:       PetscCall(VecDotBegin(W, Dold, &DolddotDold));
237:       PetscCall(VecDotBegin(W, D, &DolddotD));
238:       PetscCall(VecDotEnd(W, Dold, &DolddotDold));
239:       PetscCall(VecDotEnd(W, D, &DolddotD));
240:       if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma * PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE;
241:     }
242:     periodic = PETSC_FALSE;
243:     if (qn->restart_type == SNES_QN_RESTART_PERIODIC) {
244:       if (i_r > qn->m - 1) periodic = PETSC_TRUE;
245:     }

247:     /* restart if either powell or periodic restart is satisfied. */
248:     restart = PETSC_FALSE;
249:     if (badstep || powell || periodic) {
250:       restart = PETSC_TRUE;
251:       if (qn->monflg) {
252:         PetscCall(PetscViewerASCIIAddTab(qn->monitor, ((PetscObject)snes)->tablevel + 2));
253:         if (powell) {
254:           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));
255:         } else {
256:           PetscCall(PetscViewerASCIIPrintf(qn->monitor, "Periodic restart! i_r = %" PetscInt_FMT "\n", i_r));
257:         }
258:         PetscCall(PetscViewerASCIISubtractTab(qn->monitor, ((PetscObject)snes)->tablevel + 2));
259:       }
260:       i_r = -1;
261:       PetscCall(MatLMVMReset(qn->B, PETSC_FALSE));
262:     }
263:   }
264:   PetscFunctionReturn(PETSC_SUCCESS);
265: }

267: static PetscErrorCode SNESSetUp_QN(SNES snes)
268: {
269:   SNES_QN *qn = (SNES_QN *)snes->data;
270:   DM       dm;
271:   PetscInt n, N;

273:   PetscFunctionBegin;
274:   if (!snes->vec_sol) {
275:     PetscCall(SNESGetDM(snes, &dm));
276:     PetscCall(DMCreateGlobalVector(dm, &snes->vec_sol));
277:   }
278:   PetscCall(SNESSetWorkVecs(snes, 3));
279:   PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));

281:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) PetscCall(SNESSetUpMatrices(snes));
282:   if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_PRECONDITIONED;

284:   /* set method defaults */
285:   if (qn->scale_type == SNES_QN_SCALE_DEFAULT) {
286:     if (qn->type == SNES_QN_BADBROYDEN) {
287:       qn->scale_type = SNES_QN_SCALE_NONE;
288:     } else {
289:       qn->scale_type = SNES_QN_SCALE_SCALAR;
290:     }
291:   }
292:   if (qn->restart_type == SNES_QN_RESTART_DEFAULT) {
293:     if (qn->type == SNES_QN_LBFGS) {
294:       qn->restart_type = SNES_QN_RESTART_POWELL;
295:     } else {
296:       qn->restart_type = SNES_QN_RESTART_PERIODIC;
297:     }
298:   }
299:   /* Set up the LMVM matrix */
300:   PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));
301:   PetscCall(VecGetLocalSize(snes->vec_sol, &n));
302:   PetscCall(VecGetSize(snes->vec_sol, &N));
303:   PetscCall(MatSetSizes(qn->B, n, n, N, N));
304:   PetscCall(MatSetUp(qn->B));
305:   PetscCall(MatLMVMReset(qn->B, PETSC_TRUE));
306:   PetscCall(MatLMVMSetHistorySize(qn->B, qn->m));
307:   PetscCall(MatLMVMAllocate(qn->B, snes->vec_sol, snes->vec_func));
308:   PetscFunctionReturn(PETSC_SUCCESS);
309: }

311: static PetscErrorCode SNESReset_QN(SNES snes)
312: {
313:   SNES_QN *qn = (SNES_QN *)snes->data;

315:   PetscFunctionBegin;
316:   PetscCall(MatDestroy(&qn->B));
317:   PetscFunctionReturn(PETSC_SUCCESS);
318: }

320: static PetscErrorCode SNESDestroy_QN(SNES snes)
321: {
322:   PetscFunctionBegin;
323:   PetscCall(SNESReset_QN(snes));
324:   PetscCall(PetscFree(snes->data));
325:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", NULL));
326:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", NULL));
327:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", NULL));
328:   PetscFunctionReturn(PETSC_SUCCESS);
329: }

331: static PetscErrorCode SNESSetFromOptions_QN(SNES snes, PetscOptionItems *PetscOptionsObject)
332: {
333:   SNES_QN          *qn = (SNES_QN *)snes->data;
334:   PetscBool         flg;
335:   SNESLineSearch    linesearch;
336:   SNESQNRestartType rtype = qn->restart_type;
337:   SNESQNScaleType   stype = qn->scale_type;
338:   SNESQNType        qtype = qn->type;

340:   PetscFunctionBegin;
341:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES QN options");
342:   PetscCall(PetscOptionsInt("-snes_qn_m", "Number of past states saved for L-BFGS methods", "SNESQN", qn->m, &qn->m, NULL));
343:   PetscCall(PetscOptionsReal("-snes_qn_powell_gamma", "Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL));
344:   PetscCall(PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", qn->monflg, &qn->monflg, NULL));
345:   PetscCall(PetscOptionsEnum("-snes_qn_scale_type", "Scaling type", "SNESQNSetScaleType", SNESQNScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg));
346:   if (flg) PetscCall(SNESQNSetScaleType(snes, stype));

348:   PetscCall(PetscOptionsEnum("-snes_qn_restart_type", "Restart type", "SNESQNSetRestartType", SNESQNRestartTypes, (PetscEnum)rtype, (PetscEnum *)&rtype, &flg));
349:   if (flg) PetscCall(SNESQNSetRestartType(snes, rtype));

351:   PetscCall(PetscOptionsEnum("-snes_qn_type", "Quasi-Newton update type", "", SNESQNTypes, (PetscEnum)qtype, (PetscEnum *)&qtype, &flg));
352:   if (flg) PetscCall(SNESQNSetType(snes, qtype));
353:   PetscOptionsHeadEnd();
354:   PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));
355:   PetscCall(MatSetFromOptions(qn->B));
356:   if (!snes->linesearch) {
357:     PetscCall(SNESGetLineSearch(snes, &linesearch));
358:     if (!((PetscObject)linesearch)->type_name) {
359:       if (qn->type == SNES_QN_LBFGS) {
360:         PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHCP));
361:       } else if (qn->type == SNES_QN_BROYDEN) {
362:         PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC));
363:       } else {
364:         PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHL2));
365:       }
366:     }
367:   }
368:   if (qn->monflg) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &qn->monitor));
369:   PetscFunctionReturn(PETSC_SUCCESS);
370: }

372: static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
373: {
374:   SNES_QN  *qn = (SNES_QN *)snes->data;
375:   PetscBool iascii;

377:   PetscFunctionBegin;
378:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
379:   if (iascii) {
380:     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]));
381:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Stored subspace size: %" PetscInt_FMT "\n", qn->m));
382:     PetscCall(MatView(qn->B, viewer));
383:   }
384:   PetscFunctionReturn(PETSC_SUCCESS);
385: }

387: /*@
388:   SNESQNSetRestartType - Sets the restart type for `SNESQN`.

390:   Logically Collective

392:   Input Parameters:
393: + snes  - the iterative context
394: - rtype - restart type, see `SNESQNRestartType`

396:   Options Database Keys:
397: + -snes_qn_restart_type <powell,periodic,none> - set the restart type
398: - -snes_qn_m <m>                               - sets the number of stored updates and the restart period for periodic

400:   Level: intermediate

402: .seealso: [](ch_snes), `SNES`, `SNESQN`, `SNESQNRestartType`, `SNES_QN_RESTART_NONE`, `SNES_QN_RESTART_POWELL`, `SNES_QN_RESTART_PERIODIC`,
403:           `SNESQNType`, `SNESQNScaleType`
404: @*/
405: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
406: {
407:   PetscFunctionBegin;
409:   PetscTryMethod(snes, "SNESQNSetRestartType_C", (SNES, SNESQNRestartType), (snes, rtype));
410:   PetscFunctionReturn(PETSC_SUCCESS);
411: }

413: /*@
414:   SNESQNSetScaleType - Sets the scaling type for the inner inverse Jacobian in `SNESQN`.

416:   Logically Collective

418:   Input Parameters:
419: + snes  - the nonlinear solver context
420: - stype - scale type, see `SNESQNScaleType`

422:   Options Database Key:
423: . -snes_qn_scale_type <diagonal,none,scalar,jacobian> - Scaling type

425:   Level: intermediate

427: .seealso: [](ch_snes), `SNES`, `SNESQN`, `SNESLineSearch`, `SNESQNScaleType`, `SNESSetJacobian()`, `SNESQNType`, `SNESQNRestartType`
428: @*/
429: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
430: {
431:   PetscFunctionBegin;
433:   PetscTryMethod(snes, "SNESQNSetScaleType_C", (SNES, SNESQNScaleType), (snes, stype));
434:   PetscFunctionReturn(PETSC_SUCCESS);
435: }

437: static PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
438: {
439:   SNES_QN *qn = (SNES_QN *)snes->data;

441:   PetscFunctionBegin;
442:   qn->scale_type = stype;
443:   if (stype == SNES_QN_SCALE_JACOBIAN) snes->usesksp = PETSC_TRUE;
444:   PetscFunctionReturn(PETSC_SUCCESS);
445: }

447: static PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
448: {
449:   SNES_QN *qn = (SNES_QN *)snes->data;

451:   PetscFunctionBegin;
452:   qn->restart_type = rtype;
453:   PetscFunctionReturn(PETSC_SUCCESS);
454: }

456: /*@
457:   SNESQNSetType - Sets the quasi-Newton variant to be used in `SNESQN`.

459:   Logically Collective

461:   Input Parameters:
462: + snes  - the iterative context
463: - qtype - variant type, see `SNESQNType`

465:   Options Database Key:
466: . -snes_qn_type <lbfgs,broyden,badbroyden> - quasi-Newton type

468:   Level: intermediate

470: .seealso: [](ch_snes), `SNESQN`, `SNES_QN_LBFGS`, `SNES_QN_BROYDEN`, `SNES_QN_BADBROYDEN`, `SNESQNType`,  `SNESQNScaleType`, `TAOLMVM`, `TAOBLMVM`
471: @*/
472: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
473: {
474:   PetscFunctionBegin;
476:   PetscTryMethod(snes, "SNESQNSetType_C", (SNES, SNESQNType), (snes, qtype));
477:   PetscFunctionReturn(PETSC_SUCCESS);
478: }

480: static PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
481: {
482:   SNES_QN *qn = (SNES_QN *)snes->data;

484:   PetscFunctionBegin;
485:   qn->type = qtype;
486:   PetscFunctionReturn(PETSC_SUCCESS);
487: }

489: /*MC
490:       SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.

492:       Options Database Keys:
493: +     -snes_qn_m <m>                                      - Number of past states saved for the L-Broyden methods.
494: .     -snes_qn_restart_type <powell,periodic,none>        - set the restart type
495: .     -snes_qn_powell_gamma                               - Angle condition for restart.
496: .     -snes_qn_powell_descent                             - Descent condition for restart.
497: .     -snes_qn_type <lbfgs,broyden,badbroyden>            - QN type
498: .     -snes_qn_scale_type <diagonal,none,scalar,jacobian> - scaling performed on inner Jacobian
499: .     -snes_linesearch_type <cp, l2, basic>               - Type of line search.
500: -     -snes_qn_monitor                                    - Monitors the quasi-newton Jacobian.

502:       Level: beginner

504:       Notes:
505:       This implements the L-BFGS, Broyden, and "Bad" Broyden algorithms for the solution of $F(x) = b$ using
506:       previous change in $F(x)$ and $x$ to form the approximate inverse Jacobian using a series of multiplicative rank-one
507:       updates.

509:       When using a nonlinear preconditioner, one has two options as to how the preconditioner is applied.  The first of
510:       these options, sequential, uses the preconditioner to generate a new solution and function and uses those at this
511:       iteration as the current iteration's values when constructing the approximate Jacobian.  The second, composed,
512:       perturbs the problem the Jacobian represents to be $P(x, b) - x = 0 $, where $P(x, b)$ is the preconditioner.

514:       Uses left nonlinear preconditioning by default.

516:       See {cite}`byrd1994representations`, {cite}`kelley95`, {cite}`brown1985experiments`, {cite}`griewank2012broyden`, {cite}`gilbert1989some`,
517:       {cite}`dener2019accelerating`, and {cite}`bruneknepleysmithtu15`

519: .seealso: [](ch_snes), `SNESQNRestartType`, `SNESQNSetRestartType()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`,
520:           `SNESQNScaleType`, `SNESQNSetScaleType()`, `SNESQNType`, `SNESQNSetType()`
521: M*/
522: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
523: {
524:   SNES_QN *qn;

526:   PetscFunctionBegin;
527:   snes->ops->setup          = SNESSetUp_QN;
528:   snes->ops->solve          = SNESSolve_QN;
529:   snes->ops->destroy        = SNESDestroy_QN;
530:   snes->ops->setfromoptions = SNESSetFromOptions_QN;
531:   snes->ops->view           = SNESView_QN;
532:   snes->ops->reset          = SNESReset_QN;

534:   snes->npcside = PC_LEFT;

536:   snes->usesnpc = PETSC_TRUE;
537:   snes->usesksp = PETSC_FALSE;

539:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

541:   PetscCall(SNESParametersInitialize(snes));
542:   PetscObjectParameterSetDefault(snes, max_funcs, 30000);
543:   PetscObjectParameterSetDefault(snes, max_its, 10000);

545:   PetscCall(PetscNew(&qn));
546:   snes->data       = (void *)qn;
547:   qn->m            = 10;
548:   qn->scaling      = 1.0;
549:   qn->monitor      = NULL;
550:   qn->monflg       = PETSC_FALSE;
551:   qn->powell_gamma = 0.9999;
552:   qn->scale_type   = SNES_QN_SCALE_DEFAULT;
553:   qn->restart_type = SNES_QN_RESTART_DEFAULT;
554:   qn->type         = SNES_QN_LBFGS;

556:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", SNESQNSetScaleType_QN));
557:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", SNESQNSetRestartType_QN));
558:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", SNESQNSetType_QN));
559:   PetscFunctionReturn(PETSC_SUCCESS);
560: }