#include <petsc/private/snesimpl.h> /*I "petscsnes.h" I*/
#include <petscdm.h>

const char *const SNESQNScaleTypes[]   = {"DEFAULT", "NONE", "SCALAR", "DIAGONAL", "JACOBIAN", "SNESQNScaleType", "SNES_QN_SCALING_", NULL};
const char *const SNESQNRestartTypes[] = {"DEFAULT", "NONE", "POWELL", "PERIODIC", "SNESQNRestartType", "SNES_QN_RESTART_", NULL};
const char *const SNESQNTypes[]        = {"LBFGS", "BROYDEN", "BADBROYDEN", "SNESQNType", "SNES_QN_", NULL};

typedef struct {
  Mat               B;      /* Quasi-Newton approximation Matrix (MATLMVM) */
  PetscInt          m;      /* The number of kept previous steps */
  PetscReal        *lambda; /* The line search history of the method */
  PetscBool         monflg;
  PetscViewer       monitor;
  PetscReal         powell_gamma; /* Powell angle restart condition */
  PetscReal         scaling;      /* scaling of H0 */
  SNESQNType        type;         /* the type of quasi-newton method used */
  SNESQNScaleType   scale_type;   /* the type of scaling used */
  SNESQNRestartType restart_type; /* determine the frequency and type of restart conditions */
} SNES_QN;

static PetscErrorCode SNESQNGetMatrix_Private(SNES snes, Mat *B)
{
  SNES_QN    *qn = (SNES_QN *)snes->data;
  const char *optionsprefix;

  PetscFunctionBegin;
  if (!qn->B) {
    PetscCall(MatCreate(PetscObjectComm((PetscObject)snes), &qn->B));
    PetscCall(SNESGetOptionsPrefix(snes, &optionsprefix));
    PetscCall(MatSetOptionsPrefix(qn->B, optionsprefix));
    PetscCall(MatAppendOptionsPrefix(qn->B, "qn_"));
    switch (qn->type) {
    case SNES_QN_BROYDEN:
      PetscCall(MatSetType(qn->B, MATLMVMBROYDEN));
      qn->scale_type = SNES_QN_SCALE_NONE;
      break;
    case SNES_QN_BADBROYDEN:
      PetscCall(MatSetType(qn->B, MATLMVMBADBROYDEN));
      qn->scale_type = SNES_QN_SCALE_NONE;
      break;
    default:
      PetscCall(MatSetType(qn->B, MATLMVMBFGS));
      switch (qn->scale_type) {
      case SNES_QN_SCALE_NONE:
      case SNES_QN_SCALE_JACOBIAN:
        PetscCall(MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_NONE));
        break;
      case SNES_QN_SCALE_SCALAR:
      case SNES_QN_SCALE_DEFAULT:
        PetscCall(MatLMVMSymBroydenSetScaleType(qn->B, MAT_LMVM_SYMBROYDEN_SCALE_SCALAR));
        break;
      default:
        break;
      }
      break;
    }
  }
  *B = qn->B;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode SNESSolve_QN(SNES snes)
{
  SNES_QN             *qn = (SNES_QN *)snes->data;
  Vec                  X, F, W;
  Vec                  Y, D, Dold;
  PetscInt             i, i_r;
  PetscReal            fnorm, xnorm, ynorm;
  SNESLineSearchReason lsreason;
  PetscBool            powell, periodic, restart;
  PetscScalar          DolddotD, DolddotDold;
  SNESConvergedReason  reason;
#if defined(PETSC_USE_INFO)
  PetscReal gnorm;
#endif

  PetscFunctionBegin;
  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);

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

  /* directions */
  W    = snes->work[0];
  D    = snes->work[1];
  Dold = snes->work[2];

  snes->reason = SNES_CONVERGED_ITERATING;

  PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
  snes->iter = 0;
  snes->norm = 0.;
  PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

  if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_PRECONDITIONED) {
    /* we need to sample the preconditioned function only once here,
       since linesearch will always use the preconditioned function */
    PetscCall(SNESApplyNPC(snes, X, NULL, F));
    PetscCall(SNESGetConvergedReason(snes->npc, &reason));
    if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
      snes->reason = SNES_DIVERGED_INNER;
      PetscFunctionReturn(PETSC_SUCCESS);
    }
  } else {
    if (!snes->vec_func_init_set) PetscCall(SNESComputeFunction(snes, X, F));
    else snes->vec_func_init_set = PETSC_FALSE;
  }
  PetscCall(VecNorm(F, NORM_2, &fnorm));
  SNESCheckFunctionDomainError(snes, fnorm);

  PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
  snes->norm = fnorm;
  PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));
  PetscCall(SNESLogConvergenceHistory(snes, fnorm, 0));
  PetscCall(SNESMonitor(snes, 0, fnorm));

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

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

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

    if (snes->npc) {
      if (snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
        PetscCall(SNESApplyNPC(snes, X, F, D));
        PetscCall(SNESGetConvergedReason(snes->npc, &reason));
        if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
          snes->reason = SNES_DIVERGED_INNER;
          PetscFunctionReturn(PETSC_SUCCESS);
        }
      } else if (snes->npcside == PC_RIGHT) {
        PetscCall(PetscLogEventBegin(SNES_NPCSolve, snes->npc, X, 0, 0));
        PetscCall(SNESSolve(snes->npc, snes->vec_rhs, X));
        PetscCall(PetscLogEventEnd(SNES_NPCSolve, snes->npc, X, 0, 0));
        PetscCall(SNESGetConvergedReason(snes->npc, &reason));
        if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
          snes->reason = SNES_DIVERGED_INNER;
          PetscFunctionReturn(PETSC_SUCCESS);
        }
        PetscCall(SNESGetNPCFunction(snes, F, &fnorm));
        PetscCall(VecCopy(F, D));
      } else {
        PetscCall(VecCopy(F, D));
      }
    } else {
      PetscCall(VecCopy(F, D));
    }

    /* scale the initial update */
    if (qn->scale_type == SNES_QN_SCALE_JACOBIAN && restart) {
      PetscCall(SNESComputeJacobian(snes, X, snes->jacobian, snes->jacobian_pre));
      SNESCheckJacobianDomainError(snes);
      PetscCall(KSPSetOperators(snes->ksp, snes->jacobian, snes->jacobian_pre));
      PetscCall(MatLMVMSetJ0KSP(qn->B, snes->ksp));
    }

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

    /* line search for lambda */
    ynorm = 1;
#if defined(PETSC_USE_INFO)
    gnorm = fnorm;
#endif
    PetscCall(VecCopy(D, Dold));
    PetscCall(SNESLineSearchApply(snes->linesearch, X, F, &fnorm, Y));
    if (snes->reason) break;
    SNESCheckLineSearchFailure(snes);
    PetscCall(SNESLineSearchGetReason(snes->linesearch, &lsreason));

    PetscCall(SNESLineSearchGetNorms(snes->linesearch, &xnorm, &fnorm, &ynorm));
    /* convergence monitoring */
    PetscCall(PetscInfo(snes, "fnorm=%18.16e, gnorm=%18.16e, ynorm=%18.16e, lssucceed=%d\n", (double)fnorm, (double)gnorm, (double)ynorm, (int)lsreason));

    PetscCall(PetscObjectSAWsTakeAccess((PetscObject)snes));
    snes->iter  = i + 1;
    snes->norm  = fnorm;
    snes->xnorm = xnorm;
    snes->ynorm = ynorm;
    PetscCall(PetscObjectSAWsGrantAccess((PetscObject)snes));

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

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

    if (snes->npc && snes->npcside == PC_LEFT && snes->functype == SNES_FUNCTION_UNPRECONDITIONED) {
      PetscCall(SNESApplyNPC(snes, X, F, D));
      PetscCall(SNESGetConvergedReason(snes->npc, &reason));
      if (reason < 0 && reason != SNES_DIVERGED_MAX_IT && reason != SNES_DIVERGED_TR_DELTA) {
        snes->reason = SNES_DIVERGED_INNER;
        PetscFunctionReturn(PETSC_SUCCESS);
      }
    } else {
      PetscCall(VecCopy(F, D));
    }

    /* restart conditions */
    powell = PETSC_FALSE;
    if (qn->restart_type == SNES_QN_RESTART_POWELL && i_r > 1) {
      /* check restart by Powell's Criterion: |F^T H_0 Fold| > powell_gamma * |Fold^T H_0 Fold| */
      if (qn->scale_type == SNES_QN_SCALE_JACOBIAN) {
        PetscCall(MatMult(snes->jacobian, Dold, W));
      } else {
        PetscCall(VecCopy(Dold, W));
      }
      PetscCall(VecDotBegin(W, Dold, &DolddotDold));
      PetscCall(VecDotBegin(W, D, &DolddotD));
      PetscCall(VecDotEnd(W, Dold, &DolddotDold));
      PetscCall(VecDotEnd(W, D, &DolddotD));
      if (PetscAbs(PetscRealPart(DolddotD)) > qn->powell_gamma * PetscAbs(PetscRealPart(DolddotDold))) powell = PETSC_TRUE;
    }
    periodic = PETSC_FALSE;
    if (qn->restart_type == SNES_QN_RESTART_PERIODIC) {
      if (i_r > qn->m - 1) periodic = PETSC_TRUE;
    }

    /* restart if either powell or periodic restart is satisfied. */
    restart = PETSC_FALSE;
    if (lsreason || powell || periodic) {
      restart = PETSC_TRUE;
      if (qn->monflg) {
        PetscCall(PetscViewerASCIIAddTab(qn->monitor, ((PetscObject)snes)->tablevel + 2));
        if (powell) {
          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));
        } else {
          PetscCall(PetscViewerASCIIPrintf(qn->monitor, "Periodic restart! i_r = %" PetscInt_FMT "\n", i_r));
        }
        PetscCall(PetscViewerASCIISubtractTab(qn->monitor, ((PetscObject)snes)->tablevel + 2));
      }
      i_r = -1;
      PetscCall(MatLMVMReset(qn->B, PETSC_FALSE));
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode SNESSetUp_QN(SNES snes)
{
  SNES_QN *qn = (SNES_QN *)snes->data;
  DM       dm;
  PetscInt n, N;

  PetscFunctionBegin;
  if (!snes->vec_sol) {
    PetscCall(SNESGetDM(snes, &dm));
    PetscCall(DMCreateGlobalVector(dm, &snes->vec_sol));
  }
  PetscCall(SNESSetWorkVecs(snes, 3));
  PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));

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

  /* set method defaults */
  if (qn->scale_type == SNES_QN_SCALE_DEFAULT) {
    if (qn->type == SNES_QN_BADBROYDEN) {
      qn->scale_type = SNES_QN_SCALE_NONE;
    } else {
      qn->scale_type = SNES_QN_SCALE_SCALAR;
    }
  }
  if (qn->restart_type == SNES_QN_RESTART_DEFAULT) {
    if (qn->type == SNES_QN_LBFGS) {
      qn->restart_type = SNES_QN_RESTART_POWELL;
    } else {
      qn->restart_type = SNES_QN_RESTART_PERIODIC;
    }
  }
  /* Set up the LMVM matrix */
  PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));
  PetscCall(VecGetLocalSize(snes->vec_sol, &n));
  PetscCall(VecGetSize(snes->vec_sol, &N));
  PetscCall(MatSetSizes(qn->B, n, n, N, N));
  PetscCall(MatSetUp(qn->B));
  PetscCall(MatLMVMReset(qn->B, PETSC_TRUE));
  PetscCall(MatLMVMSetHistorySize(qn->B, qn->m));
  PetscCall(MatLMVMAllocate(qn->B, snes->vec_sol, snes->vec_func));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBegin;
  PetscCall(MatDestroy(&qn->B));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode SNESDestroy_QN(SNES snes)
{
  PetscFunctionBegin;
  PetscCall(SNESReset_QN(snes));
  PetscCall(PetscFree(snes->data));
  PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", NULL));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode SNESSetFromOptions_QN(SNES snes, PetscOptionItems PetscOptionsObject)
{
  SNES_QN          *qn = (SNES_QN *)snes->data;
  PetscBool         flg;
  SNESLineSearch    linesearch;
  SNESQNRestartType rtype = qn->restart_type;
  SNESQNScaleType   stype = qn->scale_type;
  SNESQNType        qtype = qn->type;

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

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

  PetscCall(PetscOptionsEnum("-snes_qn_type", "Quasi-Newton update type", "", SNESQNTypes, (PetscEnum)qtype, (PetscEnum *)&qtype, &flg));
  if (flg) PetscCall(SNESQNSetType(snes, qtype));
  PetscOptionsHeadEnd();
  PetscCall(SNESQNGetMatrix_Private(snes, &qn->B));
  PetscCall(MatSetFromOptions(qn->B));
  if (!snes->linesearch) {
    PetscCall(SNESGetLineSearch(snes, &linesearch));
    if (!((PetscObject)linesearch)->type_name) {
      if (qn->type == SNES_QN_LBFGS) {
        PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHCP));
      } else if (qn->type == SNES_QN_BROYDEN) {
        PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC));
      } else {
        PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHSECANT));
      }
    }
  }
  if (qn->monflg) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes), &qn->monitor));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

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

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

  Logically Collective

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

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

  Level: intermediate

.seealso: [](ch_snes), `SNES`, `SNESQN`, `SNESQNRestartType`, `SNES_QN_RESTART_NONE`, `SNES_QN_RESTART_POWELL`, `SNES_QN_RESTART_PERIODIC`,
          `SNESQNType`, `SNESQNScaleType`
@*/
PetscErrorCode SNESQNSetRestartType(SNES snes, SNESQNRestartType rtype)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
  PetscTryMethod(snes, "SNESQNSetRestartType_C", (SNES, SNESQNRestartType), (snes, rtype));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  Logically Collective

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

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

  Level: intermediate

.seealso: [](ch_snes), `SNES`, `SNESQN`, `SNESLineSearch`, `SNESQNScaleType`, `SNESSetJacobian()`, `SNESQNType`, `SNESQNRestartType`
@*/
PetscErrorCode SNESQNSetScaleType(SNES snes, SNESQNScaleType stype)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
  PetscTryMethod(snes, "SNESQNSetScaleType_C", (SNES, SNESQNScaleType), (snes, stype));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBegin;
  qn->scale_type = stype;
  if (stype == SNES_QN_SCALE_JACOBIAN) snes->usesksp = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBegin;
  qn->restart_type = rtype;
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  Logically Collective

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

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

  Level: intermediate

.seealso: [](ch_snes), `SNESQN`, `SNES_QN_LBFGS`, `SNES_QN_BROYDEN`, `SNES_QN_BADBROYDEN`, `SNESQNType`,  `SNESQNScaleType`, `TAOLMVM`, `TAOBLMVM`
@*/
PetscErrorCode SNESQNSetType(SNES snes, SNESQNType qtype)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(snes, SNES_CLASSID, 1);
  PetscTryMethod(snes, "SNESQNSetType_C", (SNES, SNESQNType), (snes, qtype));
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

  PetscFunctionBegin;
  qn->type = qtype;
  PetscFunctionReturn(PETSC_SUCCESS);
}

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

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

      Level: beginner

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

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

      Uses left nonlinear preconditioning by default.

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

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

  PetscFunctionBegin;
  snes->ops->setup          = SNESSetUp_QN;
  snes->ops->solve          = SNESSolve_QN;
  snes->ops->destroy        = SNESDestroy_QN;
  snes->ops->setfromoptions = SNESSetFromOptions_QN;
  snes->ops->view           = SNESView_QN;
  snes->ops->reset          = SNESReset_QN;

  snes->npcside = PC_LEFT;

  snes->usesnpc = PETSC_TRUE;
  snes->usesksp = PETSC_FALSE;

  snes->alwayscomputesfinalresidual = PETSC_TRUE;

  PetscCall(SNESParametersInitialize(snes));
  PetscObjectParameterSetDefault(snes, max_funcs, 30000);
  PetscObjectParameterSetDefault(snes, max_its, 10000);

  PetscCall(PetscNew(&qn));
  snes->data       = (void *)qn;
  qn->m            = 10;
  qn->scaling      = 1.0;
  qn->monitor      = NULL;
  qn->monflg       = PETSC_FALSE;
  qn->powell_gamma = 0.9999;
  qn->scale_type   = SNES_QN_SCALE_DEFAULT;
  qn->restart_type = SNES_QN_RESTART_DEFAULT;
  qn->type         = SNES_QN_LBFGS;

  PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetScaleType_C", SNESQNSetScaleType_QN));
  PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetRestartType_C", SNESQNSetRestartType_QN));
  PetscCall(PetscObjectComposeFunction((PetscObject)snes, "SNESQNSetType_C", SNESQNSetType_QN));
  PetscFunctionReturn(PETSC_SUCCESS);
}