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 SNESSolve_QN(SNES snes)
 24: {
 25:   SNES_QN             *qn = (SNES_QN *)snes->data;
 26:   Vec                  X, Xold;
 27:   Vec                  F, W;
 28:   Vec                  Y, D, Dold;
 29:   PetscInt             i, i_r;
 30:   PetscReal            fnorm, xnorm, ynorm, gnorm;
 31:   SNESLineSearchReason lssucceed;
 32:   PetscBool            badstep, powell, periodic;
 33:   PetscScalar          DolddotD, DolddotDold;
 34:   SNESConvergedReason  reason;

 36:   /* basically just a regular newton's method except for the application of the Jacobian */


 40:   PetscCitationsRegister(SNESCitation, &SNEScite);
 41:   F    = snes->vec_func;       /* residual vector */
 42:   Y    = snes->vec_sol_update; /* search direction generated by J^-1D*/
 43:   W    = snes->work[3];
 44:   X    = snes->vec_sol; /* solution vector */
 45:   Xold = snes->work[0];

 47:   /* directions generated by the preconditioned problem with F_pre = F or x - M(x, b) */
 48:   D    = snes->work[1];
 49:   Dold = snes->work[2];

 51:   snes->reason = SNES_CONVERGED_ITERATING;

 53:   PetscObjectSAWsTakeAccess((PetscObject)snes);
 54:   snes->iter = 0;
 55:   snes->norm = 0.;
 56:   PetscObjectSAWsGrantAccess((PetscObject)snes);

 58:   if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
 59:     SNESApplyNPC(snes, X, NULL, F);
 60:     SNESGetConvergedReason(snes->npc, &reason);
 61:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 62:       snes->reason = SNES_DIVERGED_INNER;
 63:       return 0;
 64:     }
 65:     VecNorm(F, NORM_2, &fnorm);
 66:   } else {
 67:     if (!snes->vec_func_init_set) {
 68:       SNESComputeFunction(snes, X, F);
 69:     } else snes->vec_func_init_set = PETSC_FALSE;

 71:     VecNorm(F, NORM_2, &fnorm);
 72:     SNESCheckFunctionNorm(snes, fnorm);
 73:   }
 74:   if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
 75:     SNESApplyNPC(snes, X, F, D);
 76:     SNESGetConvergedReason(snes->npc, &reason);
 77:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 78:       snes->reason = SNES_DIVERGED_INNER;
 79:       return 0;
 80:     }
 81:   } else {
 82:     VecCopy(F, D);
 83:   }

 85:   PetscObjectSAWsTakeAccess((PetscObject)snes);
 86:   snes->norm = fnorm;
 87:   PetscObjectSAWsGrantAccess((PetscObject)snes);
 88:   SNESLogConvergenceHistory(snes, fnorm, 0);
 89:   SNESMonitor(snes, 0, fnorm);

 91:   /* test convergence */
 92:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
 93:   if (snes->reason) return 0;

 95:   if (snes->npc && snes->npcside == PC_RIGHT) {
 96:     PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, 0, 0);
 97:     SNESSolve(snes->npc, snes->vec_rhs, X);
 98:     PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, 0, 0);
 99:     SNESGetConvergedReason(snes->npc, &reason);
100:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
101:       snes->reason = SNES_DIVERGED_INNER;
102:       return 0;
103:     }
104:     SNESGetNPCFunction(snes, F, &fnorm);
105:     VecCopy(F, D);
106:   }

108:   /* general purpose update */
109:   PetscTryTypeMethod(snes, update, snes->iter);

111:   /* scale the initial update */
112:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
113:     SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
114:     SNESCheckJacobianDomainerror(snes);
115:     KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre);
116:     MatLMVMSetJ0KSP(qn->B, snes->ksp);
117:   }

119:   for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
120:     /* update QN approx and calculate step */
121:     MatLMVMUpdate(qn->B, X, D);
122:     MatSolve(qn->B, D, Y);

124:     /* line search for lambda */
125:     ynorm = 1;
126:     gnorm = fnorm;
127:     VecCopy(D, Dold);
128:     VecCopy(X, Xold);
129:     SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y);
130:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
131:     SNESLineSearchGetReason(snes->linesearch, &lssucceed);
132:     SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm);
133:     badstep = PETSC_FALSE;
134:     if (lssucceed) {
135:       if (++snes->numFailures >= snes->maxFailures) {
136:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
137:         break;
138:       }
139:       badstep = PETSC_TRUE;
140:     }

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

145:     if (snes->npc && snes->npcside == PC_RIGHT) {
146:       PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, 0, 0);
147:       SNESSolve(snes->npc, snes->vec_rhs, X);
148:       PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, 0, 0);
149:       SNESGetConvergedReason(snes->npc, &reason);
150:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
151:         snes->reason = SNES_DIVERGED_INNER;
152:         return 0;
153:       }
154:       SNESGetNPCFunction(snes, F, &fnorm);
155:     }

157:     SNESSetIterationNumber(snes, i + 1);
158:     snes->norm  = fnorm;
159:     snes->xnorm = xnorm;
160:     snes->ynorm = ynorm;

162:     SNESLogConvergenceHistory(snes, snes->norm, snes->iter);
163:     SNESMonitor(snes, snes->iter, snes->norm);

165:     /* set parameter for default relative tolerance convergence test */
166:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
167:     if (snes->reason) return 0;
168:     if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
169:       SNESApplyNPC(snes, X, F, D);
170:       SNESGetConvergedReason(snes->npc, &reason);
171:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
172:         snes->reason = SNES_DIVERGED_INNER;
173:         return 0;
174:       }
175:     } else {
176:       VecCopy(F, D);
177:     }

179:     /* general purpose update */
180:     PetscTryTypeMethod(snes, update, snes->iter);

182:     /* restart conditions */
183:     powell = PETSC_FALSE;
184:     if (qn->restart_type == SNES_QN_RESTART_POWELL && i_r > 1) {
185:       /* check restart by Powell's Criterion: |F^T H_0 Fold| > powell_gamma * |Fold^T H_0 Fold| */
186:       if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
187:         MatMult(snes->jacobian, Dold, W);
188:       } else {
189:         VecCopy(Dold, W);
190:       }
191:       VecDotBegin(W, Dold, &DolddotDold);
192:       VecDotBegin(W, D, &DolddotD);
193:       VecDotEnd(W, Dold, &DolddotDold);
194:       VecDotEnd(W, D, &DolddotD);
195:       if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma * PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE;
196:     }
197:     periodic = PETSC_FALSE;
198:     if (qn->restart_type == SNES_QN_RESTART_PERIODIC) {
199:       if (i_r > qn->m - 1) periodic = PETSC_TRUE;
200:     }
201:     /* restart if either powell or periodic restart is satisfied. */
202:     if (badstep || powell || periodic) {
203:       if (qn->monflg) {
204:         PetscViewerASCIIAddTab(qn->monitor, ((PetscObject)snes)->tablevel + 2);
205:         if (powell) {
206:           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);
207:         } else {
208:           PetscViewerASCIIPrintf(qn->monitor, "Periodic restart! i_r = %" PetscInt_FMT "\n", i_r);
209:         }
210:         PetscViewerASCIISubtractTab(qn->monitor, ((PetscObject)snes)->tablevel + 2);
211:       }
212:       i_r = -1;
213:       if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
214:         SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre);
215:         SNESCheckJacobianDomainerror(snes);
216:       }
217:       MatLMVMReset(qn->B, PETSC_FALSE);
218:     }
219:   }
220:   if (i == snes->max_its) {
221:     PetscInfo(snes, "Maximum number of iterations has been reached: %" PetscInt_FMT "\n", snes->max_its);
222:     if (!snes->reason) snes->reason = SNES_DIVERGED_MAX_IT;
223:   }
224:   return 0;
225: }

227: static PetscErrorCode SNESSetUp_QN(SNES snes)
228: {
229:   SNES_QN *qn = (SNES_QN *)snes->data;
230:   DM       dm;
231:   PetscInt n, N;


234:   if (!snes->vec_sol) {
235:     SNESGetDM(snes, &dm);
236:     DMCreateGlobalVector(dm, &snes->vec_sol);
237:   }
238:   SNESSetWorkVecs(snes, 4);

240:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) SNESSetUpMatrices(snes);
241:   if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED;

243:   /* set method defaults */
244:   if (qn->scale_type == SNES_QN_SCALE_DEFAULT) {
245:     if (qn->type == SNES_QN_BADBROYDEN) {
246:       qn->scale_type = SNES_QN_SCALE_NONE;
247:     } else {
248:       qn->scale_type = SNES_QN_SCALE_SCALAR;
249:     }
250:   }
251:   if (qn->restart_type == SNES_QN_RESTART_DEFAULT) {
252:     if (qn->type == SNES_QN_LBFGS) {
253:       qn->restart_type = SNES_QN_RESTART_POWELL;
254:     } else {
255:       qn->restart_type = SNES_QN_RESTART_PERIODIC;
256:     }
257:   }

259:   /* Set up the LMVM matrix */
260:   switch (qn->type) {
261:   case SNES_QN_BROYDEN:
262:     MatSetType(qn->B, MATLMVMBROYDEN);
263:     qn->scale_type = SNES_QN_SCALE_NONE;
264:     break;
265:   case SNES_QN_BADBROYDEN:
266:     MatSetType(qn->B, MATLMVMBADBROYDEN);
267:     qn->scale_type = SNES_QN_SCALE_NONE;
268:     break;
269:   default:
270:     MatSetType(qn->B, MATLMVMBFGS);
271:     switch (qn->scale_type) {
272:     case SNES_QN_SCALE_NONE:
273:       MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_NONE);
274:       break;
275:     case SNES_QN_SCALE_SCALAR:
276:       MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_SCALAR);
277:       break;
278:     case SNES_QN_SCALE_JACOBIAN:
279:       MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_USER);
280:       break;
281:     case SNES_QN_SCALE_DIAGONAL:
282:     case SNES_QN_SCALE_DEFAULT:
283:     default:
284:       break;
285:     }
286:     break;
287:   }
288:   VecGetLocalSize(snes->vec_sol, &n);
289:   VecGetSize(snes->vec_sol, &N);
290:   MatSetSizes(qn->B, n, n, N, N);
291:   MatSetUp(qn->B);
292:   MatLMVMReset(qn->B, PETSC_TRUE);
293:   MatLMVMSetHistorySize(qn->B, qn->m);
294:   MatLMVMAllocate(qn->B, snes->vec_sol, snes->vec_func);
295:   return 0;
296: }

298: static PetscErrorCode SNESReset_QN(SNES snes)
299: {
300:   SNES_QN *qn;

302:   if (snes->data) {
303:     qn = (SNES_QN *)snes->data;
304:     MatDestroy(&qn->B);
305:   }
306:   return 0;
307: }

309: static PetscErrorCode SNESDestroy_QN(SNES snes)
310: {
311:   SNESReset_QN(snes);
312:   PetscFree(snes->data);
313:   PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", NULL);
314:   PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", NULL);
315:   PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", NULL);
316:   return 0;
317: }

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;

328:   PetscOptionsHeadBegin(PetscOptionsObject, "SNES QN options");
329:   PetscOptionsInt("-snes_qn_m", "Number of past states saved for L-BFGS methods", "SNESQN", qn->m, &qn->m, NULL);
330:   PetscOptionsReal("-snes_qn_powell_gamma", "Powell angle tolerance", "SNESQN", qn->powell_gamma, &qn->powell_gamma, NULL);
331:   PetscOptionsBool("-snes_qn_monitor", "Monitor for the QN methods", "SNESQN", qn->monflg, &qn->monflg, NULL);
332:   PetscOptionsEnum("-snes_qn_scale_type", "Scaling type", "SNESQNSetScaleType", SNESQNScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg);
333:   if (flg) SNESQNSetScaleType(snes, stype);

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

338:   PetscOptionsEnum("-snes_qn_type", "Quasi-Newton update type", "", SNESQNTypes, (PetscEnum)qtype, (PetscEnum *)&qtype, &flg);
339:   if (flg) SNESQNSetType(snes, qtype);
340:   MatSetFromOptions(qn->B);
341:   PetscOptionsHeadEnd();
342:   if (!snes->linesearch) {
343:     SNESGetLineSearch(snes, &linesearch);
344:     if (!((PetscObject)linesearch)->type_name) {
345:       if (qn->type == SNES_QN_LBFGS) {
346:         SNESLineSearchSetType(linesearch, SNESLINESEARCHCP);
347:       } else if (qn->type == SNES_QN_BROYDEN) {
348:         SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC);
349:       } else {
350:         SNESLineSearchSetType(linesearch, SNESLINESEARCHL2);
351:       }
352:     }
353:   }
354:   if (qn->monflg) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &qn->monitor);
355:   return 0;
356: }

358: static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
359: {
360:   SNES_QN  *qn = (SNES_QN *)snes->data;
361:   PetscBool iascii;

363:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
364:   if (iascii) {
365:     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]);
366:     PetscViewerASCIIPrintf(viewer, "  Stored subspace size: %" PetscInt_FMT "\n", qn->m);
367:   }
368:   return 0;
369: }

371: /*@
372:     SNESQNSetRestartType - Sets the restart type for `SNESQN`.

374:     Logically Collective

376:     Input Parameters:
377: +   snes - the iterative context
378: -   rtype - restart type

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

384:     Level: intermediate

386:     `SNESQNRestartType`s:
387: +   `SNES_QN_RESTART_NONE` - never restart
388: .   `SNES_QN_RESTART_POWELL` - restart based upon descent criteria
389: -   `SNES_QN_RESTART_PERIODIC` - restart after a fixed number of iterations

391: .seealso: `SNESQN`, `SNESQNRestartType`, `SNES_QN_RESTART_NONE`, `SNES_QN_RESTART_POWELL`, `SNES_QN_RESTART_PERIODIC`
392: @*/
393: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
394: {
396:   PetscTryMethod(snes, "SNESQNSetRestartType_C", (SNES, SNESQNRestartType), (snes, rtype));
397:   return 0;
398: }

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

403:     Logically Collective

405:     Input Parameters:
406: +   snes - the nonlinear solver context
407: -   stype - scale type

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

412:     Level: intermediate

414:     `SNESQNScaleType`s:
415: +   `SNES_QN_SCALE_NONE` - don't scale the problem
416: .   `SNES_QN_SCALE_SCALAR` - use Shanno scaling
417: .   `SNES_QN_SCALE_DIAGONAL` - scale with a diagonalized BFGS formula (see Gilbert and Lemarechal 1989), available
418: -   `SNES_QN_SCALE_JACOBIAN` - scale by solving a linear system coming from the Jacobian you provided with SNESSetJacobian() computed at the first iteration
419:                              of QN and at ever restart.

421: .seealso: `SNES`, `SNESQN`, `SNESLineSearch`, `SNESQNScaleType`, `SNESSetJacobian()`
422: @*/

424: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
425: {
427:   PetscTryMethod(snes, "SNESQNSetScaleType_C", (SNES, SNESQNScaleType), (snes, stype));
428:   return 0;
429: }

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

435:   qn->scale_type = stype;
436:   if (stype == SNES_QN_SCALE_JACOBIAN) snes->usesksp = PETSC_TRUE;
437:   return 0;
438: }

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

444:   qn->restart_type = rtype;
445:   return 0;
446: }

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

451:     Logically Collective

453:     Input Parameters:
454: +   snes - the iterative context
455: -   qtype - variant type

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

460:     Level: beginner

462:     `SNESQNType`s:
463: +   `SNES_QN_LBFGS` - LBFGS variant
464: .   `SNES_QN_BROYDEN` - Broyden variant
465: -   `SNES_QN_BADBROYDEN` - Bad Broyden variant

467: .seealso: `SNESQN`, `SNES_QN_LBFGS`, `SNES_QN_BROYDEN`, `SNES_QN_BADBROYDEN`, `SNESQNType`, `TAOLMVM`, `TAOBLMVM`
468: @*/

470: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
471: {
473:   PetscTryMethod(snes, "SNESQNSetType_C", (SNES, SNESQNType), (snes, qtype));
474:   return 0;
475: }

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

481:   qn->type = qtype;
482:   return 0;
483: }

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

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

498:       References:
499: +   * -   Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995.
500: .   * -   R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi Newton Matrices and their use in Limited Memory Methods,
501:       Technical Report, Northwestern University, June 1992.
502: .   * -   Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker, Experiments with Quasi-Newton Methods in Solving Stiff ODE
503:       Systems, SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985.
504: .   * -   Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
505:        SIAM Review, 57(4), 2015
506: .   * -   Griewank, Andreas. "Broyden updating, the good and the bad!." Doc. Math (2012): 301-315.
507: .   * -   Gilbert, Jean Charles, and Claude Lemar{\'e}chal. "Some numerical experiments with variable-storage quasi-Newton algorithms."
508:       Mathematical programming 45.1-3 (1989): 407-435.
509: -   * -   Dener A., Munson T. "Accelerating Limited-Memory Quasi-Newton Convergence for Large-Scale Optimization"
510:       Computational Science - ICCS 2019. ICCS 2019. Lecture Notes in Computer Science, vol 11538. Springer, Cham

512:       Level: beginner

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

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

524:       Uses left nonlinear preconditioning by default.

526: .seealso: `SNESQNRestartType`, `SNESQNSetRestartType()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`,
527:           `SNESQNScaleType`, `SNESQNSetScaleType()`, `SNESQNSetType`, `SNESQNSetType ()`
528: M*/
529: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
530: {
531:   SNES_QN    *qn;
532:   const char *optionsprefix;

534:   snes->ops->setup          = SNESSetUp_QN;
535:   snes->ops->solve          = SNESSolve_QN;
536:   snes->ops->destroy        = SNESDestroy_QN;
537:   snes->ops->setfromoptions = SNESSetFromOptions_QN;
538:   snes->ops->view           = SNESView_QN;
539:   snes->ops->reset          = SNESReset_QN;

541:   snes->npcside = PC_LEFT;

543:   snes->usesnpc = PETSC_TRUE;
544:   snes->usesksp = PETSC_FALSE;

546:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

548:   if (!snes->tolerancesset) {
549:     snes->max_funcs = 30000;
550:     snes->max_its   = 10000;
551:   }

553:   PetscNew(&qn);
554:   snes->data       = (void *)qn;
555:   qn->m            = 10;
556:   qn->scaling      = 1.0;
557:   qn->monitor      = NULL;
558:   qn->monflg       = PETSC_FALSE;
559:   qn->powell_gamma = 0.9999;
560:   qn->scale_type   = SNES_QN_SCALE_DEFAULT;
561:   qn->restart_type = SNES_QN_RESTART_DEFAULT;
562:   qn->type         = SNES_QN_LBFGS;

564:   MatCreate(PetscObjectComm((PetscObject)snes), &qn->B);
565:   SNESGetOptionsPrefix(snes, &optionsprefix);
566:   MatSetOptionsPrefix(qn->B, optionsprefix);

568:   PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", SNESQNSetScaleType_QN);
569:   PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", SNESQNSetRestartType_QN);
570:   PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", SNESQNSetType_QN);
571:   return 0;
572: }