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:   (*snes->ops->converged)(snes,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:   if (snes->ops->update) {
110:     (*snes->ops->update)(snes, snes->iter);
111:   }

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

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

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

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

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

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

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

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

180:     /* general purpose update */
181:     if (snes->ops->update) {
182:       (*snes->ops->update)(snes, snes->iter);
183:     }

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

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


237:   if (!snes->vec_sol) {
238:     SNESGetDM(snes,&dm);
239:     DMCreateGlobalVector(dm,&snes->vec_sol);
240:   }
241:   SNESSetWorkVecs(snes,4);

243:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
244:     SNESSetUpMatrices(snes);
245:   }
246:   if (snes->npcside== PC_LEFT && snes->functype == SNES_FUNCTION_DEFAULT) {snes->functype = SNES_FUNCTION_UNPRECONDITIONED;}

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

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

303: static PetscErrorCode SNESReset_QN(SNES snes)
304: {
305:   SNES_QN        *qn;

307:   if (snes->data) {
308:     qn = (SNES_QN*)snes->data;
309:     MatDestroy(&qn->B);
310:   }
311:   return 0;
312: }

314: static PetscErrorCode SNESDestroy_QN(SNES snes)
315: {
316:   SNESReset_QN(snes);
317:   PetscFree(snes->data);
318:   PetscObjectComposeFunction((PetscObject)snes,"",NULL);
319:   return 0;
320: }

322: static PetscErrorCode SNESSetFromOptions_QN(PetscOptionItems *PetscOptionsObject,SNES snes)
323: {

325:   SNES_QN           *qn    = (SNES_QN*)snes->data;
326:   PetscBool         flg;
327:   SNESLineSearch    linesearch;
328:   SNESQNRestartType rtype = qn->restart_type;
329:   SNESQNScaleType   stype = qn->scale_type;
330:   SNESQNType        qtype = qn->type;

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

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

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

364: static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
365: {
366:   SNES_QN        *qn    = (SNES_QN*)snes->data;
367:   PetscBool      iascii;

369:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
370:   if (iascii) {
371:     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]);
372:     PetscViewerASCIIPrintf(viewer,"  Stored subspace size: %D\n", qn->m);
373:   }
374:   return 0;
375: }

377: /*@
378:     SNESQNSetRestartType - Sets the restart type for SNESQN.

380:     Logically Collective on SNES

382:     Input Parameters:
383: +   snes - the iterative context
384: -   rtype - restart type

386:     Options Database:
387: +   -snes_qn_restart_type <powell,periodic,none> - set the restart type
388: -   -snes_qn_m <m> - sets the number of stored updates and the restart period for periodic

390:     Level: intermediate

392:     SNESQNRestartTypes:
393: +   SNES_QN_RESTART_NONE - never restart
394: .   SNES_QN_RESTART_POWELL - restart based upon descent criteria
395: -   SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations

397: @*/
398: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
399: {
401:   PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));
402:   return 0;
403: }

405: /*@
406:     SNESQNSetScaleType - Sets the scaling type for the inner inverse Jacobian in SNESQN.

408:     Logically Collective on SNES

410:     Input Parameters:
411: +   snes - the iterative context
412: -   stype - scale type

414:     Options Database:
415: .   -snes_qn_scale_type <diagonal,none,scalar,jacobian> - Scaling type

417:     Level: intermediate

419:     SNESQNScaleTypes:
420: +   SNES_QN_SCALE_NONE - don't scale the problem
421: .   SNES_QN_SCALE_SCALAR - use Shanno scaling
422: .   SNES_QN_SCALE_DIAGONAL - scale with a diagonalized BFGS formula (see Gilbert and Lemarechal 1989), available
423: -   SNES_QN_SCALE_JACOBIAN - scale by solving a linear system coming from the Jacobian you provided with SNESSetJacobian() computed at the first iteration
424:                              of QN and at ever restart.

426: .seealso: SNES, SNESQN, SNESLineSearch, SNESQNScaleType, SNESSetJacobian()
427: @*/

429: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
430: {
432:   PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));
433:   return 0;
434: }

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

440:   qn->scale_type = stype;
441:   if (stype == SNES_QN_SCALE_JACOBIAN) snes->usesksp = PETSC_TRUE;
442:   return 0;
443: }

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

449:   qn->restart_type = rtype;
450:   return 0;
451: }

453: /*@
454:     SNESQNSetType - Sets the quasi-Newton variant to be used in SNESQN.

456:     Logically Collective on SNES

458:     Input Parameters:
459: +   snes - the iterative context
460: -   qtype - variant type

462:     Options Database:
463: .   -snes_qn_type <lbfgs,broyden,badbroyden> - quasi-Newton type

465:     Level: beginner

467:     SNESQNTypes:
468: +   SNES_QN_LBFGS - LBFGS variant
469: .   SNES_QN_BROYDEN - Broyden variant
470: -   SNES_QN_BADBROYDEN - Bad Broyden variant

472: @*/

474: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
475: {
477:   PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));
478:   return 0;
479: }

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

485:   qn->type = qtype;
486:   return 0;
487: }

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

493:       Options Database:

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

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

530:       Level: beginner

532: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR

534: M*/
535: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
536: {
537:   SNES_QN        *qn;
538:   const char     *optionsprefix;

540:   snes->ops->setup          = SNESSetUp_QN;
541:   snes->ops->solve          = SNESSolve_QN;
542:   snes->ops->destroy        = SNESDestroy_QN;
543:   snes->ops->setfromoptions = SNESSetFromOptions_QN;
544:   snes->ops->view           = SNESView_QN;
545:   snes->ops->reset          = SNESReset_QN;

547:   snes->npcside= PC_LEFT;

549:   snes->usesnpc = PETSC_TRUE;
550:   snes->usesksp = PETSC_FALSE;

552:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

554:   if (!snes->tolerancesset) {
555:     snes->max_funcs = 30000;
556:     snes->max_its   = 10000;
557:   }

559:   PetscNewLog(snes,&qn);
560:   snes->data          = (void*) qn;
561:   qn->m               = 10;
562:   qn->scaling         = 1.0;
563:   qn->monitor         = NULL;
564:   qn->monflg          = PETSC_FALSE;
565:   qn->powell_gamma    = 0.9999;
566:   qn->scale_type      = SNES_QN_SCALE_DEFAULT;
567:   qn->restart_type    = SNES_QN_RESTART_DEFAULT;
568:   qn->type            = SNES_QN_LBFGS;

570:   MatCreate(PetscObjectComm((PetscObject)snes), &qn->B);
571:   SNESGetOptionsPrefix(snes, &optionsprefix);
572:   MatSetOptionsPrefix(qn->B, optionsprefix);

574:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);
575:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);
576:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);
577:   return 0;
578: }