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 SNESQNGetMatrix_Private(SNES snes, Mat *B)
     24: {
     25:   SNES_QN    *qn = (SNES_QN *)snes->data;
     26:   const char *optionsprefix;

     28:   PetscFunctionBegin;
     29:   if (!qn->B) {
     30:     PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &qn->B));
     31:     PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
     32:     PetscCall(MatSetOptionsPrefix(qn->B, optionsprefix));
     33:     PetscCall(MatAppendOptionsPrefix(qn->B, "qn_"));
     34:     switch (qn->type) {
     35:     case SNES_QN_BROYDEN:
     36:       PetscCall(MatSetType(qn->B, MATLMVMBROYDEN));
     37:       qn->scale_type = SNES_QN_SCALE_NONE;
     38:       break;
     39:     case SNES_QN_BADBROYDEN:
     40:       PetscCall(MatSetType(qn->B, MATLMVMBADBROYDEN));
     41:       qn->scale_type = SNES_QN_SCALE_NONE;
     42:       break;
     43:     default:
     44:       PetscCall(MatSetType(qn->B, MATLMVMBFGS));
     45:       switch (qn->scale_type) {
     46:       case SNES_QN_SCALE_NONE:
     47:       case SNES_QN_SCALE_JACOBIAN:
     48:         PetscCall(MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_NONE));
     49:         break;
     50:       case SNES_QN_SCALE_SCALAR:
     51:       case SNES_QN_SCALE_DEFAULT:
     52:         PetscCall(MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_SCALAR));
     53:         break;
     54:       default:
     55:         break;
     56:       }
     57:       break;
     58:     }
     59:   }
     60:   *B = qn->B;
     61:   PetscFunctionReturn(PETSC_SUCCESS);
     62: }

     64: static PetscErrorCode SNESSolve_QN(SNES snes)
     65: {
     66:   SNES_QN             *qn = (SNES_QN *)snes->data;
     67:   Vec                  X, F, W;
     68:   Vec                  Y, D, Dold;
     69:   PetscInt             i, i_r;
     70:   PetscReal            fnorm, xnorm, ynorm;
     71:   SNESLineSearchReason lssucceed;
     72:   PetscBool            badstep, powell, periodic, restart;
     73:   PetscScalar          DolddotD, DolddotDold;
     74:   SNESConvergedReason  reason;
     75: #if defined(PETSC_USE_INFO)
     76:   PetscReal gnorm;
     77: #endif

     79:   PetscFunctionBegin;
     80:   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);

     82:   PetscCall(PetscCitationsRegister(SNESCitation, &SNEScite));
     83:   F = snes->vec_func;       /* residual vector */
     84:   Y = snes->vec_sol_update; /* search direction generated by J^-1D*/
     85:   X = snes->vec_sol;        /* solution vector */

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

     92:   snes->reason = SNES_CONVERGED_ITERATING;

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

     99:   if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
    100:     /* we need to sample the preconditioned function only once here,
    101:        since linesearch will always use the preconditioned function */
    102:     PetscCall(SNESApplyNPC(snes, X, NULL, F));
    103:     PetscCall(SNESGetConvergedReason(snes->npc, &reason));
    104:     if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
    105:       snes->reason = SNES_DIVERGED_INNER;
    106:       PetscFunctionReturn(PETSC_SUCCESS);
    107:     }
    108:   } else {
    109:     if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
    110:     else snes->vec_func_init_set = PETSC_FALSE;
    111:   }
    112:   PetscCall(VecNorm(F, NORM_2, &fnorm));
    113:   SNESCheckFunctionNorm(snes, fnorm);

    115:   PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    116:   snes->norm = fnorm;
    117:   PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
    118:   PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
    119:   PetscCall(SNESMonitor(snes, 0, fnorm));

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

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

    129:     /* general purpose update */
    130:     PetscTryTypeMethod(snes, update, snes->iter);

    132:     if (snes->npc) {
    133:       if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
    134:         PetscCall(SNESApplyNPC(snes, X, F, D));
    135:         PetscCall(SNESGetConvergedReason(snes->npc, &reason));
    136:         if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
    137:           snes->reason = SNES_DIVERGED_INNER;
    138:           PetscFunctionReturn(PETSC_SUCCESS);
    139:         }
    140:       } else if (snes->npcside == PC_RIGHT) {
    141:         PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, 0, 0));
    142:         PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
    143:         PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, 0, 0));
    144:         PetscCall(SNESGetConvergedReason(snes->npc, &reason));
    145:         if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
    146:           snes->reason = SNES_DIVERGED_INNER;
    147:           PetscFunctionReturn(PETSC_SUCCESS);
    148:         }
    149:         PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
    150:         PetscCall(VecCopy(F, D));
    151:       } else {
    152:         PetscCall(VecCopy(F, D));
    153:       }
    154:     } else {
    155:       PetscCall(VecCopy(F, D));
    156:     }

    158:     /* scale the initial update */
    159:     if (qn->scale_type == SNES_QN_SCALE_JACOBIAN && restart) {
    160:       PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
    161:       SNESCheckJacobianDomainerror(snes);
    162:       PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
    163:       PetscCall(MatLMVMSetJ0KSP(qn->B, snes->ksp));
    164:     }

    166:     /* update QN approx and calculate step */
    167:     PetscCall(MatLMVMUpdate(qn->B, X, D));
    168:     PetscCall(MatSolve(qn->B, D, Y));
    169:     PetscCall(VecDot(F, Y, &gdx));
    170:     if (PetscAbsScalar(gdx) <= 0) {
    171:       /* Step is not descent or solve was not successful
    172:          Use steepest descent direction */
    173:       PetscCall(VecCopy(F, Y));
    174:       PetscCall(MatLMVMReset(qn->B, PETSC_FALSE));
    175:     }

    177:     /* line search for lambda */
    178:     ynorm = 1;
    179: #if defined(PETSC_USE_INFO)
    180:     gnorm = fnorm;
    181: #endif
    182:     PetscCall(VecCopy(D, Dold));
    183:     PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y));
    184:     if (snes->reason == SNES_DIVERGED_FUNCTION_COUNT) break;
    185:     PetscCall(SNESLineSearchGetReason(snes->linesearch, &lssucceed));
    186:     PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm));
    187:     badstep = PETSC_FALSE;
    188:     if (lssucceed) {
    189:       if (++snes->numFailures >= snes->maxFailures) {
    190:         snes->reason = SNES_DIVERGED_LINE_SEARCH;
    191:         break;
    192:       }
    193:       badstep = PETSC_TRUE;
    194:     }

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

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

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

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

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

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

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

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

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

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

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

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

    312:   PetscFunctionBegin;
    313:   PetscCall(MatDestroy(&qn->B));
    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:   PetscOptionsHeadEnd();
    351:   PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));
    352:   PetscCall(MatSetFromOptions(qn->B));
    353:   if (!snes->linesearch) {
    354:     PetscCall(SNESGetLineSearch(snes, &linesearch));
    355:     if (!((PetscObject)linesearch)->type_name) {
    356:       if (qn->type == SNES_QN_LBFGS) {
    357:         PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHCP));
    358:       } else if (qn->type == SNES_QN_BROYDEN) {
    359:         PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC));
    360:       } else {
    361:         PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHSECANT));
    362:       }
    363:     }
    364:   }
    365:   if (qn->monflg) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &qn->monitor));
    366:   PetscFunctionReturn(PETSC_SUCCESS);
    367: }

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

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

    384: /*@
    385:   SNESQNSetRestartType - Sets the restart type for `SNESQN`.

    387:   Logically Collective

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

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

    397:   Level: intermediate

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

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

    413:   Logically Collective

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

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

    422:   Level: intermediate

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

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

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

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

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

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

    456:   Logically Collective

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

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

    465:   Level: intermediate

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

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

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

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

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

    499:       Level: beginner

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

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

    511:       Uses left nonlinear preconditioning by default.

    513:       See {cite}`byrd1994representations`, {cite}`kelley95`, {cite}`brown1985experiments`, {cite}`griewank2012broyden`, {cite}`gilbert1989some`,
    514:       {cite}`dener2019accelerating`, and {cite}`bruneknepleysmithtu15`

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

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

    531:   snes->npcside = PC_LEFT;

    533:   snes->usesnpc = PETSC_TRUE;
    534:   snes->usesksp = PETSC_FALSE;

    536:   snes->alwayscomputesfinalresidual = PETSC_TRUE;

    538:   PetscCall(SNESParametersInitialize(snes));
    539:   PetscObjectParameterSetDefault(snes, max_funcs, 30000);
    540:   PetscObjectParameterSetDefault(snes, max_its, 10000);

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

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