Actual source code: qn.c

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

  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};

  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;

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

 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: }

 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 lssucceed;
 70:   PetscBool            badstep, powell, periodic, restart;
 71:   PetscScalar          DolddotD, DolddotDold;
 72:   SNESConvergedReason  reason;
 73: #if defined(PETSC_USE_INFO)
 74:   PetscReal gnorm;
 75: #endif

 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);

 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 */

 85:   /* directions */
 86:   W    = snes->work[0];
 87:   D    = snes->work[1];
 88:   Dold = snes->work[2];

 90:   snes->reason = SNES_CONVERGED_ITERATING;

 92:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
 93:   snes->iter = 0;
 94:   snes->norm = 0.;
 95:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

 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:   SNESCheckFunctionNorm(snes, fnorm);

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));

119:   /* test convergence */
120:   PetscCall(SNESConverged(snes, 0, 0.0, 0.0, fnorm));
121:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

123:   restart = PETSC_TRUE;
124:   for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
125:     PetscScalar gdx;

127:     /* general purpose update */
128:     PetscTryTypeMethod(snes, update, snes->iter);

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:     }

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:     }

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:     }

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 == SNES_DIVERGED_FUNCTION_COUNT) break;
183:     PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed));
184:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm));
185:     badstep = PETSC_FALSE;
186:     if (lssucceed) {
187:       if (++snes->numFailures >= snes->maxFailures) {
188:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
189:         break;
190:       }
191:       badstep = PETSC_TRUE;
192:     }

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

197:     PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
198:     snes->iter  = i + 1;
199:     snes->norm  = fnorm;
200:     snes->xnorm = xnorm;
201:     snes->ynorm = ynorm;
202:     PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

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

206:     /* set parameter for default relative tolerance convergence test */
207:     PetscCall(SNESConverged(snes, snes->iter, xnorm, ynorm, fnorm));
208:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));
209:     if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

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

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

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

262: static PetscErrorCode SNESSetUp_QN(SNES snes)
263: {
264:   SNES_QN *qn = (SNES_QN *)snes->data;
265:   DM       dm;
266:   PetscInt n, N;

268:   PetscFunctionBegin;
269:   if (!snes->vec_sol) {
270:     PetscCall(SNESGetDM(snes, &dm));
271:     PetscCall(DMCreateGlobalVector(dm, &snes->vec_sol));
272:   }
273:   PetscCall(SNESSetWorkVecs(snes, 3));
274:   PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));

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

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

306: static PetscErrorCode SNESReset_QN(SNES snes)
307: {
308:   SNES_QN *qn = (SNES_QN *)snes->data;

310:   PetscFunctionBegin;
311:   PetscCall(MatDestroy(&qn->B));
312:   PetscFunctionReturn(PETSC_SUCCESS);
313: }

315: static PetscErrorCode SNESDestroy_QN(SNES snes)
316: {
317:   PetscFunctionBegin;
318:   PetscCall(SNESReset_QN(snes));
319:   PetscCall(PetscFree(snes->data));
320:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", NULL));
321:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", NULL));
322:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", NULL));
323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }

326: static PetscErrorCode SNESSetFromOptions_QN(SNES snes, PetscOptionItems PetscOptionsObject)
327: {
328:   SNES_QN          *qn = (SNES_QN *)snes->data;
329:   PetscBool         flg;
330:   SNESLineSearch    linesearch;
331:   SNESQNRestartType rtype = qn->restart_type;
332:   SNESQNScaleType   stype = qn->scale_type;
333:   SNESQNType        qtype = qn->type;

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

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

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

367: static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
368: {
369:   SNES_QN  *qn = (SNES_QN *)snes->data;
370:   PetscBool isascii;

372:   PetscFunctionBegin;
373:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
374:   if (isascii) {
375:     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]));
376:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Stored subspace size: %" PetscInt_FMT "\n", qn->m));
377:     PetscCall(MatView(qn->B, viewer));
378:   }
379:   PetscFunctionReturn(PETSC_SUCCESS);
380: }

382: /*@
383:   SNESQNSetRestartType - Sets the restart type for `SNESQN`.

385:   Logically Collective

387:   Input Parameters:
388: + snes  - the iterative context
389: - rtype - restart type, see `SNESQNRestartType`

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

395:   Level: intermediate

397: .seealso: [](ch_snes), `SNES`, `SNESQN`, `SNESQNRestartType`, `SNES_QN_RESTART_NONE`, `SNES_QN_RESTART_POWELL`, `SNES_QN_RESTART_PERIODIC`,
398:           `SNESQNType`, `SNESQNScaleType`
399: @*/
400: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
401: {
402:   PetscFunctionBegin;
404:   PetscTryMethod(snes, "SNESQNSetRestartType_C", (SNES, SNESQNRestartType), (snes, rtype));
405:   PetscFunctionReturn(PETSC_SUCCESS);
406: }

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

411:   Logically Collective

413:   Input Parameters:
414: + snes  - the nonlinear solver context
415: - stype - scale type, see `SNESQNScaleType`

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

420:   Level: intermediate

422: .seealso: [](ch_snes), `SNES`, `SNESQN`, `SNESLineSearch`, `SNESQNScaleType`, `SNESSetJacobian()`, `SNESQNType`, `SNESQNRestartType`
423: @*/
424: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
425: {
426:   PetscFunctionBegin;
428:   PetscTryMethod(snes, "SNESQNSetScaleType_C", (SNES, SNESQNScaleType), (snes, stype));
429:   PetscFunctionReturn(PETSC_SUCCESS);
430: }

432: static PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
433: {
434:   SNES_QN *qn = (SNES_QN *)snes->data;

436:   PetscFunctionBegin;
437:   qn->scale_type = stype;
438:   if (stype == SNES_QN_SCALE_JACOBIAN) snes->usesksp = PETSC_TRUE;
439:   PetscFunctionReturn(PETSC_SUCCESS);
440: }

442: static PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
443: {
444:   SNES_QN *qn = (SNES_QN *)snes->data;

446:   PetscFunctionBegin;
447:   qn->restart_type = rtype;
448:   PetscFunctionReturn(PETSC_SUCCESS);
449: }

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

454:   Logically Collective

456:   Input Parameters:
457: + snes  - the iterative context
458: - qtype - variant type, see `SNESQNType`

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

463:   Level: intermediate

465: .seealso: [](ch_snes), `SNESQN`, `SNES_QN_LBFGS`, `SNES_QN_BROYDEN`, `SNES_QN_BADBROYDEN`, `SNESQNType`,  `SNESQNScaleType`, `TAOLMVM`, `TAOBLMVM`
466: @*/
467: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
468: {
469:   PetscFunctionBegin;
471:   PetscTryMethod(snes, "SNESQNSetType_C", (SNES, SNESQNType), (snes, qtype));
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: static PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
476: {
477:   SNES_QN *qn = (SNES_QN *)snes->data;

479:   PetscFunctionBegin;
480:   qn->type = qtype;
481:   PetscFunctionReturn(PETSC_SUCCESS);
482: }

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

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

497:       Level: beginner

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

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

509:       Uses left nonlinear preconditioning by default.

511:       See {cite}`byrd1994representations`, {cite}`kelley95`, {cite}`brown1985experiments`, {cite}`griewank2012broyden`, {cite}`gilbert1989some`,
512:       {cite}`dener2019accelerating`, and {cite}`bruneknepleysmithtu15`

514: .seealso: [](ch_snes), `SNESQNRestartType`, `SNESQNSetRestartType()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`,
515:           `SNESQNScaleType`, `SNESQNSetScaleType()`, `SNESQNType`, `SNESQNSetType()`
516: M*/
517: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
518: {
519:   SNES_QN *qn;

521:   PetscFunctionBegin;
522:   snes->ops->setup          = SNESSetUp_QN;
523:   snes->ops->solve          = SNESSolve_QN;
524:   snes->ops->destroy        = SNESDestroy_QN;
525:   snes->ops->setfromoptions = SNESSetFromOptions_QN;
526:   snes->ops->view           = SNESView_QN;
527:   snes->ops->reset          = SNESReset_QN;

529:   snes->npcside = PC_LEFT;

531:   snes->usesnpc = PETSC_TRUE;
532:   snes->usesksp = PETSC_FALSE;

534:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

536:   PetscCall(SNESParametersInitialize(snes));
537:   PetscObjectParameterSetDefault(snes, max_funcs, 30000);
538:   PetscObjectParameterSetDefault(snes, max_its, 10000);

540:   PetscCall(PetscNew(&qn));
541:   snes->data       = (void *)qn;
542:   qn->m            = 10;
543:   qn->scaling      = 1.0;
544:   qn->monitor      = NULL;
545:   qn->monflg       = PETSC_FALSE;
546:   qn->powell_gamma = 0.9999;
547:   qn->scale_type   = SNES_QN_SCALE_DEFAULT;
548:   qn->restart_type = SNES_QN_RESTART_DEFAULT;
549:   qn->type         = SNES_QN_LBFGS;

551:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", SNESQNSetScaleType_QN));
552:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", SNESQNSetRestartType_QN));
553:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", SNESQNSetType_QN));
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }