#include <../src/ts/impls/implicit/glle/glle.h> /*I   "petscts.h"   I*/
#include <petscdm.h>
#include <petscblaslapack.h>

static const char       *TSGLLEErrorDirections[] = {"FORWARD", "BACKWARD", "TSGLLEErrorDirection", "TSGLLEERROR_", NULL};
static PetscFunctionList TSGLLEList;
static PetscFunctionList TSGLLEAcceptList;
static PetscBool         TSGLLEPackageInitialized;
static PetscBool         TSGLLERegisterAllCalled;

/* This function is pure */
static PetscScalar Factorial(PetscInt n)
{
  PetscInt i;
  if (n < 12) { /* Can compute with 32-bit integers */
    PetscInt f = 1;
    for (i = 2; i <= n; i++) f *= i;
    return (PetscScalar)f;
  } else {
    PetscScalar f = 1.;
    for (i = 2; i <= n; i++) f *= (PetscScalar)i;
    return f;
  }
}

/* This function is pure */
static PetscScalar CPowF(PetscScalar c, PetscInt p)
{
  return PetscPowRealInt(PetscRealPart(c), p) / Factorial(p);
}

static PetscErrorCode TSGLLEGetVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage)
{
  TS_GLLE *gl = (TS_GLLE *)ts->data;

  PetscFunctionBegin;
  if (Z) {
    if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSGLLE_Z", Z));
    else *Z = gl->Z;
  }
  if (Ydotstage) {
    if (dm && dm != ts->dm) PetscCall(DMGetNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage));
    else *Ydotstage = gl->Ydot[gl->stage];
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLERestoreVecs(TS ts, DM dm, Vec *Z, Vec *Ydotstage)
{
  PetscFunctionBegin;
  if (Z) {
    if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLLE_Z", Z));
  }
  if (Ydotstage) {
    if (dm && dm != ts->dm) PetscCall(DMRestoreNamedGlobalVector(dm, "TSGLLE_Ydot", Ydotstage));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMCoarsenHook_TSGLLE(DM fine, DM coarse, PetscCtx ctx)
{
  PetscFunctionBegin;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMRestrictHook_TSGLLE(DM fine, Mat restrct, Vec rscale, Mat inject, DM coarse, PetscCtx ctx)
{
  TS  ts = (TS)ctx;
  Vec Ydot, Ydot_c;

  PetscFunctionBegin;
  PetscCall(TSGLLEGetVecs(ts, fine, NULL, &Ydot));
  PetscCall(TSGLLEGetVecs(ts, coarse, NULL, &Ydot_c));
  PetscCall(MatRestrict(restrct, Ydot, Ydot_c));
  PetscCall(VecPointwiseMult(Ydot_c, rscale, Ydot_c));
  PetscCall(TSGLLERestoreVecs(ts, fine, NULL, &Ydot));
  PetscCall(TSGLLERestoreVecs(ts, coarse, NULL, &Ydot_c));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMSubDomainHook_TSGLLE(DM dm, DM subdm, PetscCtx ctx)
{
  PetscFunctionBegin;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode DMSubDomainRestrictHook_TSGLLE(DM dm, VecScatter gscat, VecScatter lscat, DM subdm, PetscCtx ctx)
{
  TS  ts = (TS)ctx;
  Vec Ydot, Ydot_s;

  PetscFunctionBegin;
  PetscCall(TSGLLEGetVecs(ts, dm, NULL, &Ydot));
  PetscCall(TSGLLEGetVecs(ts, subdm, NULL, &Ydot_s));

  PetscCall(VecScatterBegin(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD));
  PetscCall(VecScatterEnd(gscat, Ydot, Ydot_s, INSERT_VALUES, SCATTER_FORWARD));

  PetscCall(TSGLLERestoreVecs(ts, dm, NULL, &Ydot));
  PetscCall(TSGLLERestoreVecs(ts, subdm, NULL, &Ydot_s));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLESchemeCreate(PetscInt p, PetscInt q, PetscInt r, PetscInt s, const PetscScalar *c, const PetscScalar *a, const PetscScalar *b, const PetscScalar *u, const PetscScalar *v, TSGLLEScheme *inscheme)
{
  TSGLLEScheme scheme;
  PetscInt     j;

  PetscFunctionBegin;
  PetscCheck(p >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Scheme order must be positive");
  PetscCheck(r >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "At least one item must be carried between steps");
  PetscCheck(s >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "At least one stage is required");
  PetscAssertPointer(inscheme, 10);
  *inscheme = NULL;
  PetscCall(PetscNew(&scheme));
  scheme->p = p;
  scheme->q = q;
  scheme->r = r;
  scheme->s = s;

  PetscCall(PetscMalloc5(s, &scheme->c, s * s, &scheme->a, r * s, &scheme->b, r * s, &scheme->u, r * r, &scheme->v));
  PetscCall(PetscArraycpy(scheme->c, c, s));
  for (j = 0; j < s * s; j++) scheme->a[j] = (PetscAbsScalar(a[j]) < 1e-12) ? 0 : a[j];
  for (j = 0; j < r * s; j++) scheme->b[j] = (PetscAbsScalar(b[j]) < 1e-12) ? 0 : b[j];
  for (j = 0; j < s * r; j++) scheme->u[j] = (PetscAbsScalar(u[j]) < 1e-12) ? 0 : u[j];
  for (j = 0; j < r * r; j++) scheme->v[j] = (PetscAbsScalar(v[j]) < 1e-12) ? 0 : v[j];

  PetscCall(PetscMalloc6(r, &scheme->alpha, r, &scheme->beta, r, &scheme->gamma, 3 * s, &scheme->phi, 3 * r, &scheme->psi, r, &scheme->stage_error));
  {
    PetscInt     i, j, k, ss = s + 2;
    PetscBLASInt m, n, one = 1, *ipiv, lwork, info, ldb;
    PetscReal    rcond, *sing, *workreal;
    PetscScalar *ImV, *H, *bmat, *workscalar, *c = scheme->c, *a = scheme->a, *b = scheme->b, *u = scheme->u, *v = scheme->v;
    PetscBLASInt rank;

    PetscCall(PetscBLASIntCast(4 * ((s + 3) * 3 + 3), &lwork));
    PetscCall(PetscMalloc7(PetscSqr(r), &ImV, 3 * s, &H, 3 * ss, &bmat, lwork, &workscalar, 5 * (3 + r), &workreal, r + s, &sing, r + s, &ipiv));

    /* column-major input */
    for (i = 0; i < r - 1; i++) {
      for (j = 0; j < r - 1; j++) ImV[i + j * r] = 1.0 * (i == j) - v[(i + 1) * r + j + 1];
    }
    /* Build right-hand side for alpha (tp - glm.B(2:end,:)*(glm.c.^(p)./factorial(p))) */
    for (i = 1; i < r; i++) {
      scheme->alpha[i] = 1. / Factorial(p + 1 - i);
      for (j = 0; j < s; j++) scheme->alpha[i] -= b[i * s + j] * CPowF(c[j], p);
    }
    PetscCall(PetscBLASIntCast(r - 1, &m));
    PetscCall(PetscBLASIntCast(r, &n));
    PetscCallBLAS("LAPACKgesv", LAPACKgesv_(&m, &one, ImV, &n, ipiv, scheme->alpha + 1, &n, &info));
    PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GESV");
    PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_MAT_LU_ZRPVT, "Bad LU factorization");

    /* Build right-hand side for beta (tp1 - glm.B(2:end,:)*(glm.c.^(p+1)./factorial(p+1)) - e.alpha) */
    for (i = 1; i < r; i++) {
      scheme->beta[i] = 1. / Factorial(p + 2 - i) - scheme->alpha[i];
      for (j = 0; j < s; j++) scheme->beta[i] -= b[i * s + j] * CPowF(c[j], p + 1);
    }
    PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->beta + 1, &n, &info));
    PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GETRS");
    PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Should not happen");

    /* Build stage_error vector
           xi = glm.c.^(p+1)/factorial(p+1) - glm.A*glm.c.^p/factorial(p) + glm.U(:,2:end)*e.alpha;
    */
    for (i = 0; i < s; i++) {
      scheme->stage_error[i] = CPowF(c[i], p + 1);
      for (j = 0; j < s; j++) scheme->stage_error[i] -= a[i * s + j] * CPowF(c[j], p);
      for (j = 1; j < r; j++) scheme->stage_error[i] += u[i * r + j] * scheme->alpha[j];
    }

    /* alpha[0] (epsilon in B,J,W 2007)
           epsilon = 1/factorial(p+1) - B(1,:)*c.^p/factorial(p) + V(1,2:end)*e.alpha;
    */
    scheme->alpha[0] = 1. / Factorial(p + 1);
    for (j = 0; j < s; j++) scheme->alpha[0] -= b[0 * s + j] * CPowF(c[j], p);
    for (j = 1; j < r; j++) scheme->alpha[0] += v[0 * r + j] * scheme->alpha[j];

    /* right-hand side for gamma (glm.B(2:end,:)*e.xi - e.epsilon*eye(s-1,1)) */
    for (i = 1; i < r; i++) {
      scheme->gamma[i] = (i == 1 ? -1. : 0) * scheme->alpha[0];
      for (j = 0; j < s; j++) scheme->gamma[i] += b[i * s + j] * scheme->stage_error[j];
    }
    PetscCallBLAS("LAPACKgetrs", LAPACKgetrs_("No transpose", &m, &one, ImV, &n, ipiv, scheme->gamma + 1, &n, &info));
    PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GETRS");
    PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Should not happen");

    /* beta[0] (rho in B,J,W 2007)
        e.rho = 1/factorial(p+2) - glm.B(1,:)*glm.c.^(p+1)/factorial(p+1) ...
            + glm.V(1,2:end)*e.beta;% - e.epsilon;
    % Note: The paper (B,J,W 2007) includes the last term in their definition
    * */
    scheme->beta[0] = 1. / Factorial(p + 2);
    for (j = 0; j < s; j++) scheme->beta[0] -= b[0 * s + j] * CPowF(c[j], p + 1);
    for (j = 1; j < r; j++) scheme->beta[0] += v[0 * r + j] * scheme->beta[j];

    /* gamma[0] (sigma in B,J,W 2007)
    *   e.sigma = glm.B(1,:)*e.xi + glm.V(1,2:end)*e.gamma;
    * */
    scheme->gamma[0] = 0.0;
    for (j = 0; j < s; j++) scheme->gamma[0] += b[0 * s + j] * scheme->stage_error[j];
    for (j = 1; j < r; j++) scheme->gamma[0] += v[0 * s + j] * scheme->gamma[j];

    /* Assemble H
    *    % " PetscInt_FMT "etermine the error estimators phi
       H = [[cpow(glm.c,p) + C*e.alpha] [cpow(glm.c,p+1) + C*e.beta] ...
               [e.xi - C*(e.gamma + 0*e.epsilon*eye(s-1,1))]]';
    % Paper has formula above without the 0, but that term must be left
    % out to satisfy the conditions they propose and to make the
    % example schemes work
    e.H = H;
    e.phi = (H \ [1 0 0;1 1 0;0 0 -1])';
    e.psi = -e.phi*C;
    * */
    for (j = 0; j < s; j++) {
      H[0 + j * 3] = CPowF(c[j], p);
      H[1 + j * 3] = CPowF(c[j], p + 1);
      H[2 + j * 3] = scheme->stage_error[j];
      for (k = 1; k < r; k++) {
        H[0 + j * 3] += CPowF(c[j], k - 1) * scheme->alpha[k];
        H[1 + j * 3] += CPowF(c[j], k - 1) * scheme->beta[k];
        H[2 + j * 3] -= CPowF(c[j], k - 1) * scheme->gamma[k];
      }
    }
    bmat[0 + 0 * ss] = 1.;
    bmat[0 + 1 * ss] = 0.;
    bmat[0 + 2 * ss] = 0.;
    bmat[1 + 0 * ss] = 1.;
    bmat[1 + 1 * ss] = 1.;
    bmat[1 + 2 * ss] = 0.;
    bmat[2 + 0 * ss] = 0.;
    bmat[2 + 1 * ss] = 0.;
    bmat[2 + 2 * ss] = -1.;
    m                = 3;
    PetscCall(PetscBLASIntCast(s, &n));
    PetscCall(PetscBLASIntCast(ss, &ldb));
    rcond = 1e-12;
#if defined(PETSC_USE_COMPLEX)
    /* ZGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, RWORK, INFO) */
    PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, workreal, &info));
#else
    /* DGELSS( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO) */
    PetscCallBLAS("LAPACKgelss", LAPACKgelss_(&m, &n, &m, H, &m, bmat, &ldb, sing, &rcond, &rank, workscalar, &lwork, &info));
#endif
    PetscCheck(info >= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "Bad argument to GELSS");
    PetscCheck(info <= 0, PETSC_COMM_SELF, PETSC_ERR_LIB, "SVD failed to converge");

    for (j = 0; j < 3; j++) {
      for (k = 0; k < s; k++) scheme->phi[k + j * s] = bmat[k + j * ss];
    }

    /* the other part of the error estimator, psi in B,J,W 2007 */
    scheme->psi[0 * r + 0] = 0.;
    scheme->psi[1 * r + 0] = 0.;
    scheme->psi[2 * r + 0] = 0.;
    for (j = 1; j < r; j++) {
      scheme->psi[0 * r + j] = 0.;
      scheme->psi[1 * r + j] = 0.;
      scheme->psi[2 * r + j] = 0.;
      for (k = 0; k < s; k++) {
        scheme->psi[0 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[0 * s + k];
        scheme->psi[1 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[1 * s + k];
        scheme->psi[2 * r + j] -= CPowF(c[k], j - 1) * scheme->phi[2 * s + k];
      }
    }
    PetscCall(PetscFree7(ImV, H, bmat, workscalar, workreal, sing, ipiv));
  }
  /* Check which properties are satisfied */
  scheme->stiffly_accurate = PETSC_TRUE;
  if (scheme->c[s - 1] != 1.) scheme->stiffly_accurate = PETSC_FALSE;
  for (j = 0; j < s; j++)
    if (a[(s - 1) * s + j] != b[j]) scheme->stiffly_accurate = PETSC_FALSE;
  for (j = 0; j < r; j++)
    if (u[(s - 1) * r + j] != v[j]) scheme->stiffly_accurate = PETSC_FALSE;
  scheme->fsal = scheme->stiffly_accurate; /* FSAL is stronger */
  for (j = 0; j < s - 1; j++)
    if (r > 1 && b[1 * s + j] != 0.) scheme->fsal = PETSC_FALSE;
  if (b[1 * s + r - 1] != 1.) scheme->fsal = PETSC_FALSE;
  for (j = 0; j < r; j++)
    if (r > 1 && v[1 * r + j] != 0.) scheme->fsal = PETSC_FALSE;

  *inscheme = scheme;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLESchemeDestroy(TSGLLEScheme sc)
{
  PetscFunctionBegin;
  PetscCall(PetscFree5(sc->c, sc->a, sc->b, sc->u, sc->v));
  PetscCall(PetscFree6(sc->alpha, sc->beta, sc->gamma, sc->phi, sc->psi, sc->stage_error));
  PetscCall(PetscFree(sc));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLEDestroy_Default(TS_GLLE *gl)
{
  PetscInt i;

  PetscFunctionBegin;
  for (i = 0; i < gl->nschemes; i++) {
    if (gl->schemes[i]) PetscCall(TSGLLESchemeDestroy(gl->schemes[i]));
  }
  PetscCall(PetscFree(gl->schemes));
  gl->nschemes = 0;
  PetscCall(PetscMemzero(gl->type_name, sizeof(gl->type_name)));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLEViewTable_Private(PetscViewer viewer, PetscInt m, PetscInt n, const PetscScalar a[], const char name[])
{
  PetscBool isascii;
  PetscInt  i, j;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
  if (isascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "%30s = [", name));
    for (i = 0; i < m; i++) {
      if (i) PetscCall(PetscViewerASCIIPrintf(viewer, "%30s   [", ""));
      PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
      for (j = 0; j < n; j++) PetscCall(PetscViewerASCIIPrintf(viewer, " %12.8g", (double)PetscRealPart(a[i * n + j])));
      PetscCall(PetscViewerASCIIPrintf(viewer, "]\n"));
      PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
    }
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLESchemeView(TSGLLEScheme sc, PetscBool view_details, PetscViewer viewer)
{
  PetscBool isascii;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
  PetscCheck(isascii, PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type %s not supported", ((PetscObject)viewer)->type_name);
  PetscCall(PetscViewerASCIIPrintf(viewer, "GL scheme p,q,r,s = %" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "\n", sc->p, sc->q, sc->r, sc->s));
  PetscCall(PetscViewerASCIIPushTab(viewer));
  PetscCall(PetscViewerASCIIPrintf(viewer, "Stiffly accurate: %s,  FSAL: %s\n", sc->stiffly_accurate ? "yes" : "no", sc->fsal ? "yes" : "no"));
  PetscCall(PetscViewerASCIIPrintf(viewer, "Leading error constants: %10.3e  %10.3e  %10.3e\n", (double)PetscRealPart(sc->alpha[0]), (double)PetscRealPart(sc->beta[0]), (double)PetscRealPart(sc->gamma[0])));
  PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->s, sc->c, "Abscissas c"));
  if (view_details) {
    PetscCall(TSGLLEViewTable_Private(viewer, sc->s, sc->s, sc->a, "A"));
    PetscCall(TSGLLEViewTable_Private(viewer, sc->r, sc->s, sc->b, "B"));
    PetscCall(TSGLLEViewTable_Private(viewer, sc->s, sc->r, sc->u, "U"));
    PetscCall(TSGLLEViewTable_Private(viewer, sc->r, sc->r, sc->v, "V"));

    PetscCall(TSGLLEViewTable_Private(viewer, 3, sc->s, sc->phi, "Error estimate phi"));
    PetscCall(TSGLLEViewTable_Private(viewer, 3, sc->r, sc->psi, "Error estimate psi"));
    PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->alpha, "Modify alpha"));
    PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->beta, "Modify beta"));
    PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->r, sc->gamma, "Modify gamma"));
    PetscCall(TSGLLEViewTable_Private(viewer, 1, sc->s, sc->stage_error, "Stage error xi"));
  }
  PetscCall(PetscViewerASCIIPopTab(viewer));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLEEstimateHigherMoments_Default(TSGLLEScheme sc, PetscReal h, Vec Ydot[], Vec Xold[], Vec hm[])
{
  PetscInt i;

  PetscFunctionBegin;
  PetscCheck(sc->r <= 64 && sc->s <= 64, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Ridiculous number of stages or items passed between stages");
  /* build error vectors*/
  for (i = 0; i < 3; i++) {
    PetscScalar phih[64];
    PetscInt    j;
    for (j = 0; j < sc->s; j++) phih[j] = sc->phi[i * sc->s + j] * h;
    PetscCall(VecZeroEntries(hm[i]));
    PetscCall(VecMAXPY(hm[i], sc->s, phih, Ydot));
    PetscCall(VecMAXPY(hm[i], sc->r, &sc->psi[i * sc->r], Xold));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLECompleteStep_Rescale(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[])
{
  PetscScalar brow[32], vrow[32];
  PetscInt    i, j, r, s;

  PetscFunctionBegin;
  /* Build the new solution from (X,Ydot) */
  r = sc->r;
  s = sc->s;
  for (i = 0; i < r; i++) {
    PetscCall(VecZeroEntries(X[i]));
    for (j = 0; j < s; j++) brow[j] = h * sc->b[i * s + j];
    PetscCall(VecMAXPY(X[i], s, brow, Ydot));
    for (j = 0; j < r; j++) vrow[j] = sc->v[i * r + j];
    PetscCall(VecMAXPY(X[i], r, vrow, Xold));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLECompleteStep_RescaleAndModify(TSGLLEScheme sc, PetscReal h, TSGLLEScheme next_sc, PetscReal next_h, Vec Ydot[], Vec Xold[], Vec X[])
{
  PetscScalar brow[32], vrow[32];
  PetscReal   ratio;
  PetscInt    i, j, p, r, s;

  PetscFunctionBegin;
  /* Build the new solution from (X,Ydot) */
  p     = sc->p;
  r     = sc->r;
  s     = sc->s;
  ratio = next_h / h;
  for (i = 0; i < r; i++) {
    PetscCall(VecZeroEntries(X[i]));
    for (j = 0; j < s; j++) {
      brow[j] = h * (PetscPowRealInt(ratio, i) * sc->b[i * s + j] + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 1)) * (+sc->alpha[i] * sc->phi[0 * s + j]) + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 2)) * (+sc->beta[i] * sc->phi[1 * s + j] + sc->gamma[i] * sc->phi[2 * s + j]));
    }
    PetscCall(VecMAXPY(X[i], s, brow, Ydot));
    for (j = 0; j < r; j++) {
      vrow[j] = (PetscPowRealInt(ratio, i) * sc->v[i * r + j] + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 1)) * (+sc->alpha[i] * sc->psi[0 * r + j]) + (PetscPowRealInt(ratio, i) - PetscPowRealInt(ratio, p + 2)) * (+sc->beta[i] * sc->psi[1 * r + j] + sc->gamma[i] * sc->psi[2 * r + j]));
    }
    PetscCall(VecMAXPY(X[i], r, vrow, Xold));
  }
  if (r < next_sc->r) {
    PetscCheck(r + 1 == next_sc->r, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Cannot accommodate jump in r greater than 1");
    PetscCall(VecZeroEntries(X[r]));
    for (j = 0; j < s; j++) brow[j] = h * PetscPowRealInt(ratio, p + 1) * sc->phi[0 * s + j];
    PetscCall(VecMAXPY(X[r], s, brow, Ydot));
    for (j = 0; j < r; j++) vrow[j] = PetscPowRealInt(ratio, p + 1) * sc->psi[0 * r + j];
    PetscCall(VecMAXPY(X[r], r, vrow, Xold));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLECreate_IRKS(TS ts)
{
  TS_GLLE *gl = (TS_GLLE *)ts->data;

  PetscFunctionBegin;
  gl->Destroy               = TSGLLEDestroy_Default;
  gl->EstimateHigherMoments = TSGLLEEstimateHigherMoments_Default;
  gl->CompleteStep          = TSGLLECompleteStep_RescaleAndModify;
  PetscCall(PetscMalloc1(10, &gl->schemes));
  gl->nschemes = 0;

  {
    /* p=1,q=1, r=s=2, A- and L-stable with error estimates of order 2 and 3
    * Listed in Butcher & Podhaisky 2006. On error estimation in general linear methods for stiff ODE.
    * irks(0.3,0,[.3,1],[1],1)
    * Note: can be made to have classical order (not stage order) 2 by replacing 0.3 with 1-sqrt(1/2)
    * but doing so would sacrifice the error estimator.
    */
    const PetscScalar c[2]    = {3. / 10., 1.};
    const PetscScalar a[2][2] = {
      {3. / 10., 0       },
      {7. / 10., 3. / 10.}
    };
    const PetscScalar b[2][2] = {
      {7. / 10., 3. / 10.},
      {0,        1       }
    };
    const PetscScalar u[2][2] = {
      {1, 0},
      {1, 0}
    };
    const PetscScalar v[2][2] = {
      {1, 0},
      {0, 0}
    };
    PetscCall(TSGLLESchemeCreate(1, 1, 2, 2, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]));
  }

  {
    /* p=q=2, r=s=3: irks(4/9,0,[1:3]/3,[0.33852],1) */
    /* http://www.math.auckland.ac.nz/~hpod/atlas/i2a.html */
    const PetscScalar c[3]    = {1. / 3., 2. / 3., 1};
    const PetscScalar a[3][3] = {
      {4. / 9.,              0,                     0      },
      {1.03750643704090e+00, 4. / 9.,               0      },
      {7.67024779410304e-01, -3.81140216918943e-01, 4. / 9.}
    };
    const PetscScalar b[3][3] = {
      {0.767024779410304,  -0.381140216918943, 4. / 9.          },
      {0.000000000000000,  0.000000000000000,  1.000000000000000},
      {-2.075048385225385, 0.621728385225383,  1.277197204924873}
    };
    const PetscScalar u[3][3] = {
      {1.0000000000000000, -0.1111111111111109, -0.0925925925925922},
      {1.0000000000000000, -0.8152842148186744, -0.4199095530877056},
      {1.0000000000000000, 0.1696709930641948,  0.0539741070314165 }
    };
    const PetscScalar v[3][3] = {
      {1.0000000000000000, 0.1696709930641948, 0.0539741070314165},
      {0.000000000000000,  0.000000000000000,  0.000000000000000 },
      {0.000000000000000,  0.176122795075129,  0.000000000000000 }
    };
    PetscCall(TSGLLESchemeCreate(2, 2, 3, 3, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]));
  }
  {
    /* p=q=3, r=s=4: irks(9/40,0,[1:4]/4,[0.3312 1.0050],[0.49541 1;1 0]) */
    const PetscScalar c[4]    = {0.25, 0.5, 0.75, 1.0};
    const PetscScalar a[4][4] = {
      {9. / 40.,             0,                     0,                 0       },
      {2.11286958887701e-01, 9. / 40.,              0,                 0       },
      {9.46338294287584e-01, -3.42942861246094e-01, 9. / 40.,          0       },
      {0.521490453970721,    -0.662474225622980,    0.490476425459734, 9. / 40.}
    };
    const PetscScalar b[4][4] = {
      {0.521490453970721,  -0.662474225622980, 0.490476425459734,  9. / 40.         },
      {0.000000000000000,  0.000000000000000,  0.000000000000000,  1.000000000000000},
      {-0.084677029310348, 1.390757514776085,  -1.568157386206001, 2.023192696767826},
      {0.465383797936408,  1.478273530625148,  -1.930836081010182, 1.644872111193354}
    };
    const PetscScalar u[4][4] = {
      {1.00000000000000000, 0.02500000000001035,  -0.02499999999999053, -0.00442708333332865},
      {1.00000000000000000, 0.06371304111232945,  -0.04032173972189845, -0.01389438413189452},
      {1.00000000000000000, -0.07839543304147778, 0.04738685705116663,  0.02032603595928376 },
      {1.00000000000000000, 0.42550734619251651,  0.10800718022400080,  -0.01726712647760034}
    };
    const PetscScalar v[4][4] = {
      {1.00000000000000000, 0.42550734619251651, 0.10800718022400080, -0.01726712647760034},
      {0.000000000000000,   0.000000000000000,   0.000000000000000,   0.000000000000000   },
      {0.000000000000000,   -1.761115796027561,  -0.521284157173780,  0.258249384305463   },
      {0.000000000000000,   -1.657693358744728,  -1.052227765232394,  0.521284157173780   }
    };
    PetscCall(TSGLLESchemeCreate(3, 3, 4, 4, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]));
  }
  {
    /* p=q=4, r=s=5:
          irks(3/11,0,[1:5]/5, [0.1715   -0.1238    0.6617],...
          [ -0.0812    0.4079    1.0000
             1.0000         0         0
             0.8270    1.0000         0])
    */
    const PetscScalar c[5]    = {0.2, 0.4, 0.6, 0.8, 1.0};
    const PetscScalar a[5][5] = {
      {2.72727272727352e-01,  0.00000000000000e+00,  0.00000000000000e+00,  0.00000000000000e+00, 0.00000000000000e+00},
      {-1.03980153733431e-01, 2.72727272727405e-01,  0.00000000000000e+00,  0.00000000000000e+00, 0.00000000000000e+00},
      {-1.58615400341492e+00, 7.44168951881122e-01,  2.72727272727309e-01,  0.00000000000000e+00, 0.00000000000000e+00},
      {-8.73658042865628e-01, 5.37884671894595e-01,  -1.63298538799523e-01, 2.72727272726996e-01, 0.00000000000000e+00},
      {2.95489397443992e-01,  -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00, 2.72727272727288e-01}
    };
    const PetscScalar b[5][5] = {
      {2.95489397443992e-01,  -1.18481693910097e+00, -6.68029812659953e-01, 1.00716687860943e+00,  2.72727272727288e-01},
      {0.00000000000000e+00,  1.11022302462516e-16,  -2.22044604925031e-16, 0.00000000000000e+00,  1.00000000000000e+00},
      {-4.05882503986005e+00, -4.00924006567769e+00, -1.38930610972481e+00, 4.45223930308488e+00,  6.32331093108427e-01},
      {8.35690179937017e+00,  -2.26640927349732e+00, 6.86647884973826e+00,  -5.22595158025740e+00, 4.50893068837431e+00},
      {1.27656267027479e+01,  2.80882153840821e+00,  8.91173096522890e+00,  -1.07936444078906e+01, 4.82534148988854e+00}
    };
    const PetscScalar u[5][5] = {
      {1.00000000000000e+00, -7.27272727273551e-02, -3.45454545454419e-02, -4.12121212119565e-03, -2.96969696964014e-04},
      {1.00000000000000e+00, 2.31252881006154e-01,  -8.29487834416481e-03, -9.07191207681020e-03, -1.70378403743473e-03},
      {1.00000000000000e+00, 1.16925777880663e+00,  3.59268562942635e-02,  -4.09013451730615e-02, -1.02411119670164e-02},
      {1.00000000000000e+00, 1.02634463704356e+00,  1.59375044913405e-01,  1.89673015035370e-03,  -4.89987231897569e-03},
      {1.00000000000000e+00, 1.27746320298021e+00,  2.37186008132728e-01,  -8.28694373940065e-02, -5.34396510196430e-02}
    };
    const PetscScalar v[5][5] = {
      {1.00000000000000e+00, 1.27746320298021e+00,  2.37186008132728e-01,  -8.28694373940065e-02, -5.34396510196430e-02},
      {0.00000000000000e+00, -1.77635683940025e-15, -1.99840144432528e-15, -9.99200722162641e-16, -3.33066907387547e-16},
      {0.00000000000000e+00, 4.37280081906924e+00,  5.49221645016377e-02,  -8.88913177394943e-02, 1.12879077989154e-01 },
      {0.00000000000000e+00, -1.22399504837280e+01, -5.21287338448645e+00, -8.03952325565291e-01, 4.60298678047147e-01 },
      {0.00000000000000e+00, -1.85178762883829e+01, -5.21411849862624e+00, -1.04283436528809e+00, 7.49030161063651e-01 }
    };
    PetscCall(TSGLLESchemeCreate(4, 4, 5, 5, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]));
  }
  {
    /* p=q=5, r=s=6;
       irks(1/3,0,[1:6]/6,...
          [-0.0489    0.4228   -0.8814    0.9021],...
          [-0.3474   -0.6617    0.6294    0.2129
            0.0044   -0.4256   -0.1427   -0.8936
           -0.8267    0.4821    0.1371   -0.2557
           -0.4426   -0.3855   -0.7514    0.3014])
    */
    const PetscScalar c[6]    = {1. / 6, 2. / 6, 3. / 6, 4. / 6, 5. / 6, 1.};
    const PetscScalar a[6][6] = {
      {3.33333333333940e-01,  0,                     0,                     0,                     0,                    0                   },
      {-8.64423857333350e-02, 3.33333333332888e-01,  0,                     0,                     0,                    0                   },
      {-2.16850174258252e+00, -2.23619072028839e+00, 3.33333333335204e-01,  0,                     0,                    0                   },
      {-4.73160970138997e+00, -3.89265344629268e+00, -2.76318716520933e-01, 3.33333333335759e-01,  0,                    0                   },
      {-6.75187540297338e+00, -7.90756533769377e+00, 7.90245051802259e-01,  -4.48352364517632e-01, 3.33333333328483e-01, 0                   },
      {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00,  -2.23969848002481e+00, 6.62807710124007e-01, 3.33333333335440e-01}
    };
    const PetscScalar b[6][6] = {
      {-4.26488287921548e+00, -1.19320395589302e+01, 3.38924509887755e+00,  -2.23969848002481e+00, 6.62807710124007e-01,  3.33333333335440e-01 },
      {-8.88178419700125e-16, 4.44089209850063e-16,  -1.54737334057131e-15, -8.88178419700125e-16, 0.00000000000000e+00,  1.00000000000001e+00 },
      {-2.87780425770651e+01, -1.13520448264971e+01, 2.62002318943161e+01,  2.56943874812797e+01,  -3.06702268304488e+01, 6.68067773510103e+00 },
      {5.47971245256474e+01,  6.80366875868284e+01,  -6.50952588861999e+01, -8.28643975339097e+01, 8.17416943896414e+01,  -1.17819043489036e+01},
      {-2.33332114788869e+02, 6.12942539462634e+01,  -4.91850135865944e+01, 1.82716844135480e+02,  -1.29788173979395e+02, 3.09968095651099e+01 },
      {-1.72049132343751e+02, 8.60194713593999e+00,  7.98154219170200e-01,  1.50371386053218e+02,  -1.18515423962066e+02, 2.50898277784663e+01 }
    };
    const PetscScalar u[6][6] = {
      {1.00000000000000e+00, -1.66666666666870e-01, -4.16666666664335e-02, -3.85802469124815e-03, -2.25051440302250e-04, -9.64506172339142e-06},
      {1.00000000000000e+00, 8.64423857327162e-02,  -4.11484912671353e-02, -1.11450903217645e-02, -1.47651050487126e-03, -1.34395070766826e-04},
      {1.00000000000000e+00, 4.57135912953434e+00,  1.06514719719137e+00,  1.33517564218007e-01,  1.11365952968659e-02,  6.12382756769504e-04 },
      {1.00000000000000e+00, 9.23391519753404e+00,  2.22431212392095e+00,  2.91823807741891e-01,  2.52058456411084e-02,  1.22800542949647e-03 },
      {1.00000000000000e+00, 1.48175480533865e+01,  3.73439117461835e+00,  5.14648336541804e-01,  4.76430038853402e-02,  2.56798515502156e-03 },
      {1.00000000000000e+00, 1.50512347758335e+01,  4.10099701165164e+00,  5.66039141003603e-01,  3.91213893800891e-02,  -2.99136269067853e-03}
    };
    const PetscScalar v[6][6] = {
      {1.00000000000000e+00, 1.50512347758335e+01,  4.10099701165164e+00,  5.66039141003603e-01,  3.91213893800891e-02,  -2.99136269067853e-03},
      {0.00000000000000e+00, -4.88498130835069e-15, -6.43929354282591e-15, -3.55271367880050e-15, -1.22124532708767e-15, -3.12250225675825e-16},
      {0.00000000000000e+00, 1.22250171233141e+01,  -1.77150760606169e+00, 3.54516769879390e-01,  6.22298845883398e-01,  2.31647447450276e-01 },
      {0.00000000000000e+00, -4.48339457331040e+01, -3.57363126641880e-01, 5.18750173123425e-01,  6.55727990241799e-02,  1.63175368287079e-01 },
      {0.00000000000000e+00, 1.37297394708005e+02,  -1.60145272991317e+00, -5.05319555199441e+00, 1.55328940390990e-01,  9.16629423682464e-01 },
      {0.00000000000000e+00, 1.05703241119022e+02,  -1.16610260983038e+00, -2.99767252773859e+00, -1.13472315553890e-01, 1.09742849254729e+00 }
    };
    PetscCall(TSGLLESchemeCreate(5, 5, 6, 6, c, *a, *b, *u, *v, &gl->schemes[gl->nschemes++]));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  TSGLLESetType - sets the class of general linear method, `TSGLLE` to use for time-stepping

  Collective

  Input Parameters:
+ ts   - the `TS` context
- type - a method

  Options Database Key:
. -ts_gl_type <type> - sets the method, use -help for a list of available method (e.g. irks)

  Level: intermediate

  Notes:
  See "petsc/include/petscts.h" for available methods (for instance)
.    TSGLLE_IRKS - Diagonally implicit methods with inherent Runge-Kutta stability (for stiff problems)

  Normally, it is best to use the `TSSetFromOptions()` command and then set the `TSGLLE` type
  from the options database rather than by using this routine.  Using the options database
  provides the user with maximum flexibility in evaluating the many different solvers.  The
  `TSGLLESetType()` routine is provided for those situations where it is necessary to set the
  timestepping solver independently of the command line or options database.  This might be the
  case, for example, when the choice of solver changes during the execution of the program, and
  the user's application is taking responsibility for choosing the appropriate method.

.seealso: [](ch_ts), `TS`, `TSGLLEType`, `TSGLLE`
@*/
PetscErrorCode TSGLLESetType(TS ts, TSGLLEType type)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
  PetscAssertPointer(type, 2);
  PetscTryMethod(ts, "TSGLLESetType_C", (TS, TSGLLEType), (ts, type));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  TSGLLESetAcceptType - sets the acceptance test for `TSGLLE`

  Logically Collective

  Input Parameters:
+ ts   - the `TS` context
- type - the type

  Options Database Key:
. -ts_gl_accept_type <type> - sets the method used to determine whether to accept or reject a step

  Level: intermediate

  Notes:
  Time integrators that need to control error must have the option to reject a time step based
  on local error estimates. This function allows different schemes to be set.

.seealso: [](ch_ts), `TS`, `TSGLLE`, `TSGLLEAcceptRegister()`, `TSGLLEAdapt`
@*/
PetscErrorCode TSGLLESetAcceptType(TS ts, TSGLLEAcceptType type)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
  PetscAssertPointer(type, 2);
  PetscTryMethod(ts, "TSGLLESetAcceptType_C", (TS, TSGLLEAcceptType), (ts, type));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@
  TSGLLEGetAdapt - gets the `TSGLLEAdapt` object from the `TS`

  Not Collective

  Input Parameter:
. ts - the `TS` context

  Output Parameter:
. adapt - the `TSGLLEAdapt` context

  Level: advanced

  Note:
  This allows the user set options on the `TSGLLEAdapt` object. Usually it is better to do this
  using the options database, so this function is rarely needed.

.seealso: [](ch_ts), `TS`, `TSGLLE`, `TSGLLEAdapt`, `TSGLLEAdaptRegister()`
@*/
PetscErrorCode TSGLLEGetAdapt(TS ts, TSGLLEAdapt *adapt)
{
  PetscFunctionBegin;
  PetscValidHeaderSpecific(ts, TS_CLASSID, 1);
  PetscAssertPointer(adapt, 2);
  PetscUseMethod(ts, "TSGLLEGetAdapt_C", (TS, TSGLLEAdapt *), (ts, adapt));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLEAccept_Always(TS ts, PetscReal tleft, PetscReal h, const PetscReal enorms[], PetscBool *accept)
{
  PetscFunctionBegin;
  *accept = PETSC_TRUE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLEUpdateWRMS(TS ts)
{
  TS_GLLE     *gl = (TS_GLLE *)ts->data;
  PetscScalar *x, *w;
  PetscInt     n, i;

  PetscFunctionBegin;
  PetscCall(VecGetArray(gl->X[0], &x));
  PetscCall(VecGetArray(gl->W, &w));
  PetscCall(VecGetLocalSize(gl->W, &n));
  for (i = 0; i < n; i++) w[i] = 1. / (gl->wrms_atol + gl->wrms_rtol * PetscAbsScalar(x[i]));
  PetscCall(VecRestoreArray(gl->X[0], &x));
  PetscCall(VecRestoreArray(gl->W, &w));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLEVecNormWRMS(TS ts, Vec X, PetscReal *nrm)
{
  TS_GLLE     *gl = (TS_GLLE *)ts->data;
  PetscScalar *x, *w;
  PetscReal    sum = 0.0, gsum;
  PetscInt     n, N, i;

  PetscFunctionBegin;
  PetscCall(VecGetArray(X, &x));
  PetscCall(VecGetArray(gl->W, &w));
  PetscCall(VecGetLocalSize(gl->W, &n));
  for (i = 0; i < n; i++) sum += PetscAbsScalar(PetscSqr(x[i] * w[i]));
  PetscCall(VecRestoreArray(X, &x));
  PetscCall(VecRestoreArray(gl->W, &w));
  PetscCallMPI(MPIU_Allreduce(&sum, &gsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
  PetscCall(VecGetSize(gl->W, &N));
  *nrm = PetscSqrtReal(gsum / (1. * N));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLESetType_GLLE(TS ts, TSGLLEType type)
{
  PetscBool same;
  TS_GLLE  *gl = (TS_GLLE *)ts->data;
  PetscErrorCode (*r)(TS);

  PetscFunctionBegin;
  if (gl->type_name[0]) {
    PetscCall(PetscStrcmp(gl->type_name, type, &same));
    if (same) PetscFunctionReturn(PETSC_SUCCESS);
    PetscCall((*gl->Destroy)(gl));
  }

  PetscCall(PetscFunctionListFind(TSGLLEList, type, &r));
  PetscCheck(r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLE type \"%s\" given", type);
  PetscCall((*r)(ts));
  PetscCall(PetscStrncpy(gl->type_name, type, sizeof(gl->type_name)));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLESetAcceptType_GLLE(TS ts, TSGLLEAcceptType type)
{
  TSGLLEAcceptFn *r;
  TS_GLLE        *gl = (TS_GLLE *)ts->data;

  PetscFunctionBegin;
  PetscCall(PetscFunctionListFind(TSGLLEAcceptList, type, &r));
  PetscCheck(r, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown TSGLLEAccept type \"%s\" given", type);
  gl->Accept = r;
  PetscCall(PetscStrncpy(gl->accept_name, type, sizeof(gl->accept_name)));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLEGetAdapt_GLLE(TS ts, TSGLLEAdapt *adapt)
{
  TS_GLLE *gl = (TS_GLLE *)ts->data;

  PetscFunctionBegin;
  if (!gl->adapt) {
    PetscCall(TSGLLEAdaptCreate(PetscObjectComm((PetscObject)ts), &gl->adapt));
    PetscCall(PetscObjectIncrementTabLevel((PetscObject)gl->adapt, (PetscObject)ts, 1));
  }
  *adapt = gl->adapt;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLEChooseNextScheme(TS ts, PetscReal h, const PetscReal hmnorm[], PetscInt *next_scheme, PetscReal *next_h, PetscBool *finish)
{
  TS_GLLE  *gl = (TS_GLLE *)ts->data;
  PetscInt  i, n, cur_p, cur, next_sc, candidates[64], orders[64];
  PetscReal errors[64], costs[64], tleft;

  PetscFunctionBegin;
  cur   = -1;
  cur_p = gl->schemes[gl->current_scheme]->p;
  tleft = ts->max_time - (ts->ptime + ts->time_step);
  for (i = 0, n = 0; i < gl->nschemes; i++) {
    TSGLLEScheme sc = gl->schemes[i];
    if (sc->p < gl->min_order || gl->max_order < sc->p) continue;
    if (sc->p == cur_p - 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[0];
    else if (sc->p == cur_p) errors[n] = PetscAbsScalar(sc->alpha[0]) * hmnorm[1];
    else if (sc->p == cur_p + 1) errors[n] = PetscAbsScalar(sc->alpha[0]) * (hmnorm[2] + hmnorm[3]);
    else continue;
    candidates[n] = i;
    orders[n]     = PetscMin(sc->p, sc->q); /* order of global truncation error */
    costs[n]      = sc->s;                  /* estimate the cost as the number of stages */
    if (i == gl->current_scheme) cur = n;
    n++;
  }
  PetscCheck(cur >= 0 && gl->nschemes > cur, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Current scheme not found in scheme list");
  PetscCall(TSGLLEAdaptChoose(gl->adapt, n, orders, errors, costs, cur, h, tleft, &next_sc, next_h, finish));
  *next_scheme = candidates[next_sc];
  PetscCall(PetscInfo(ts, "Adapt chose scheme %" PetscInt_FMT " (%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ") with step size %6.2e, finish=%s\n", *next_scheme, gl->schemes[*next_scheme]->p, gl->schemes[*next_scheme]->q,
                      gl->schemes[*next_scheme]->r, gl->schemes[*next_scheme]->s, (double)*next_h, PetscBools[*finish]));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSGLLEGetMaxSizes(TS ts, PetscInt *max_r, PetscInt *max_s)
{
  TS_GLLE *gl = (TS_GLLE *)ts->data;

  PetscFunctionBegin;
  *max_r = gl->schemes[gl->nschemes - 1]->r;
  *max_s = gl->schemes[gl->nschemes - 1]->s;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSSolve_GLLE(TS ts)
{
  TS_GLLE            *gl = (TS_GLLE *)ts->data;
  PetscInt            i, k, its, lits, max_r, max_s;
  PetscBool           final_step, finish;
  SNESConvergedReason snesreason;

  PetscFunctionBegin;
  PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));

  PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
  PetscCall(VecCopy(ts->vec_sol, gl->X[0]));
  for (i = 1; i < max_r; i++) PetscCall(VecZeroEntries(gl->X[i]));
  PetscCall(TSGLLEUpdateWRMS(ts));

  if (0) {
    /* Find consistent initial data for DAE */
    gl->stage_time = ts->ptime + ts->time_step;
    gl->scoeff     = 1.;
    gl->stage      = 0;

    PetscCall(VecCopy(ts->vec_sol, gl->Z));
    PetscCall(VecCopy(ts->vec_sol, gl->Y));
    PetscCall(SNESSolve(ts->snes, NULL, gl->Y));
    PetscCall(SNESGetIterationNumber(ts->snes, &its));
    PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
    PetscCall(SNESGetConvergedReason(ts->snes, &snesreason));

    ts->snes_its += its;
    ts->ksp_its += lits;
    if (snesreason < 0 && ts->max_snes_failures != PETSC_UNLIMITED && ++ts->num_snes_failures >= ts->max_snes_failures) {
      ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
      PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures));
      PetscFunctionReturn(PETSC_SUCCESS);
    }
  }

  PetscCheck(gl->current_scheme >= 0, PETSC_COMM_SELF, PETSC_ERR_ORDER, "A starting scheme has not been provided");

  for (k = 0, final_step = PETSC_FALSE, finish = PETSC_FALSE; k < ts->max_steps && !finish; k++) {
    PetscInt           j, r, s, next_scheme = 0, rejections;
    PetscReal          h, hmnorm[4], enorm[3], next_h;
    PetscBool          accept;
    const PetscScalar *c, *a, *u;
    Vec               *X, *Ydot, Y;
    TSGLLEScheme       scheme = gl->schemes[gl->current_scheme];

    r    = scheme->r;
    s    = scheme->s;
    c    = scheme->c;
    a    = scheme->a;
    u    = scheme->u;
    h    = ts->time_step;
    X    = gl->X;
    Ydot = gl->Ydot;
    Y    = gl->Y;

    if (ts->ptime > ts->max_time) break;

    /*
      We only call PreStep at the start of each STEP, not each STAGE.  This is because it is
      possible to fail (have to restart a step) after multiple stages.
    */
    PetscCall(TSPreStep(ts));

    rejections = 0;
    while (1) {
      for (i = 0; i < s; i++) {
        PetscScalar shift;
        gl->scoeff     = 1. / PetscRealPart(a[i * s + i]);
        shift          = gl->scoeff / ts->time_step;
        gl->stage      = i;
        gl->stage_time = ts->ptime + PetscRealPart(c[i]) * h;

        /*
        * Stage equation: Y = h A Y' + U X
        * We assume that A is lower-triangular so that we can solve the stages (Y,Y') sequentially
        * Build the affine vector z_i = -[1/(h a_ii)](h sum_j a_ij y'_j + sum_j u_ij x_j)
        * Then y'_i = z + 1/(h a_ii) y_i
        */
        PetscCall(VecZeroEntries(gl->Z));
        for (j = 0; j < r; j++) PetscCall(VecAXPY(gl->Z, -shift * u[i * r + j], X[j]));
        for (j = 0; j < i; j++) PetscCall(VecAXPY(gl->Z, -shift * h * a[i * s + j], Ydot[j]));
        /* Note: Z is used within function evaluation, Ydot = Z + shift*Y */

        /* Compute an estimate of Y to start Newton iteration */
        if (gl->extrapolate) {
          if (i == 0) {
            /* Linear extrapolation on the first stage */
            PetscCall(VecWAXPY(Y, c[i] * h, X[1], X[0]));
          } else {
            /* Linear extrapolation from the last stage */
            PetscCall(VecAXPY(Y, (c[i] - c[i - 1]) * h, Ydot[i - 1]));
          }
        } else if (i == 0) { /* Directly use solution from the last step, otherwise reuse the last stage (do nothing) */
          PetscCall(VecCopy(X[0], Y));
        }

        /* Solve this stage (Ydot[i] is computed during function evaluation) */
        PetscCall(SNESSolve(ts->snes, NULL, Y));
        PetscCall(SNESGetIterationNumber(ts->snes, &its));
        PetscCall(SNESGetLinearSolveIterations(ts->snes, &lits));
        PetscCall(SNESGetConvergedReason(ts->snes, &snesreason));
        ts->snes_its += its;
        ts->ksp_its += lits;
        if (snesreason < 0 && ts->max_snes_failures != PETSC_UNLIMITED && ++ts->num_snes_failures >= ts->max_snes_failures) {
          ts->reason = TS_DIVERGED_NONLINEAR_SOLVE;
          PetscCall(PetscInfo(ts, "Step=%" PetscInt_FMT ", nonlinear solve failures %" PetscInt_FMT " greater than current TS allowed, stopping solve\n", ts->steps, ts->num_snes_failures));
          PetscFunctionReturn(PETSC_SUCCESS);
        }
      }

      gl->stage_time = ts->ptime + ts->time_step;

      PetscCall((*gl->EstimateHigherMoments)(scheme, h, Ydot, gl->X, gl->himom));
      /* hmnorm[i] = h^{p+i}x^{(p+i)} with i=0,1,2; hmnorm[3] = h^{p+2}(dx'/dx) x^{(p+1)} */
      for (i = 0; i < 3; i++) PetscCall(TSGLLEVecNormWRMS(ts, gl->himom[i], &hmnorm[i + 1]));
      enorm[0] = PetscRealPart(scheme->alpha[0]) * hmnorm[1];
      enorm[1] = PetscRealPart(scheme->beta[0]) * hmnorm[2];
      enorm[2] = PetscRealPart(scheme->gamma[0]) * hmnorm[3];
      PetscCall((*gl->Accept)(ts, ts->max_time - gl->stage_time, h, enorm, &accept));
      if (accept) goto accepted;
      rejections++;
      PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " (t=%g) not accepted, rejections=%" PetscInt_FMT "\n", k, (double)gl->stage_time, rejections));
      if (rejections > gl->max_step_rejections) break;
      /*
        There are lots of reasons why a step might be rejected, including solvers not converging and other factors that
        TSGLLEChooseNextScheme does not support.  Additionally, the error estimates may be very screwed up, so I'm not
        convinced that it's safe to just compute a new error estimate using the same interface as the current adaptor
        (the adaptor interface probably has to change).  Here we make an arbitrary and naive choice.  This assumes that
        steps were written in Nordsieck form.  The "correct" method would be to re-complete the previous time step with
        the correct "next" step size.  It is unclear to me whether the present ad-hoc method of rescaling X is stable.
      */
      h *= 0.5;
      for (i = 1; i < scheme->r; i++) PetscCall(VecScale(X[i], PetscPowRealInt(0.5, i)));
    }
    SETERRQ(PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Time step %" PetscInt_FMT " (t=%g) not accepted after %" PetscInt_FMT " failures", k, (double)gl->stage_time, rejections);

  accepted:
    /* This term is not error, but it *would* be the leading term for a lower order method */
    PetscCall(TSGLLEVecNormWRMS(ts, gl->X[scheme->r - 1], &hmnorm[0]));
    /* Correct scaling so that these are equivalent to norms of the Nordsieck vectors */

    PetscCall(PetscInfo(ts, "Last moment norm %10.2e, estimated error norms %10.2e %10.2e %10.2e\n", (double)hmnorm[0], (double)enorm[0], (double)enorm[1], (double)enorm[2]));
    if (!final_step) {
      PetscCall(TSGLLEChooseNextScheme(ts, h, hmnorm, &next_scheme, &next_h, &final_step));
    } else {
      /* Dummy values to complete the current step in a consistent manner */
      next_scheme = gl->current_scheme;
      next_h      = h;
      finish      = PETSC_TRUE;
    }

    X        = gl->Xold;
    gl->Xold = gl->X;
    gl->X    = X;
    PetscCall((*gl->CompleteStep)(scheme, h, gl->schemes[next_scheme], next_h, Ydot, gl->Xold, gl->X));

    PetscCall(TSGLLEUpdateWRMS(ts));

    /* Post the solution for the user, we could avoid this copy with a small bit of cleverness */
    PetscCall(VecCopy(gl->X[0], ts->vec_sol));
    ts->ptime += h;
    ts->steps++;

    PetscCall(TSPostEvaluate(ts));
    PetscCall(TSPostStep(ts));
    PetscCall(TSMonitor(ts, ts->steps, ts->ptime, ts->vec_sol));

    gl->current_scheme = next_scheme;
    ts->time_step      = next_h;
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSReset_GLLE(TS ts)
{
  TS_GLLE *gl = (TS_GLLE *)ts->data;
  PetscInt max_r, max_s;

  PetscFunctionBegin;
  if (gl->setupcalled) {
    PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
    PetscCall(VecDestroyVecs(max_r, &gl->Xold));
    PetscCall(VecDestroyVecs(max_r, &gl->X));
    PetscCall(VecDestroyVecs(max_s, &gl->Ydot));
    PetscCall(VecDestroyVecs(3, &gl->himom));
    PetscCall(VecDestroy(&gl->W));
    PetscCall(VecDestroy(&gl->Y));
    PetscCall(VecDestroy(&gl->Z));
  }
  gl->setupcalled = PETSC_FALSE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSDestroy_GLLE(TS ts)
{
  TS_GLLE *gl = (TS_GLLE *)ts->data;

  PetscFunctionBegin;
  PetscCall(TSReset_GLLE(ts));
  if (ts->dm) {
    PetscCall(DMCoarsenHookRemove(ts->dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts));
    PetscCall(DMSubDomainHookRemove(ts->dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts));
  }
  if (gl->adapt) PetscCall(TSGLLEAdaptDestroy(&gl->adapt));
  if (gl->Destroy) PetscCall((*gl->Destroy)(gl));
  PetscCall(PetscFree(ts->data));
  PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", NULL));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
    This defines the nonlinear equation that is to be solved with SNES
    g(x) = f(t,x,z+shift*x) = 0
*/
static PetscErrorCode SNESTSFormFunction_GLLE(SNES snes, Vec x, Vec f, TS ts)
{
  TS_GLLE *gl = (TS_GLLE *)ts->data;
  Vec      Z, Ydot;
  DM       dm, dmsave;

  PetscFunctionBegin;
  PetscCall(SNESGetDM(snes, &dm));
  PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot));
  PetscCall(VecWAXPY(Ydot, gl->scoeff / ts->time_step, x, Z));
  dmsave = ts->dm;
  ts->dm = dm;
  PetscCall(TSComputeIFunction(ts, gl->stage_time, x, Ydot, f, PETSC_FALSE));
  ts->dm = dmsave;
  PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode SNESTSFormJacobian_GLLE(SNES snes, Vec x, Mat A, Mat B, TS ts)
{
  TS_GLLE *gl = (TS_GLLE *)ts->data;
  Vec      Z, Ydot;
  DM       dm, dmsave;

  PetscFunctionBegin;
  PetscCall(SNESGetDM(snes, &dm));
  PetscCall(TSGLLEGetVecs(ts, dm, &Z, &Ydot));
  dmsave = ts->dm;
  ts->dm = dm;
  /* gl->Xdot will have already been computed in SNESTSFormFunction_GLLE */
  PetscCall(TSComputeIJacobian(ts, gl->stage_time, x, gl->Ydot[gl->stage], gl->scoeff / ts->time_step, A, B, PETSC_FALSE));
  ts->dm = dmsave;
  PetscCall(TSGLLERestoreVecs(ts, dm, &Z, &Ydot));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSSetUp_GLLE(TS ts)
{
  TS_GLLE *gl = (TS_GLLE *)ts->data;
  PetscInt max_r, max_s;
  DM       dm;

  PetscFunctionBegin;
  gl->setupcalled = PETSC_TRUE;
  PetscCall(TSGLLEGetMaxSizes(ts, &max_r, &max_s));
  PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->X));
  PetscCall(VecDuplicateVecs(ts->vec_sol, max_r, &gl->Xold));
  PetscCall(VecDuplicateVecs(ts->vec_sol, max_s, &gl->Ydot));
  PetscCall(VecDuplicateVecs(ts->vec_sol, 3, &gl->himom));
  PetscCall(VecDuplicate(ts->vec_sol, &gl->W));
  PetscCall(VecDuplicate(ts->vec_sol, &gl->Y));
  PetscCall(VecDuplicate(ts->vec_sol, &gl->Z));

  /* Default acceptance tests and adaptivity */
  if (!gl->Accept) PetscCall(TSGLLESetAcceptType(ts, TSGLLEACCEPT_ALWAYS));
  if (!gl->adapt) PetscCall(TSGLLEGetAdapt(ts, &gl->adapt));

  if (gl->current_scheme < 0) {
    PetscInt i;
    for (i = 0;; i++) {
      if (gl->schemes[i]->p == gl->start_order) break;
      PetscCheck(i + 1 != gl->nschemes, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No schemes available with requested start order %" PetscInt_FMT, i);
    }
    gl->current_scheme = i;
  }
  PetscCall(TSGetDM(ts, &dm));
  PetscCall(DMCoarsenHookAdd(dm, DMCoarsenHook_TSGLLE, DMRestrictHook_TSGLLE, ts));
  PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_TSGLLE, DMSubDomainRestrictHook_TSGLLE, ts));
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSSetFromOptions_GLLE(TS ts, PetscOptionItems PetscOptionsObject)
{
  TS_GLLE *gl         = (TS_GLLE *)ts->data;
  char     tname[256] = TSGLLE_IRKS, completef[256] = "rescale-and-modify";

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "General Linear ODE solver options");
  {
    PetscBool flg;
    PetscCall(PetscOptionsFList("-ts_gl_type", "Type of GL method", "TSGLLESetType", TSGLLEList, gl->type_name[0] ? gl->type_name : tname, tname, sizeof(tname), &flg));
    if (flg || !gl->type_name[0]) PetscCall(TSGLLESetType(ts, tname));
    PetscCall(PetscOptionsInt("-ts_gl_max_step_rejections", "Maximum number of times to attempt a step", "None", gl->max_step_rejections, &gl->max_step_rejections, NULL));
    PetscCall(PetscOptionsInt("-ts_gl_max_order", "Maximum order to try", "TSGLLESetMaxOrder", gl->max_order, &gl->max_order, NULL));
    PetscCall(PetscOptionsInt("-ts_gl_min_order", "Minimum order to try", "TSGLLESetMinOrder", gl->min_order, &gl->min_order, NULL));
    PetscCall(PetscOptionsInt("-ts_gl_start_order", "Initial order to try", "TSGLLESetMinOrder", gl->start_order, &gl->start_order, NULL));
    PetscCall(PetscOptionsEnum("-ts_gl_error_direction", "Which direction to look when estimating error", "TSGLLESetErrorDirection", TSGLLEErrorDirections, (PetscEnum)gl->error_direction, (PetscEnum *)&gl->error_direction, NULL));
    PetscCall(PetscOptionsBool("-ts_gl_extrapolate", "Extrapolate stage solution from previous solution (sometimes unstable)", "TSGLLESetExtrapolate", gl->extrapolate, &gl->extrapolate, NULL));
    PetscCall(PetscOptionsReal("-ts_gl_atol", "Absolute tolerance", "TSGLLESetTolerances", gl->wrms_atol, &gl->wrms_atol, NULL));
    PetscCall(PetscOptionsReal("-ts_gl_rtol", "Relative tolerance", "TSGLLESetTolerances", gl->wrms_rtol, &gl->wrms_rtol, NULL));
    PetscCall(PetscOptionsString("-ts_gl_complete", "Method to use for completing the step", "none", completef, completef, sizeof(completef), &flg));
    if (flg) {
      PetscBool match1, match2;
      PetscCall(PetscStrcmp(completef, "rescale", &match1));
      PetscCall(PetscStrcmp(completef, "rescale-and-modify", &match2));
      if (match1) gl->CompleteStep = TSGLLECompleteStep_Rescale;
      else if (match2) gl->CompleteStep = TSGLLECompleteStep_RescaleAndModify;
      else SETERRQ(PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_UNKNOWN_TYPE, "%s", completef);
    }
    {
      char type[256] = TSGLLEACCEPT_ALWAYS;
      PetscCall(PetscOptionsFList("-ts_gl_accept_type", "Method to use for determining whether to accept a step", "TSGLLESetAcceptType", TSGLLEAcceptList, gl->accept_name[0] ? gl->accept_name : type, type, sizeof(type), &flg));
      if (flg || !gl->accept_name[0]) PetscCall(TSGLLESetAcceptType(ts, type));
    }
    {
      TSGLLEAdapt adapt;
      PetscCall(TSGLLEGetAdapt(ts, &adapt));
      PetscCall(TSGLLEAdaptSetFromOptions(adapt, PetscOptionsObject));
    }
  }
  PetscOptionsHeadEnd();
  PetscFunctionReturn(PETSC_SUCCESS);
}

static PetscErrorCode TSView_GLLE(TS ts, PetscViewer viewer)
{
  TS_GLLE  *gl = (TS_GLLE *)ts->data;
  PetscInt  i;
  PetscBool isascii, details;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
  if (isascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "  min order %" PetscInt_FMT ", max order %" PetscInt_FMT ", current order %" PetscInt_FMT "\n", gl->min_order, gl->max_order, gl->schemes[gl->current_scheme]->p));
    PetscCall(PetscViewerASCIIPrintf(viewer, "  Error estimation: %s\n", TSGLLEErrorDirections[gl->error_direction]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "  Extrapolation: %s\n", gl->extrapolate ? "yes" : "no"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "  Acceptance test: %s\n", gl->accept_name[0] ? gl->accept_name : "(not yet set)"));
    PetscCall(PetscViewerASCIIPushTab(viewer));
    PetscCall(TSGLLEAdaptView(gl->adapt, viewer));
    PetscCall(PetscViewerASCIIPopTab(viewer));
    PetscCall(PetscViewerASCIIPrintf(viewer, "  type: %s\n", gl->type_name[0] ? gl->type_name : "(not yet set)"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "Schemes within family (%" PetscInt_FMT "):\n", gl->nschemes));
    details = PETSC_FALSE;
    PetscCall(PetscOptionsGetBool(((PetscObject)ts)->options, ((PetscObject)ts)->prefix, "-ts_gl_view_detailed", &details, NULL));
    PetscCall(PetscViewerASCIIPushTab(viewer));
    for (i = 0; i < gl->nschemes; i++) PetscCall(TSGLLESchemeView(gl->schemes[i], details, viewer));
    if (gl->View) PetscCall((*gl->View)(gl, viewer));
    PetscCall(PetscViewerASCIIPopTab(viewer));
  }
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  TSGLLERegister -  adds a `TSGLLE` implementation

  Not Collective, No Fortran Support

  Input Parameters:
+ sname    - name of user-defined general linear scheme
- function - routine to create method context

  Level: advanced

  Note:
  `TSGLLERegister()` may be called multiple times to add several user-defined families.

  Example Usage:
.vb
  TSGLLERegister("my_scheme", MySchemeCreate);
.ve

  Then, your scheme can be chosen with the procedural interface via
.vb
  TSGLLESetType(ts, "my_scheme")
.ve
  or at runtime via the option
.vb
  -ts_gl_type my_scheme
.ve

.seealso: [](ch_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`
@*/
PetscErrorCode TSGLLERegister(const char sname[], PetscErrorCode (*function)(TS))
{
  PetscFunctionBegin;
  PetscCall(TSGLLEInitializePackage());
  PetscCall(PetscFunctionListAdd(&TSGLLEList, sname, function));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  TSGLLEAcceptRegister -  adds a `TSGLLE` acceptance scheme

  Not Collective

  Input Parameters:
+ sname    - name of user-defined acceptance scheme
- function - routine to create method context, see `TSGLLEAcceptFn` for the calling sequence

  Level: advanced

  Note:
  `TSGLLEAcceptRegister()` may be called multiple times to add several user-defined families.

  Example Usage:
.vb
  TSGLLEAcceptRegister("my_scheme", MySchemeCreate);
.ve

  Then, your scheme can be chosen with the procedural interface via
.vb
  TSGLLESetAcceptType(ts, "my_scheme")
.ve
  or at runtime via the option `-ts_gl_accept_type my_scheme`

.seealso: [](ch_ts), `TSGLLE`, `TSGLLEType`, `TSGLLERegisterAll()`, `TSGLLEAcceptFn`
@*/
PetscErrorCode TSGLLEAcceptRegister(const char sname[], TSGLLEAcceptFn *function)
{
  PetscFunctionBegin;
  PetscCall(PetscFunctionListAdd(&TSGLLEAcceptList, sname, function));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  TSGLLERegisterAll - Registers all of the general linear methods in `TSGLLE`

  Not Collective

  Level: advanced

.seealso: [](ch_ts), `TSGLLE`, `TSGLLERegisterDestroy()`
@*/
PetscErrorCode TSGLLERegisterAll(void)
{
  PetscFunctionBegin;
  if (TSGLLERegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
  TSGLLERegisterAllCalled = PETSC_TRUE;

  PetscCall(TSGLLERegister(TSGLLE_IRKS, TSGLLECreate_IRKS));
  PetscCall(TSGLLEAcceptRegister(TSGLLEACCEPT_ALWAYS, TSGLLEAccept_Always));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  TSGLLEInitializePackage - This function initializes everything in the `TSGLLE` package. It is called
  from `TSInitializePackage()`.

  Level: developer

.seealso: [](ch_ts), `PetscInitialize()`, `TSInitializePackage()`, `TSGLLEFinalizePackage()`
@*/
PetscErrorCode TSGLLEInitializePackage(void)
{
  PetscFunctionBegin;
  if (TSGLLEPackageInitialized) PetscFunctionReturn(PETSC_SUCCESS);
  TSGLLEPackageInitialized = PETSC_TRUE;
  PetscCall(TSGLLERegisterAll());
  PetscCall(PetscRegisterFinalize(TSGLLEFinalizePackage));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*@C
  TSGLLEFinalizePackage - This function destroys everything in the `TSGLLE` package. It is
  called from `PetscFinalize()`.

  Level: developer

.seealso: [](ch_ts), `PetscFinalize()`, `TSGLLEInitializePackage()`, `TSInitializePackage()`
@*/
PetscErrorCode TSGLLEFinalizePackage(void)
{
  PetscFunctionBegin;
  PetscCall(PetscFunctionListDestroy(&TSGLLEList));
  PetscCall(PetscFunctionListDestroy(&TSGLLEAcceptList));
  TSGLLEPackageInitialized = PETSC_FALSE;
  TSGLLERegisterAllCalled  = PETSC_FALSE;
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*MC
  TSGLLE - DAE solver using implicit General Linear methods {cite}`butcher_2007` {cite}`butcher2016numerical`

  Options Database Keys:
+  -ts_gl_type <type> - the class of general linear method (irks)
.  -ts_gl_rtol <tol>  - relative error
.  -ts_gl_atol <tol>  - absolute error
.  -ts_gl_min_order <p> - minimum order method to consider (default=1)
.  -ts_gl_max_order <p> - maximum order method to consider (default=3)
.  -ts_gl_start_order <p> - order of starting method (default=1)
.  -ts_gl_complete <method> - method to use for completing the step (rescale-and-modify or rescale)
-  -ts_adapt_type <method> - adaptive controller to use (none step both)

  Level: beginner

  Notes:
  These methods contain Runge-Kutta and multistep schemes as special cases. These special cases
  have some fundamental limitations. For example, diagonally implicit Runge-Kutta cannot have
  stage order greater than 1 which limits their applicability to very stiff systems.
  Meanwhile, multistep methods cannot be A-stable for order greater than 2 and BDF are not
  0-stable for order greater than 6. GL methods can be A- and L-stable with arbitrarily high
  stage order and reliable error estimates for both 1 and 2 orders higher to facilitate
  adaptive step sizes and adaptive order schemes. All this is possible while preserving a
  singly diagonally implicit structure.

  This integrator can be applied to DAE.

  Diagonally implicit general linear (DIGL) methods are a generalization of diagonally implicit
  Runge-Kutta (DIRK). They are represented by the tableau

.vb
  A  |  U
  -------
  B  |  V
.ve

  combined with a vector c of abscissa. "Diagonally implicit" means that $A$ is lower
  triangular. A step of the general method reads

  $$
  \begin{align*}
  [ Y ] = [A  U] [  Y'   ] \\
  [X^k] = [B  V] [X^{k-1}]
  \end{align*}
  $$

  where Y is the multivector of stage values, $Y'$ is the multivector of stage derivatives, $X^k$
  is the Nordsieck vector of the solution at step $k$. The Nordsieck vector consists of the first
  $r$ moments of the solution, given by

  $$
  X = [x_0,x_1,...,x_{r-1}] = [x, h x', h^2 x'', ..., h^{r-1} x^{(r-1)} ]
  $$

  If $A$ is lower triangular, we can solve the stages $(Y, Y')$ sequentially

  $$
  y_i = h \sum_{j=0}^{s-1} (a_{ij} y'_j) + \sum_{j=0}^{r-1} u_{ij} x_j, \, \,   i=0,...,{s-1}
  $$

  and then construct the pieces to carry to the next step

  $$
  xx_i = h \sum_{j=0}^{s-1} b_{ij} y'_j  + \sum_{j=0}^{r-1} v_{ij} x_j,  \, \,  i=0,...,{r-1}
  $$

  Note that when the equations are cast in implicit form, we are using the stage equation to
  define $y'_i$ in terms of $y_i$ and known stuff ($y_j$ for $j<i$ and $x_j$ for all $j$).

  Error estimation

  At present, the most attractive GL methods for stiff problems are singly diagonally implicit
  schemes which posses Inherent Runge-Kutta Stability (`TSIRKS`).  These methods have $r=s$, the
  number of items passed between steps is equal to the number of stages.  The order and
  stage-order are one less than the number of stages.  We use the error estimates in the 2007
  paper which provide the following estimates

  $$
  \begin{align*}
  h^{p+1} X^{(p+1)}          = \phi_0^T Y' + [0 \psi_0^T] Xold \\
  h^{p+2} X^{(p+2)}          = \phi_1^T Y' + [0 \psi_1^T] Xold \\
  h^{p+2} (dx'/dx) X^{(p+1)} = \phi_2^T Y' + [0 \psi_2^T] Xold
  \end{align*}
  $$

  These estimates are accurate to $ O(h^{p+3})$.

  Changing the step size

  Uses the generalized "rescale and modify" scheme, see equation (4.5) of {cite}`butcher_2007`.

.seealso: [](ch_ts), `TSCreate()`, `TS`, `TSSetType()`, `TSType`
M*/
PETSC_EXTERN PetscErrorCode TSCreate_GLLE(TS ts)
{
  TS_GLLE *gl;

  PetscFunctionBegin;
  PetscCall(TSGLLEInitializePackage());

  PetscCall(PetscNew(&gl));
  ts->data = (void *)gl;

  ts->ops->reset          = TSReset_GLLE;
  ts->ops->destroy        = TSDestroy_GLLE;
  ts->ops->view           = TSView_GLLE;
  ts->ops->setup          = TSSetUp_GLLE;
  ts->ops->solve          = TSSolve_GLLE;
  ts->ops->setfromoptions = TSSetFromOptions_GLLE;
  ts->ops->snesfunction   = SNESTSFormFunction_GLLE;
  ts->ops->snesjacobian   = SNESTSFormJacobian_GLLE;

  ts->usessnes = PETSC_TRUE;

  gl->max_step_rejections = 1;
  gl->min_order           = 1;
  gl->max_order           = 3;
  gl->start_order         = 1;
  gl->current_scheme      = -1;
  gl->extrapolate         = PETSC_FALSE;

  gl->wrms_atol = 1e-8;
  gl->wrms_rtol = 1e-5;

  PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetType_C", &TSGLLESetType_GLLE));
  PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLESetAcceptType_C", &TSGLLESetAcceptType_GLLE));
  PetscCall(PetscObjectComposeFunction((PetscObject)ts, "TSGLLEGetAdapt_C", &TSGLLEGetAdapt_GLLE));
  PetscFunctionReturn(PETSC_SUCCESS);
}