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;
 31:   SNESLineSearchReason lssucceed;
 32:   PetscBool            badstep, powell, periodic;
 33:   PetscScalar          DolddotD, DolddotDold;
 34:   SNESConvergedReason  reason;
 35: #if defined(PETSC_USE_INFO)
 36:   PetscReal gnorm;
 37: #endif

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

 41:   PetscFunctionBegin;
 42:   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);

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

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

 55:   snes->reason = SNES_CONVERGED_ITERATING;

 57:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
 58:   snes->iter = 0;
 59:   snes->norm = 0.;
 60:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

 62:   if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
 63:     PetscCall(SNESApplyNPC(snes, X, NULL, F));
 64:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
 65:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 66:       snes->reason = SNES_DIVERGED_INNER;
 67:       PetscFunctionReturn(PETSC_SUCCESS);
 68:     }
 69:     PetscCall(VecNorm(F, NORM_2, &fnorm));
 70:   } else {
 71:     if (!snes->vec_func_init_set) {
 72:       PetscCall(SNESComputeFunction(snes, X, F));
 73:     } else snes->vec_func_init_set = PETSC_FALSE;

 75:     PetscCall(VecNorm(F, NORM_2, &fnorm));
 76:     SNESCheckFunctionNorm(snes, fnorm);
 77:   }
 78:   if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
 79:     PetscCall(SNESApplyNPC(snes, X, F, D));
 80:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
 81:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
 82:       snes->reason = SNES_DIVERGED_INNER;
 83:       PetscFunctionReturn(PETSC_SUCCESS);
 84:     }
 85:   } else {
 86:     PetscCall(VecCopy(F, D));
 87:   }

 89:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
 90:   snes->norm = fnorm;
 91:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
 92:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
 93:   PetscCall(SNESMonitor(snes, 0, fnorm));

 95:   /* test convergence */
 96:   PetscUseTypeMethod(snes, converged, 0, 0.0, 0.0, fnorm, &snes->reason, snes->cnvP);
 97:   if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);

 99:   if (snes->npc && snes->npcside == PC_RIGHT) {
100:     PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, 0, 0));
101:     PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
102:     PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, 0, 0));
103:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
104:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
105:       snes->reason = SNES_DIVERGED_INNER;
106:       PetscFunctionReturn(PETSC_SUCCESS);
107:     }
108:     PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
109:     PetscCall(VecCopy(F, D));
110:   }

112:   /* general purpose update */
113:   PetscTryTypeMethod(snes, update, snes->iter);

115:   /* scale the initial update */
116:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
117:     PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
118:     SNESCheckJacobianDomainerror(snes);
119:     PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
120:     PetscCall(MatLMVMSetJ0KSP(qn->B, snes->ksp));
121:   }

123:   for (i = 0, i_r = 0; i < snes->max_its; i++, i_r++) {
124:     /* update QN approx and calculate step */
125:     PetscCall(MatLMVMUpdate(qn->B, X, D));
126:     PetscCall(MatSolve(qn->B, D, Y));

128:     /* line search for lambda */
129:     ynorm = 1;
130: #if defined(PETSC_USE_INFO)
131:     gnorm = fnorm;
132: #endif
133:     PetscCall(VecCopy(D, Dold));
134:     PetscCall(VecCopy(X, Xold));
135:     PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y));
136:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
137:     PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed));
138:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm));
139:     badstep = PETSC_FALSE;
140:     if (lssucceed) {
141:       if (++snes->numFailures >= snes->maxFailures) {
142:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
143:         break;
144:       }
145:       badstep = PETSC_TRUE;
146:     }

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

151:     if (snes->npc && snes->npcside == PC_RIGHT) {
152:       PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, 0, 0));
153:       PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
154:       PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, 0, 0));
155:       PetscCall(SNESGetConvergedReason(snes->npc, &reason));
156:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
157:         snes->reason = SNES_DIVERGED_INNER;
158:         PetscFunctionReturn(PETSC_SUCCESS);
159:       }
160:       PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
161:     }

163:     PetscCall(SNESSetIterationNumber(snes, i + 1));
164:     snes->norm  = fnorm;
165:     snes->xnorm = xnorm;
166:     snes->ynorm = ynorm;

168:     PetscCall(SNESLogConvergenceHistory(snes, snes->norm, snes->iter));
169:     PetscCall(SNESMonitor(snes, snes->iter, snes->norm));

171:     /* set parameter for default relative tolerance convergence test */
172:     PetscUseTypeMethod(snes, converged, snes->iter, xnorm, ynorm, fnorm, &snes->reason, snes->cnvP);
173:     if (snes->reason) PetscFunctionReturn(PETSC_SUCCESS);
174:     if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
175:       PetscCall(SNESApplyNPC(snes, X, F, D));
176:       PetscCall(SNESGetConvergedReason(snes->npc, &reason));
177:       if (reason < 0 && reason != SNES_DIVERGED_MAX_IT) {
178:         snes->reason = SNES_DIVERGED_INNER;
179:         PetscFunctionReturn(PETSC_SUCCESS);
180:       }
181:     } else {
182:       PetscCall(VecCopy(F, D));
183:     }

185:     /* general purpose update */
186:     PetscTryTypeMethod(snes, update, snes->iter);

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

233: static PetscErrorCode SNESSetUp_QN(SNES snes)
234: {
235:   SNES_QN *qn = (SNES_QN *)snes->data;
236:   DM       dm;
237:   PetscInt n, N;

239:   PetscFunctionBegin;

241:   if (!snes->vec_sol) {
242:     PetscCall(SNESGetDM(snes, &dm));
243:     PetscCall(DMCreateGlobalVector(dm, &snes->vec_sol));
244:   }
245:   PetscCall(SNESSetWorkVecs(snes, 4));

247:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) PetscCall(SNESSetUpMatrices(snes));
248:   if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) snes->functype = SNES_FUNCTION_UNPRECONDITIONED;

250:   /* set method defaults */
251:   if (qn->scale_type == SNES_QN_SCALE_DEFAULT) {
252:     if (qn->type == SNES_QN_BADBROYDEN) {
253:       qn->scale_type = SNES_QN_SCALE_NONE;
254:     } else {
255:       qn->scale_type = SNES_QN_SCALE_SCALAR;
256:     }
257:   }
258:   if (qn->restart_type == SNES_QN_RESTART_DEFAULT) {
259:     if (qn->type == SNES_QN_LBFGS) {
260:       qn->restart_type = SNES_QN_RESTART_POWELL;
261:     } else {
262:       qn->restart_type = SNES_QN_RESTART_PERIODIC;
263:     }
264:   }

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

305: static PetscErrorCode SNESReset_QN(SNES snes)
306: {
307:   SNES_QN *qn;

309:   PetscFunctionBegin;
310:   if (snes->data) {
311:     qn = (SNES_QN *)snes->data;
312:     PetscCall(MatDestroy(&qn->B));
313:   }
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }

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

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

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

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

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

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

373:   PetscFunctionBegin;
374:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
375:   if (iascii) {
376:     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]));
377:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Stored subspace size: %" PetscInt_FMT "\n", qn->m));
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

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:     `SNESQNRestartType`s:
398: +   `SNES_QN_RESTART_NONE` - never restart
399: .   `SNES_QN_RESTART_POWELL` - restart based upon descent criteria
400: -   `SNES_QN_RESTART_PERIODIC` - restart after a fixed number of iterations

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

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

415:     Logically Collective

417:     Input Parameters:
418: +   snes - the nonlinear solver context
419: -   stype - scale type

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

424:     Level: intermediate

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

433: .seealso: `SNES`, `SNESQN`, `SNESLineSearch`, `SNESQNScaleType`, `SNESSetJacobian()`
434: @*/

436: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
437: {
438:   PetscFunctionBegin;
440:   PetscTryMethod(snes, "SNESQNSetScaleType_C", (SNES, SNESQNScaleType), (snes, stype));
441:   PetscFunctionReturn(PETSC_SUCCESS);
442: }

444: PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
445: {
446:   SNES_QN *qn = (SNES_QN *)snes->data;

448:   PetscFunctionBegin;
449:   qn->scale_type = stype;
450:   if (stype == SNES_QN_SCALE_JACOBIAN) snes->usesksp = PETSC_TRUE;
451:   PetscFunctionReturn(PETSC_SUCCESS);
452: }

454: PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
455: {
456:   SNES_QN *qn = (SNES_QN *)snes->data;

458:   PetscFunctionBegin;
459:   qn->restart_type = rtype;
460:   PetscFunctionReturn(PETSC_SUCCESS);
461: }

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

466:     Logically Collective

468:     Input Parameters:
469: +   snes - the iterative context
470: -   qtype - variant type

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

475:     Level: beginner

477:     `SNESQNType`s:
478: +   `SNES_QN_LBFGS` - LBFGS variant
479: .   `SNES_QN_BROYDEN` - Broyden variant
480: -   `SNES_QN_BADBROYDEN` - Bad Broyden variant

482: .seealso: `SNESQN`, `SNES_QN_LBFGS`, `SNES_QN_BROYDEN`, `SNES_QN_BADBROYDEN`, `SNESQNType`, `TAOLMVM`, `TAOBLMVM`
483: @*/

485: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
486: {
487:   PetscFunctionBegin;
489:   PetscTryMethod(snes, "SNESQNSetType_C", (SNES, SNESQNType), (snes, qtype));
490:   PetscFunctionReturn(PETSC_SUCCESS);
491: }

493: PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
494: {
495:   SNES_QN *qn = (SNES_QN *)snes->data;

497:   PetscFunctionBegin;
498:   qn->type = qtype;
499:   PetscFunctionReturn(PETSC_SUCCESS);
500: }

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

505:       Options Database Keys:
506: +     -snes_qn_m <m> - Number of past states saved for the L-Broyden methods.
507: .     -snes_qn_restart_type <powell,periodic,none> - set the restart type
508: .     -snes_qn_powell_gamma - Angle condition for restart.
509: .     -snes_qn_powell_descent - Descent condition for restart.
510: .     -snes_qn_type <lbfgs,broyden,badbroyden> - QN type
511: .     -snes_qn_scale_type <diagonal,none,scalar,jacobian> - scaling performed on inner Jacobian
512: .     -snes_linesearch_type <cp, l2, basic> - Type of line search.
513: -     -snes_qn_monitor - Monitors the quasi-newton Jacobian.

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

529:       Level: beginner

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

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

541:       Uses left nonlinear preconditioning by default.

543: .seealso: `SNESQNRestartType`, `SNESQNSetRestartType()`, `SNESCreate()`, `SNES`, `SNESSetType()`, `SNESNEWTONLS`, `SNESNEWTONTR`,
544:           `SNESQNScaleType`, `SNESQNSetScaleType()`, `SNESQNSetType`, `SNESQNSetType()`
545: M*/
546: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
547: {
548:   SNES_QN    *qn;
549:   const char *optionsprefix;

551:   PetscFunctionBegin;
552:   snes->ops->setup          = SNESSetUp_QN;
553:   snes->ops->solve          = SNESSolve_QN;
554:   snes->ops->destroy        = SNESDestroy_QN;
555:   snes->ops->setfromoptions = SNESSetFromOptions_QN;
556:   snes->ops->view           = SNESView_QN;
557:   snes->ops->reset          = SNESReset_QN;

559:   snes->npcside = PC_LEFT;

561:   snes->usesnpc = PETSC_TRUE;
562:   snes->usesksp = PETSC_FALSE;

564:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

566:   if (!snes->tolerancesset) {
567:     snes->max_funcs = 30000;
568:     snes->max_its   = 10000;
569:   }

571:   PetscCall(PetscNew(&qn));
572:   snes->data       = (void *)qn;
573:   qn->m            = 10;
574:   qn->scaling      = 1.0;
575:   qn->monitor      = NULL;
576:   qn->monflg       = PETSC_FALSE;
577:   qn->powell_gamma = 0.9999;
578:   qn->scale_type   = SNES_QN_SCALE_DEFAULT;
579:   qn->restart_type = SNES_QN_RESTART_DEFAULT;
580:   qn->type         = SNES_QN_LBFGS;

582:   PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &qn->B));
583:   PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
584:   PetscCall(MatSetOptionsPrefix(qn->B, optionsprefix));

586:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", SNESQNSetScaleType_QN));
587:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", SNESQNSetRestartType_QN));
588:   PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", SNESQNSetType_QN));
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }