Actual source code: qn.c

petsc-3.13.1 2020-05-02
Report Typos and Errors
  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_",0};
  7: const char *const SNESQNRestartTypes[] =      {"DEFAULT","NONE","POWELL","PERIODIC","SNESQNRestartType","SNES_QN_RESTART_",0};
  8: const char *const SNESQNTypes[] =             {"LBFGS","BROYDEN","BADBROYDEN","SNESQNType","SNES_QN_",0};

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

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

 40:   if (snes->xl || snes->xu || snes->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE, "SNES solver %s does not support bounds", ((PetscObject)snes)->type_name);

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

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

 53:   snes->reason = SNES_CONVERGED_ITERATING;

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

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

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

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

 93:   /* test convergence */
 94:   (*snes->ops->converged)(snes,0,0.0,0.0,fnorm,&snes->reason,snes->cnvP);
 95:   if (snes->reason) return(0);

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

110:   /* general purpose update */
111:   if (snes->ops->update) {
112:     (*snes->ops->update)(snes, snes->iter);
113:   }

115:   /* scale the initial update */
116:   if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
117:     SNESComputeJacobian(snes,X,snes->jacobian,snes->jacobian_pre);
118:     SNESCheckJacobianDomainerror(snes);
119:     KSPSetOperators(snes->ksp,snes->jacobian,snes->jacobian_pre);
120:     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:     MatLMVMUpdate(qn->B, X, D);
126:     MatSolve(qn->B, D, Y);

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

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

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

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

165:     SNESLogConvergenceHistory(snes,snes->norm,snes->iter);
166:     SNESMonitor(snes,snes->iter,snes->norm);

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

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

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

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


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

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

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

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

307: static PetscErrorCode SNESReset_QN(SNES snes)
308: {
310:   SNES_QN        *qn;

313:   if (snes->data) {
314:     qn = (SNES_QN*)snes->data;
315:     MatDestroy(&qn->B);
316:   }
317:   return(0);
318: }

320: static PetscErrorCode SNESDestroy_QN(SNES snes)
321: {

325:   SNESReset_QN(snes);
326:   PetscFree(snes->data);
327:   PetscObjectComposeFunction((PetscObject)snes,"",NULL);
328:   return(0);
329: }

331: static PetscErrorCode SNESSetFromOptions_QN(PetscOptionItems *PetscOptionsObject,SNES snes)
332: {

334:   PetscErrorCode    ierr;
335:   SNES_QN           *qn    = (SNES_QN*)snes->data;
336:   PetscBool         flg;
337:   SNESLineSearch    linesearch;
338:   SNESQNRestartType rtype = qn->restart_type;
339:   SNESQNScaleType   stype = qn->scale_type;
340:   SNESQNType        qtype = qn->type;

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

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

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

375: static PetscErrorCode SNESView_QN(SNES snes, PetscViewer viewer)
376: {
377:   SNES_QN        *qn    = (SNES_QN*)snes->data;
378:   PetscBool      iascii;

382:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
383:   if (iascii) {
384:     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]);
385:     PetscViewerASCIIPrintf(viewer,"  Stored subspace size: %D\n", qn->m);
386:   }
387:   return(0);
388: }

390: /*@
391:     SNESQNSetRestartType - Sets the restart type for SNESQN.

393:     Logically Collective on SNES

395:     Input Parameters:
396: +   snes - the iterative context
397: -   rtype - restart type

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

403:     Level: intermediate

405:     SNESQNRestartTypes:
406: +   SNES_QN_RESTART_NONE - never restart
407: .   SNES_QN_RESTART_POWELL - restart based upon descent criteria
408: -   SNES_QN_RESTART_PERIODIC - restart after a fixed number of iterations

410: @*/
411: PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
412: {

417:   PetscTryMethod(snes,"SNESQNSetRestartType_C",(SNES,SNESQNRestartType),(snes,rtype));
418:   return(0);
419: }

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

424:     Logically Collective on SNES

426:     Input Parameters:
427: +   snes - the iterative context
428: -   stype - scale type

430:     Options Database:
431: .   -snes_qn_scale_type <diagonal,none,scalar,jacobian>

433:     Level: intermediate

435:     SNESQNScaleTypes:
436: +   SNES_QN_SCALE_NONE - don't scale the problem
437: .   SNES_QN_SCALE_SCALAR - use shanno scaling
438: .   SNES_QN_SCALE_DIAGONAL - scale with a diagonalized BFGS formula (see Gilbert and Lemarechal 1989), available 
439: -   SNES_QN_SCALE_JACOBIAN - scale by solving a linear system coming from the Jacobian you provided with SNESSetJacobian() computed at the first iteration
440:                              of QN and at ever restart.

442: .seealso: SNES, SNESQN, SNESLineSearch, SNESQNScaleType, SNESetJacobian()
443: @*/

445: PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
446: {

451:   PetscTryMethod(snes,"SNESQNSetScaleType_C",(SNES,SNESQNScaleType),(snes,stype));
452:   return(0);
453: }

455: PetscErrorCode SNESQNSetScaleType_QN(SNES snes, SNESQNScaleType stype)
456: {
457:   SNES_QN *qn = (SNES_QN*)snes->data;

460:   qn->scale_type = stype;
461:   return(0);
462: }

464: PetscErrorCode SNESQNSetRestartType_QN(SNES snes, SNESQNRestartType rtype)
465: {
466:   SNES_QN *qn = (SNES_QN*)snes->data;

469:   qn->restart_type = rtype;
470:   return(0);
471: }

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

476:     Logically Collective on SNES

478:     Input Parameters:
479: +   snes - the iterative context
480: -   qtype - variant type

482:     Options Database:
483: .   -snes_qn_type <lbfgs,broyden,badbroyden>

485:     Level: beginner

487:     SNESQNTypes:
488: +   SNES_QN_LBFGS - LBFGS variant
489: .   SNES_QN_BROYDEN - Broyden variant
490: -   SNES_QN_BADBROYDEN - Bad Broyden variant

492: @*/

494: PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
495: {

500:   PetscTryMethod(snes,"SNESQNSetType_C",(SNES,SNESQNType),(snes,qtype));
501:   return(0);
502: }

504: PetscErrorCode SNESQNSetType_QN(SNES snes, SNESQNType qtype)
505: {
506:   SNES_QN *qn = (SNES_QN*)snes->data;

509:   qn->type = qtype;
510:   return(0);
511: }

513: /* -------------------------------------------------------------------------- */
514: /*MC
515:       SNESQN - Limited-Memory Quasi-Newton methods for the solution of nonlinear systems.

517:       Options Database:

519: +     -snes_qn_m <m> - Number of past states saved for the L-Broyden methods.
520: +     -snes_qn_restart_type <powell,periodic,none> - set the restart type
521: .     -snes_qn_powell_gamma - Angle condition for restart.
522: .     -snes_qn_powell_descent - Descent condition for restart.
523: .     -snes_qn_type <lbfgs,broyden,badbroyden> - QN type
524: .     -snes_qn_scale_type <diagonal,none,scalar,jacobian> - scaling performed on inner Jacobian
525: .     -snes_linesearch_type <cp, l2, basic> - Type of line search.
526: -     -snes_qn_monitor - Monitors the quasi-newton Jacobian.

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

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

538:       Uses left nonlinear preconditioning by default.

540:       References:
541: +   1. -   Kelley, C.T., Iterative Methods for Linear and Nonlinear Equations, Chapter 8, SIAM, 1995.
542: .   2. -   R. Byrd, J. Nocedal, R. Schnabel, Representations of Quasi Newton Matrices and their use in Limited Memory Methods,
543:       Technical Report, Northwestern University, June 1992.
544: .   3. -   Peter N. Brown, Alan C. Hindmarsh, Homer F. Walker, Experiments with Quasi-Newton Methods in Solving Stiff ODE
545:       Systems, SIAM J. Sci. Stat. Comput. Vol 6(2), April 1985.
546: .   4. -   Peter R. Brune, Matthew G. Knepley, Barry F. Smith, and Xuemin Tu, "Composing Scalable Nonlinear Algebraic Solvers",
547:        SIAM Review, 57(4), 2015
548: .   5. -   Griewank, Andreas. "Broyden updating, the good and the bad!." Doc. Math (2012): 301-315.
549: .   6. -   Gilbert, Jean Charles, and Claude Lemaréchal. "Some numerical experiments with variable-storage quasi-Newton algorithms." 
550:       Mathematical programming 45.1-3 (1989): 407-435.
551: -   7. -   Dener A., Munson T. "Accelerating Limited-Memory Quasi-Newton Convergence for Large-Scale Optimization" 
552:       Computational Science – ICCS 2019. ICCS 2019. Lecture Notes in Computer Science, vol 11538. Springer, Cham

554:       Level: beginner

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

558: M*/
559: PETSC_EXTERN PetscErrorCode SNESCreate_QN(SNES snes)
560: {
562:   SNES_QN        *qn;
563:   const char     *optionsprefix;

566:   snes->ops->setup          = SNESSetUp_QN;
567:   snes->ops->solve          = SNESSolve_QN;
568:   snes->ops->destroy        = SNESDestroy_QN;
569:   snes->ops->setfromoptions = SNESSetFromOptions_QN;
570:   snes->ops->view           = SNESView_QN;
571:   snes->ops->reset          = SNESReset_QN;

573:   snes->npcside= PC_LEFT;

575:   snes->usesnpc = PETSC_TRUE;
576:   snes->usesksp = PETSC_FALSE;

578:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

580:   if (!snes->tolerancesset) {
581:     snes->max_funcs = 30000;
582:     snes->max_its   = 10000;
583:   }

585:   PetscNewLog(snes,&qn);
586:   snes->data          = (void*) qn;
587:   qn->m               = 10;
588:   qn->scaling         = 1.0;
589:   qn->monitor         = NULL;
590:   qn->monflg          = PETSC_FALSE;
591:   qn->powell_gamma    = 0.9999;
592:   qn->scale_type      = SNES_QN_SCALE_DEFAULT;
593:   qn->restart_type    = SNES_QN_RESTART_DEFAULT;
594:   qn->type            = SNES_QN_LBFGS;

596:   MatCreate(PetscObjectComm((PetscObject)snes), &qn->B);
597:   SNESGetOptionsPrefix(snes, &optionsprefix);
598:   MatSetOptionsPrefix(qn->B, optionsprefix);

600:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetScaleType_C",SNESQNSetScaleType_QN);
601:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetRestartType_C",SNESQNSetRestartType_QN);
602:   PetscObjectComposeFunction((PetscObject)snes,"SNESQNSetType_C",SNESQNSetType_QN);
603:   return(0);
604: }