diff --git a/bindings/sundials4py/sundomeigest/sundomeigest_power_generated.hpp b/bindings/sundials4py/sundomeigest/sundomeigest_power_generated.hpp
index 27b613dd62..5928fb972e 100644
--- a/bindings/sundials4py/sundomeigest/sundomeigest_power_generated.hpp
+++ b/bindings/sundials4py/sundomeigest/sundomeigest_power_generated.hpp
@@ -33,6 +33,10 @@ m.def(
   nb::arg("q"), nb::arg("max_iters"), nb::arg("rel_tol"), nb::arg("sunctx"),
   "nb::keep_alive<0, 4>()", nb::keep_alive<0, 4>());
 
+m.def("SUNDomEigEstimator_SetRhsAtLinearizationPoint_Power",
+      SUNDomEigEstimator_SetRhsAtLinearizationPoint_Power, nb::arg("DEE"),
+      nb::arg("Fyt"));
+
 m.def("SUNDomEigEstimator_SetIsReal_Power", SUNDomEigEstimator_SetIsReal_Power,
       nb::arg("DEE"), nb::arg("real"));
 // #ifdef __cplusplus
diff --git a/include/sundials/sundials_domeigestimator.h b/include/sundials/sundials_domeigestimator.h
index a7b4be5d6a..e824a20ac8 100644
--- a/include/sundials/sundials_domeigestimator.h
+++ b/include/sundials/sundials_domeigestimator.h
@@ -52,6 +52,7 @@ struct SUNDomEigEstimator_Ops_
   SUNErrCode (*setrhs)(SUNDomEigEstimator, void*, SUNRhsFn);
   SUNErrCode (*setrhslinearizationpoint)(SUNDomEigEstimator, sunrealtype,
                                          N_Vector);
+  SUNErrCode (*setrhsatlinearizationpoint)(SUNDomEigEstimator, N_Vector);
   SUNErrCode (*setoptions)(SUNDomEigEstimator DEE, const char* Did,
                            const char* file_name, int argc, char* argv[]);
   SUNErrCode (*setmaxiters)(SUNDomEigEstimator, long int);
@@ -101,6 +102,10 @@ SUNDIALS_EXPORT
 SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint(SUNDomEigEstimator DEE,
                                                        sunrealtype t, N_Vector v);
 
+SUNDIALS_EXPORT
+SUNErrCode SUNDomEigEstimator_SetRhsAtLinearizationPoint(SUNDomEigEstimator DEE,
+                                                         N_Vector Fyt);
+
 SUNDIALS_EXPORT
 SUNErrCode SUNDomEigEstimator_SetOptions(SUNDomEigEstimator DEE,
                                          const char* Did, const char* file_name,
diff --git a/include/sundomeigest/sundomeigest_arnoldi.h b/include/sundomeigest/sundomeigest_arnoldi.h
index b46dec25b8..5f01c6e1f3 100644
--- a/include/sundomeigest/sundomeigest_arnoldi.h
+++ b/include/sundomeigest/sundomeigest_arnoldi.h
@@ -43,12 +43,13 @@ struct SUNDomEigEstimatorContent_Arnoldi_
 
   /* Krylov subspace vectors */
   N_Vector* V;
-  N_Vector q, rhs_linY, Fy, work;
+  N_Vector rhs_linY, Fy, work;
 
   int kry_dim;                  /* Krylov subspace dimension */
   int num_warmups;              /* Number of preprocessing iterations */
   long int num_iters;           /* Number of iterations in last Estimate call */
   sunbooleantype warmup_to_tol; /* Whether to use warmup iterations */
+  sunbooleantype Fy_is_current; /* Flag to track if Fy is current */
   sunrealtype tol_warmup;       /* Tolerance for warmup iterations */
   sunrealtype rhs_linT;         /* Time value for linearization point */
 
@@ -90,6 +91,10 @@ SUNDIALS_EXPORT
 SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(
   SUNDomEigEstimator DEE, sunrealtype t, N_Vector v);
 
+SUNDIALS_EXPORT
+SUNErrCode SUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(
+  SUNDomEigEstimator DEE, N_Vector Fyt);
+
 SUNDIALS_EXPORT
 SUNErrCode SUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(SUNDomEigEstimator DEE,
                                                             int num_iters);
diff --git a/include/sundomeigest/sundomeigest_power.h b/include/sundomeigest/sundomeigest_power.h
index d9df90215d..cc629c0835 100644
--- a/include/sundomeigest/sundomeigest_power.h
+++ b/include/sundomeigest/sundomeigest_power.h
@@ -41,7 +41,7 @@ struct SUNDomEigEstimatorContent_Power_
   SUNATimesFn ATimes; /* User provided ATimes function */
   void* ATdata;       /* ATimes function data*/
 
-  N_Vector V, q, q_prev, rhs_linY, Fy, work; /* workspace vectors */
+  N_Vector V, Av, v_prev, rhs_linY, Fy, work; /* workspace vectors */
 
   int num_warmups;      /* Number of preprocessing iterations */
   long int max_iters;   /* Maximum number of power iterations */
@@ -57,7 +57,8 @@ struct SUNDomEigEstimatorContent_Power_
   void* rhs_data;   /* RHS function data */
   long int nfevals; /* Number of RHS evaluations */
 
-  sunbooleantype is_complex; /* Flag for complex eigenvalue request */
+  sunbooleantype is_complex;    /* Flag for complex eigenvalue request */
+  sunbooleantype Fy_is_current; /* Flag to track if Fy is current */
 };
 
 typedef struct SUNDomEigEstimatorContent_Power_* SUNDomEigEstimatorContent_Power;
@@ -99,6 +100,10 @@ SUNDIALS_EXPORT
 SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint_Power(
   SUNDomEigEstimator DEE, sunrealtype t, N_Vector v);
 
+SUNDIALS_EXPORT
+SUNErrCode SUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(
+  SUNDomEigEstimator DEE, N_Vector Fyt);
+
 SUNDIALS_EXPORT
 SUNErrCode SUNDomEigEstimator_SetIsReal_Power(SUNDomEigEstimator DEE,
                                               sunbooleantype real);
diff --git a/src/arkode/arkode_lsrkstep.c b/src/arkode/arkode_lsrkstep.c
index d1b4616391..8ffcff3028 100644
--- a/src/arkode/arkode_lsrkstep.c
+++ b/src/arkode/arkode_lsrkstep.c
@@ -2890,6 +2890,29 @@ int lsrkStep_ComputeNewDomEig(ARKodeMem ark_mem, ARKodeLSRKStepMem step_mem)
 
   if (step_mem->DEE != NULL)
   {
+    retval = SUNDomEigEstimator_SetRhsLinearizationPoint(step_mem->DEE,
+                                                         ark_mem->tn,
+                                                         ark_mem->yn);
+    if (retval != SUN_SUCCESS)
+    {
+      arkProcessError(ark_mem, ARK_DEE_FAIL, __LINE__, __func__, __FILE__,
+                      "SUNDomEigEstimator_SetRhsLinearizationPoint failed");
+      return ARK_DEE_FAIL;
+    }
+
+    /* Pass the current RHS vector fn to the dominant eigenvalue estimator */
+    if (ark_mem->fn_is_current == SUNTRUE)
+    {
+      retval = SUNDomEigEstimator_SetRhsAtLinearizationPoint(step_mem->DEE,
+                                                             ark_mem->fn);
+      if (retval != SUN_SUCCESS)
+      {
+        arkProcessError(ark_mem, ARK_DEE_FAIL, __LINE__, __func__, __FILE__,
+                        "SUNDomEigEstimator_SetRhsAtLinearizationPoint failed");
+        return ARK_DEE_FAIL;
+      }
+    }
+
     retval = SUNDomEigEstimator_Estimate(step_mem->DEE, &step_mem->lambdaR,
                                          &step_mem->lambdaI);
     step_mem->dom_eig_num_evals++;
@@ -2910,6 +2933,14 @@ int lsrkStep_ComputeNewDomEig(ARKodeMem ark_mem, ARKodeLSRKStepMem step_mem)
     }
     step_mem->num_dee_iters += num_iters;
 
+    retval = SUNDomEigEstimator_GetNumRhsEvals(step_mem->DEE, &(step_mem->nfeDQ));
+    if (retval != SUN_SUCCESS)
+    {
+      arkProcessError(ark_mem, ARK_DEE_FAIL, __LINE__, __func__, __FILE__,
+                      "SUNDomEigEstimator_GetNumRhsEvals failed");
+      return ARK_DEE_FAIL;
+    }
+
     /* After the first call to SUNDomEigEstimator_Estimate, the number of warmups is set to
        num_warmups, this allows the successive calls to
        SUNDomEigEstimator_Estimate to use a different number of warmups. */
@@ -3230,89 +3261,6 @@ int lsrkStep_legendre_P_complex(int s, sunrealtype zR, sunrealtype zI,
   return ARK_SUCCESS;
 }
 
-/*---------------------------------------------------------------
-  lsrkStep_DQJtimes:
-
-  This routine generates a difference quotient approximation to
-  the Jacobian-vector product f_y(t,y) * v. The approximation is
-  Jv = [f(y + v*sig) - f(y)]/sig, where sig = 1 / ||v||_WRMS,
-  i.e. the WRMS norm of v*sig is 1.
-  ---------------------------------------------------------------*/
-int lsrkStep_DQJtimes(void* arkode_mem, N_Vector v, N_Vector Jv)
-{
-  ARKodeMem ark_mem;
-  ARKodeLSRKStepMem step_mem;
-
-  sunrealtype sig, siginv;
-  int iter, retval;
-
-  /* access ARKodeLSRKStepMem structure */
-  retval = lsrkStep_AccessARKODEStepMem(arkode_mem, __func__, &ark_mem,
-                                        &step_mem);
-  if (retval != ARK_SUCCESS) { return retval; }
-
-  sunrealtype t = ark_mem->tn;
-  N_Vector y    = ark_mem->yn;
-  N_Vector work = ark_mem->tempv3;
-
-  /* Compute RHS function, if necessary. */
-  if ((!ark_mem->fn_is_current && ark_mem->initsetup) ||
-      (step_mem->step_nst != ark_mem->nst))
-  {
-    /* call the user-supplied pre-RHS function (if supplied) */
-    if (ark_mem->PreRhsFn)
-    {
-      retval = ark_mem->PreRhsFn(t, y, ark_mem->user_data);
-      if (retval != 0) { return ARK_PRERHSFN_FAIL; }
-    }
-
-    retval = step_mem->fe(t, y, ark_mem->fn, ark_mem->user_data);
-    step_mem->nfeDQ++;
-    if (retval != ARK_SUCCESS)
-    {
-      SUNLogExtraDebugVec(ARK_LOGGER, "DomEig JvTimes RHS", ark_mem->fn,
-                          "F_n(:) =");
-      SUNLogInfo(ARK_LOGGER, "DomEig JvTimes",
-                 "status = failed rhs eval, retval = %i", retval);
-      return (ARK_RHSFUNC_FAIL);
-    }
-    ark_mem->fn_is_current = SUNTRUE;
-  }
-
-  /* Initialize perturbation to 1/||v|| */
-  sig = ONE / N_VWrmsNorm(v, ark_mem->ewt);
-
-  for (iter = 0; iter < MAX_DQITERS; iter++)
-  {
-    /* Set work = y + sig*v */
-    N_VLinearSum(sig, v, ONE, y, work);
-
-    /* call the user-supplied pre-RHS function (if supplied) */
-    if (ark_mem->PreRhsFn)
-    {
-      retval = ark_mem->PreRhsFn(t, work, ark_mem->user_data);
-      if (retval != 0) { return ARK_PRERHSFN_FAIL; }
-    }
-    /* Set Jv = f(tn, y+sig*v) */
-    retval = step_mem->fe(t, work, Jv, ark_mem->user_data);
-    step_mem->nfeDQ++;
-    if (retval == 0) { break; }
-    if (retval < 0) { return (-1); }
-
-    /* If f failed recoverably, shrink sig and retry */
-    sig *= SUN_RCONST(0.25);
-  }
-
-  /* If retval still isn't 0, return with a recoverable failure */
-  if (retval > 0) { return (+1); }
-
-  /* Replace Jv by (Jv - fn)/sig */
-  siginv = ONE / sig;
-  N_VLinearSum(siginv, Jv, -siginv, ark_mem->fn, Jv);
-
-  return ARK_SUCCESS;
-}
-
 /*===============================================================
   EOF
   ===============================================================*/
diff --git a/src/arkode/arkode_lsrkstep_impl.h b/src/arkode/arkode_lsrkstep_impl.h
index 47e002e721..088fb54714 100644
--- a/src/arkode/arkode_lsrkstep_impl.h
+++ b/src/arkode/arkode_lsrkstep_impl.h
@@ -254,7 +254,6 @@ int lsrkStep_cheb_T_complex(int s, sunrealtype zR, sunrealtype zI,
                             sunrealtype* TsR, sunrealtype* TsI);
 int lsrkStep_legendre_P_complex(int s, sunrealtype zR, sunrealtype zI,
                                 sunrealtype* PsR, sunrealtype* PsI);
-int lsrkStep_DQJtimes(void* arkode_mem, N_Vector v, N_Vector Jv);
 
 /*===============================================================
   Reusable LSRKStep Error Messages
diff --git a/src/arkode/arkode_lsrkstep_io.c b/src/arkode/arkode_lsrkstep_io.c
index eb1d6a1820..5f72be15ab 100644
--- a/src/arkode/arkode_lsrkstep_io.c
+++ b/src/arkode/arkode_lsrkstep_io.c
@@ -553,12 +553,12 @@ int LSRKStepSetDomEigEstimator(void* arkode_mem, SUNDomEigEstimator DEE)
   /* Attach the DEE pointer to the step memory */
   step_mem->DEE = DEE;
 
-  /* Set the ATimes function for the DEE with A_data = arkode_mem */
-  retval = SUNDomEigEstimator_SetATimes(DEE, arkode_mem, lsrkStep_DQJtimes);
+  /* Set the RHS function for the DEE with user_data */
+  retval = SUNDomEigEstimator_SetRhs(DEE, ark_mem->user_data, step_mem->fe);
   if (retval != SUN_SUCCESS)
   {
     arkProcessError(ark_mem, ARK_DEE_FAIL, __LINE__, __func__, __FILE__,
-                    "SUNDomEigEstimator_SetATimes failed");
+                    "SUNDomEigEstimator_SetRhs failed");
     return ARK_DEE_FAIL;
   }
 
diff --git a/src/sundials/fmod_int32/fsundials_core_mod.c b/src/sundials/fmod_int32/fsundials_core_mod.c
index 189a5533fb..f3af918f59 100644
--- a/src/sundials/fmod_int32/fsundials_core_mod.c
+++ b/src/sundials/fmod_int32/fsundials_core_mod.c
@@ -4026,6 +4026,20 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint(SUNDomEigEstim
 }
 
 
+SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint(SUNDomEigEstimator farg1, N_Vector farg2) {
+  int fresult ;
+  SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
+  N_Vector arg2 = (N_Vector) 0 ;
+  SUNErrCode result;
+  
+  arg1 = (SUNDomEigEstimator)(farg1);
+  arg2 = (N_Vector)(farg2);
+  result = (SUNErrCode)SUNDomEigEstimator_SetRhsAtLinearizationPoint(arg1,arg2);
+  fresult = (SUNErrCode)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetMaxIters(SUNDomEigEstimator farg1, long const *farg2) {
   int fresult ;
   SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
diff --git a/src/sundials/fmod_int32/fsundials_core_mod.f90 b/src/sundials/fmod_int32/fsundials_core_mod.f90
index e8fe90677d..750402806a 100644
--- a/src/sundials/fmod_int32/fsundials_core_mod.f90
+++ b/src/sundials/fmod_int32/fsundials_core_mod.f90
@@ -707,6 +707,7 @@ module fsundials_core_mod
   type(C_FUNPTR), public :: setatimes
   type(C_FUNPTR), public :: setrhs
   type(C_FUNPTR), public :: setrhslinearizationpoint
+  type(C_FUNPTR), public :: setrhsatlinearizationpoint
   type(C_FUNPTR), public :: setoptions
   type(C_FUNPTR), public :: setmaxiters
   type(C_FUNPTR), public :: setnumpreprocessiters
@@ -733,6 +734,7 @@ module fsundials_core_mod
  public :: FSUNDomEigEstimator_SetATimes
  public :: FSUNDomEigEstimator_SetRhs
  public :: FSUNDomEigEstimator_SetRhsLinearizationPoint
+ public :: FSUNDomEigEstimator_SetRhsAtLinearizationPoint
  public :: FSUNDomEigEstimator_SetMaxIters
  public :: FSUNDomEigEstimator_SetNumPreprocessIters
  public :: FSUNDomEigEstimator_SetRelTol
@@ -3103,6 +3105,15 @@ function swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint(farg1, farg2, farg3)
 integer(C_INT) :: fresult
 end function
 
+function swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint(farg1, farg2) &
+bind(C, name="_wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+type(C_PTR), value :: farg1
+type(C_PTR), value :: farg2
+integer(C_INT) :: fresult
+end function
+
 function swigc_FSUNDomEigEstimator_SetMaxIters(farg1, farg2) &
 bind(C, name="_wrap_FSUNDomEigEstimator_SetMaxIters") &
 result(fresult)
@@ -7566,6 +7577,22 @@ function FSUNDomEigEstimator_SetRhsLinearizationPoint(dee, t, v) &
 swig_result = fresult
 end function
 
+function FSUNDomEigEstimator_SetRhsAtLinearizationPoint(dee, fyt) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT) :: swig_result
+type(SUNDomEigEstimator), target, intent(inout) :: dee
+type(N_Vector), target, intent(inout) :: fyt
+integer(C_INT) :: fresult 
+type(C_PTR) :: farg1 
+type(C_PTR) :: farg2 
+
+farg1 = c_loc(dee)
+farg2 = c_loc(fyt)
+fresult = swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint(farg1, farg2)
+swig_result = fresult
+end function
+
 function FSUNDomEigEstimator_SetMaxIters(dee, max_iters) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
diff --git a/src/sundials/fmod_int64/fsundials_core_mod.c b/src/sundials/fmod_int64/fsundials_core_mod.c
index a3981c5914..ad280e89f6 100644
--- a/src/sundials/fmod_int64/fsundials_core_mod.c
+++ b/src/sundials/fmod_int64/fsundials_core_mod.c
@@ -4026,6 +4026,20 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint(SUNDomEigEstim
 }
 
 
+SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint(SUNDomEigEstimator farg1, N_Vector farg2) {
+  int fresult ;
+  SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
+  N_Vector arg2 = (N_Vector) 0 ;
+  SUNErrCode result;
+  
+  arg1 = (SUNDomEigEstimator)(farg1);
+  arg2 = (N_Vector)(farg2);
+  result = (SUNErrCode)SUNDomEigEstimator_SetRhsAtLinearizationPoint(arg1,arg2);
+  fresult = (SUNErrCode)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetMaxIters(SUNDomEigEstimator farg1, long const *farg2) {
   int fresult ;
   SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
diff --git a/src/sundials/fmod_int64/fsundials_core_mod.f90 b/src/sundials/fmod_int64/fsundials_core_mod.f90
index bac979e5b2..3c0ee80e1a 100644
--- a/src/sundials/fmod_int64/fsundials_core_mod.f90
+++ b/src/sundials/fmod_int64/fsundials_core_mod.f90
@@ -707,6 +707,7 @@ module fsundials_core_mod
   type(C_FUNPTR), public :: setatimes
   type(C_FUNPTR), public :: setrhs
   type(C_FUNPTR), public :: setrhslinearizationpoint
+  type(C_FUNPTR), public :: setrhsatlinearizationpoint
   type(C_FUNPTR), public :: setoptions
   type(C_FUNPTR), public :: setmaxiters
   type(C_FUNPTR), public :: setnumpreprocessiters
@@ -733,6 +734,7 @@ module fsundials_core_mod
  public :: FSUNDomEigEstimator_SetATimes
  public :: FSUNDomEigEstimator_SetRhs
  public :: FSUNDomEigEstimator_SetRhsLinearizationPoint
+ public :: FSUNDomEigEstimator_SetRhsAtLinearizationPoint
  public :: FSUNDomEigEstimator_SetMaxIters
  public :: FSUNDomEigEstimator_SetNumPreprocessIters
  public :: FSUNDomEigEstimator_SetRelTol
@@ -3103,6 +3105,15 @@ function swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint(farg1, farg2, farg3)
 integer(C_INT) :: fresult
 end function
 
+function swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint(farg1, farg2) &
+bind(C, name="_wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+type(C_PTR), value :: farg1
+type(C_PTR), value :: farg2
+integer(C_INT) :: fresult
+end function
+
 function swigc_FSUNDomEigEstimator_SetMaxIters(farg1, farg2) &
 bind(C, name="_wrap_FSUNDomEigEstimator_SetMaxIters") &
 result(fresult)
@@ -7566,6 +7577,22 @@ function FSUNDomEigEstimator_SetRhsLinearizationPoint(dee, t, v) &
 swig_result = fresult
 end function
 
+function FSUNDomEigEstimator_SetRhsAtLinearizationPoint(dee, fyt) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT) :: swig_result
+type(SUNDomEigEstimator), target, intent(inout) :: dee
+type(N_Vector), target, intent(inout) :: fyt
+integer(C_INT) :: fresult 
+type(C_PTR) :: farg1 
+type(C_PTR) :: farg2 
+
+farg1 = c_loc(dee)
+farg2 = c_loc(fyt)
+fresult = swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint(farg1, farg2)
+swig_result = fresult
+end function
+
 function FSUNDomEigEstimator_SetMaxIters(dee, max_iters) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
diff --git a/src/sundials/sundials_domeigestimator.c b/src/sundials/sundials_domeigestimator.c
index e8e92851c7..fbf5c977f8 100644
--- a/src/sundials/sundials_domeigestimator.c
+++ b/src/sundials/sundials_domeigestimator.c
@@ -223,6 +223,20 @@ SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint(SUNDomEigEstimator DEE,
   return (ier);
 }
 
+SUNErrCode SUNDomEigEstimator_SetRhsAtLinearizationPoint(SUNDomEigEstimator DEE,
+                                                         N_Vector Fyt)
+{
+  SUNErrCode ier;
+  SUNDIALS_MARK_FUNCTION_BEGIN(getSUNProfiler(DEE));
+  if (DEE->ops->setrhsatlinearizationpoint)
+  {
+    ier = DEE->ops->setrhsatlinearizationpoint(DEE, Fyt);
+  }
+  else { ier = SUN_SUCCESS; }
+  SUNDIALS_MARK_FUNCTION_END(getSUNProfiler(DEE));
+  return (ier);
+}
+
 SUNErrCode SUNDomEigEstimator_SetOptions(SUNDomEigEstimator DEE,
                                          const char* Did, const char* file_name,
                                          int argc, char* argv[])
diff --git a/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.c b/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.c
index 37994780aa..d960c81b9b 100644
--- a/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.c
+++ b/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.c
@@ -372,30 +372,6 @@ SWIGEXPORT void * _wrap_SUNDomEigEstimatorContent_Arnoldi__V_get(SwigClassWrappe
 }
 
 
-SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__q_set(SwigClassWrapper const *farg1, N_Vector farg2) {
-  struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
-  N_Vector arg2 = (N_Vector) 0 ;
-  
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::q", return );
-  arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr);
-  arg2 = (N_Vector)(farg2);
-  if (arg1) (arg1)->q = arg2;
-}
-
-
-SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Arnoldi__q_get(SwigClassWrapper const *farg1) {
-  N_Vector fresult ;
-  struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
-  N_Vector result;
-  
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::q", return 0);
-  arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr);
-  result = (N_Vector) ((arg1)->q);
-  fresult = result;
-  return fresult;
-}
-
-
 SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set(SwigClassWrapper const *farg1, N_Vector farg2) {
   struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
   N_Vector arg2 = (N_Vector) 0 ;
@@ -564,6 +540,30 @@ SWIGEXPORT int _wrap_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get(SwigCl
 }
 
 
+SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set(SwigClassWrapper const *farg1, int const *farg2) {
+  struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
+  int arg2 ;
+  
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::Fy_is_current", return );
+  arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr);
+  arg2 = (int)(*farg2);
+  if (arg1) (arg1)->Fy_is_current = arg2;
+}
+
+
+SWIGEXPORT int _wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get(SwigClassWrapper const *farg1) {
+  int fresult ;
+  struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
+  int result;
+  
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::Fy_is_current", return 0);
+  arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr);
+  result = (int) ((arg1)->Fy_is_current);
+  fresult = (int)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_set(SwigClassWrapper const *farg1, double const *farg2) {
   struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
   sunrealtype arg2 ;
@@ -971,6 +971,20 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(SUNDom
 }
 
 
+SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(SUNDomEigEstimator farg1, N_Vector farg2) {
+  int fresult ;
+  SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
+  N_Vector arg2 = (N_Vector) 0 ;
+  SUNErrCode result;
+  
+  arg1 = (SUNDomEigEstimator)(farg1);
+  arg2 = (N_Vector)(farg2);
+  result = (SUNErrCode)SUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(arg1,arg2);
+  fresult = (SUNErrCode)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(SUNDomEigEstimator farg1, int const *farg2) {
   int fresult ;
   SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
diff --git a/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.f90 b/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.f90
index 8edafe1f65..4a7318c844 100644
--- a/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.f90
+++ b/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.f90
@@ -46,8 +46,6 @@ module fsundomeigest_arnoldi_mod
   procedure :: get_ATdata => swigf_SUNDomEigEstimatorContent_Arnoldi__ATdata_get
   procedure :: set_V => swigf_SUNDomEigEstimatorContent_Arnoldi__V_set
   procedure :: get_V => swigf_SUNDomEigEstimatorContent_Arnoldi__V_get
-  procedure :: set_q => swigf_SUNDomEigEstimatorContent_Arnoldi__q_set
-  procedure :: get_q => swigf_SUNDomEigEstimatorContent_Arnoldi__q_get
   procedure :: set_rhs_linY => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set
   procedure :: get_rhs_linY => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_get
   procedure :: set_Fy => swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_set
@@ -62,6 +60,8 @@ module fsundomeigest_arnoldi_mod
   procedure :: get_num_iters => swigf_SUNDomEigEstimatorContent_Arnoldi__num_iters_get
   procedure :: set_warmup_to_tol => swigf_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_set
   procedure :: get_warmup_to_tol => swigf_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get
+  procedure :: set_Fy_is_current => swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set
+  procedure :: get_Fy_is_current => swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get
   procedure :: set_tol_warmup => swigf_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_set
   procedure :: get_tol_warmup => swigf_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_get
   procedure :: set_rhs_linT => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set
@@ -99,6 +99,7 @@ module fsundomeigest_arnoldi_mod
  public :: FSUNDomEigEstimator_SetATimes_Arnoldi
  public :: FSUNDomEigEstimator_SetRhs_Arnoldi
  public :: FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi
+ public :: FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi
  public :: FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi
  public :: FSUNDomEigEstimator_SetRelTol_Arnoldi
  public :: FSUNDomEigEstimator_SetInitialGuess_Arnoldi
@@ -163,23 +164,6 @@ function swigc_SUNDomEigEstimatorContent_Arnoldi__V_get(farg1) &
 type(C_PTR) :: fresult
 end function
 
-subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__q_set(farg1, farg2) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__q_set")
-use, intrinsic :: ISO_C_BINDING
-import :: swigclasswrapper
-type(SwigClassWrapper) :: farg1
-type(C_PTR), value :: farg2
-end subroutine
-
-function swigc_SUNDomEigEstimatorContent_Arnoldi__q_get(farg1) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__q_get") &
-result(fresult)
-use, intrinsic :: ISO_C_BINDING
-import :: swigclasswrapper
-type(SwigClassWrapper) :: farg1
-type(C_PTR) :: fresult
-end function
-
 subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set(farg1, farg2) &
 bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set")
 use, intrinsic :: ISO_C_BINDING
@@ -299,6 +283,23 @@ function swigc_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get(farg1) &
 integer(C_INT) :: fresult
 end function
 
+subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set(farg1, farg2) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set")
+use, intrinsic :: ISO_C_BINDING
+import :: swigclasswrapper
+type(SwigClassWrapper) :: farg1
+integer(C_INT), intent(in) :: farg2
+end subroutine
+
+function swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get(farg1) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+import :: swigclasswrapper
+type(SwigClassWrapper) :: farg1
+integer(C_INT) :: fresult
+end function
+
 subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_set(farg1, farg2) &
 bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_set")
 use, intrinsic :: ISO_C_BINDING
@@ -583,6 +584,15 @@ function swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(farg1, farg2
 integer(C_INT) :: fresult
 end function
 
+function swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(farg1, farg2) &
+bind(C, name="_wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+type(C_PTR), value :: farg1
+type(C_PTR), value :: farg2
+integer(C_INT) :: fresult
+end function
+
 function swigc_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(farg1, farg2) &
 bind(C, name="_wrap_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi") &
 result(fresult)
@@ -752,31 +762,6 @@ function swigf_SUNDomEigEstimatorContent_Arnoldi__V_get(self) &
 swig_result = fresult
 end function
 
-subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__q_set(self, q)
-use, intrinsic :: ISO_C_BINDING
-class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
-type(N_Vector), target, intent(inout) :: q
-type(SwigClassWrapper) :: farg1 
-type(C_PTR) :: farg2 
-
-farg1 = self%swigdata
-farg2 = c_loc(q)
-call swigc_SUNDomEigEstimatorContent_Arnoldi__q_set(farg1, farg2)
-end subroutine
-
-function swigf_SUNDomEigEstimatorContent_Arnoldi__q_get(self) &
-result(swig_result)
-use, intrinsic :: ISO_C_BINDING
-type(N_Vector), pointer :: swig_result
-class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
-type(C_PTR) :: fresult 
-type(SwigClassWrapper) :: farg1 
-
-farg1 = self%swigdata
-fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__q_get(farg1)
-call c_f_pointer(fresult, swig_result)
-end function
-
 subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set(self, rhs_liny)
 use, intrinsic :: ISO_C_BINDING
 class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
@@ -952,6 +937,31 @@ function swigf_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get(self) &
 swig_result = fresult
 end function
 
+subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set(self, fy_is_current)
+use, intrinsic :: ISO_C_BINDING
+class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
+integer(C_INT), intent(in) :: fy_is_current
+type(SwigClassWrapper) :: farg1 
+integer(C_INT) :: farg2 
+
+farg1 = self%swigdata
+farg2 = fy_is_current
+call swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set(farg1, farg2)
+end subroutine
+
+function swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get(self) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT) :: swig_result
+class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
+integer(C_INT) :: fresult 
+type(SwigClassWrapper) :: farg1 
+
+farg1 = self%swigdata
+fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get(farg1)
+swig_result = fresult
+end function
+
 subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_set(self, tol_warmup)
 use, intrinsic :: ISO_C_BINDING
 class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
@@ -1390,6 +1400,22 @@ function FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(dee, t, v) &
 swig_result = fresult
 end function
 
+function FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(dee, fyt) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT) :: swig_result
+type(SUNDomEigEstimator), target, intent(inout) :: dee
+type(N_Vector), target, intent(inout) :: fyt
+integer(C_INT) :: fresult 
+type(C_PTR) :: farg1 
+type(C_PTR) :: farg2 
+
+farg1 = c_loc(dee)
+farg2 = c_loc(fyt)
+fresult = swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(farg1, farg2)
+swig_result = fresult
+end function
+
 function FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(dee, num_iters) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
diff --git a/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.c b/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.c
index 6cc33c49a7..b81d4dde90 100644
--- a/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.c
+++ b/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.c
@@ -372,30 +372,6 @@ SWIGEXPORT void * _wrap_SUNDomEigEstimatorContent_Arnoldi__V_get(SwigClassWrappe
 }
 
 
-SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__q_set(SwigClassWrapper const *farg1, N_Vector farg2) {
-  struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
-  N_Vector arg2 = (N_Vector) 0 ;
-  
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::q", return );
-  arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr);
-  arg2 = (N_Vector)(farg2);
-  if (arg1) (arg1)->q = arg2;
-}
-
-
-SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Arnoldi__q_get(SwigClassWrapper const *farg1) {
-  N_Vector fresult ;
-  struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
-  N_Vector result;
-  
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::q", return 0);
-  arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr);
-  result = (N_Vector) ((arg1)->q);
-  fresult = result;
-  return fresult;
-}
-
-
 SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set(SwigClassWrapper const *farg1, N_Vector farg2) {
   struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
   N_Vector arg2 = (N_Vector) 0 ;
@@ -564,6 +540,30 @@ SWIGEXPORT int _wrap_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get(SwigCl
 }
 
 
+SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set(SwigClassWrapper const *farg1, int const *farg2) {
+  struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
+  int arg2 ;
+  
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::Fy_is_current", return );
+  arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr);
+  arg2 = (int)(*farg2);
+  if (arg1) (arg1)->Fy_is_current = arg2;
+}
+
+
+SWIGEXPORT int _wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get(SwigClassWrapper const *farg1) {
+  int fresult ;
+  struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
+  int result;
+  
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::Fy_is_current", return 0);
+  arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr);
+  result = (int) ((arg1)->Fy_is_current);
+  fresult = (int)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_set(SwigClassWrapper const *farg1, double const *farg2) {
   struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ;
   sunrealtype arg2 ;
@@ -971,6 +971,20 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(SUNDom
 }
 
 
+SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(SUNDomEigEstimator farg1, N_Vector farg2) {
+  int fresult ;
+  SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
+  N_Vector arg2 = (N_Vector) 0 ;
+  SUNErrCode result;
+  
+  arg1 = (SUNDomEigEstimator)(farg1);
+  arg2 = (N_Vector)(farg2);
+  result = (SUNErrCode)SUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(arg1,arg2);
+  fresult = (SUNErrCode)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(SUNDomEigEstimator farg1, int const *farg2) {
   int fresult ;
   SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
diff --git a/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.f90 b/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.f90
index 85861d0a3c..17f9171dbd 100644
--- a/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.f90
+++ b/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.f90
@@ -46,8 +46,6 @@ module fsundomeigest_arnoldi_mod
   procedure :: get_ATdata => swigf_SUNDomEigEstimatorContent_Arnoldi__ATdata_get
   procedure :: set_V => swigf_SUNDomEigEstimatorContent_Arnoldi__V_set
   procedure :: get_V => swigf_SUNDomEigEstimatorContent_Arnoldi__V_get
-  procedure :: set_q => swigf_SUNDomEigEstimatorContent_Arnoldi__q_set
-  procedure :: get_q => swigf_SUNDomEigEstimatorContent_Arnoldi__q_get
   procedure :: set_rhs_linY => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set
   procedure :: get_rhs_linY => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_get
   procedure :: set_Fy => swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_set
@@ -62,6 +60,8 @@ module fsundomeigest_arnoldi_mod
   procedure :: get_num_iters => swigf_SUNDomEigEstimatorContent_Arnoldi__num_iters_get
   procedure :: set_warmup_to_tol => swigf_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_set
   procedure :: get_warmup_to_tol => swigf_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get
+  procedure :: set_Fy_is_current => swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set
+  procedure :: get_Fy_is_current => swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get
   procedure :: set_tol_warmup => swigf_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_set
   procedure :: get_tol_warmup => swigf_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_get
   procedure :: set_rhs_linT => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set
@@ -99,6 +99,7 @@ module fsundomeigest_arnoldi_mod
  public :: FSUNDomEigEstimator_SetATimes_Arnoldi
  public :: FSUNDomEigEstimator_SetRhs_Arnoldi
  public :: FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi
+ public :: FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi
  public :: FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi
  public :: FSUNDomEigEstimator_SetRelTol_Arnoldi
  public :: FSUNDomEigEstimator_SetInitialGuess_Arnoldi
@@ -163,23 +164,6 @@ function swigc_SUNDomEigEstimatorContent_Arnoldi__V_get(farg1) &
 type(C_PTR) :: fresult
 end function
 
-subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__q_set(farg1, farg2) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__q_set")
-use, intrinsic :: ISO_C_BINDING
-import :: swigclasswrapper
-type(SwigClassWrapper) :: farg1
-type(C_PTR), value :: farg2
-end subroutine
-
-function swigc_SUNDomEigEstimatorContent_Arnoldi__q_get(farg1) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__q_get") &
-result(fresult)
-use, intrinsic :: ISO_C_BINDING
-import :: swigclasswrapper
-type(SwigClassWrapper) :: farg1
-type(C_PTR) :: fresult
-end function
-
 subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set(farg1, farg2) &
 bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set")
 use, intrinsic :: ISO_C_BINDING
@@ -299,6 +283,23 @@ function swigc_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get(farg1) &
 integer(C_INT) :: fresult
 end function
 
+subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set(farg1, farg2) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set")
+use, intrinsic :: ISO_C_BINDING
+import :: swigclasswrapper
+type(SwigClassWrapper) :: farg1
+integer(C_INT), intent(in) :: farg2
+end subroutine
+
+function swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get(farg1) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+import :: swigclasswrapper
+type(SwigClassWrapper) :: farg1
+integer(C_INT) :: fresult
+end function
+
 subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_set(farg1, farg2) &
 bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_set")
 use, intrinsic :: ISO_C_BINDING
@@ -583,6 +584,15 @@ function swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(farg1, farg2
 integer(C_INT) :: fresult
 end function
 
+function swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(farg1, farg2) &
+bind(C, name="_wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+type(C_PTR), value :: farg1
+type(C_PTR), value :: farg2
+integer(C_INT) :: fresult
+end function
+
 function swigc_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(farg1, farg2) &
 bind(C, name="_wrap_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi") &
 result(fresult)
@@ -752,31 +762,6 @@ function swigf_SUNDomEigEstimatorContent_Arnoldi__V_get(self) &
 swig_result = fresult
 end function
 
-subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__q_set(self, q)
-use, intrinsic :: ISO_C_BINDING
-class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
-type(N_Vector), target, intent(inout) :: q
-type(SwigClassWrapper) :: farg1 
-type(C_PTR) :: farg2 
-
-farg1 = self%swigdata
-farg2 = c_loc(q)
-call swigc_SUNDomEigEstimatorContent_Arnoldi__q_set(farg1, farg2)
-end subroutine
-
-function swigf_SUNDomEigEstimatorContent_Arnoldi__q_get(self) &
-result(swig_result)
-use, intrinsic :: ISO_C_BINDING
-type(N_Vector), pointer :: swig_result
-class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
-type(C_PTR) :: fresult 
-type(SwigClassWrapper) :: farg1 
-
-farg1 = self%swigdata
-fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__q_get(farg1)
-call c_f_pointer(fresult, swig_result)
-end function
-
 subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set(self, rhs_liny)
 use, intrinsic :: ISO_C_BINDING
 class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
@@ -952,6 +937,31 @@ function swigf_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get(self) &
 swig_result = fresult
 end function
 
+subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set(self, fy_is_current)
+use, intrinsic :: ISO_C_BINDING
+class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
+integer(C_INT), intent(in) :: fy_is_current
+type(SwigClassWrapper) :: farg1 
+integer(C_INT) :: farg2 
+
+farg1 = self%swigdata
+farg2 = fy_is_current
+call swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_set(farg1, farg2)
+end subroutine
+
+function swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get(self) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT) :: swig_result
+class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
+integer(C_INT) :: fresult 
+type(SwigClassWrapper) :: farg1 
+
+farg1 = self%swigdata
+fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_is_current_get(farg1)
+swig_result = fresult
+end function
+
 subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__tol_warmup_set(self, tol_warmup)
 use, intrinsic :: ISO_C_BINDING
 class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self
@@ -1390,6 +1400,22 @@ function FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(dee, t, v) &
 swig_result = fresult
 end function
 
+function FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(dee, fyt) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT) :: swig_result
+type(SUNDomEigEstimator), target, intent(inout) :: dee
+type(N_Vector), target, intent(inout) :: fyt
+integer(C_INT) :: fresult 
+type(C_PTR) :: farg1 
+type(C_PTR) :: farg2 
+
+farg1 = c_loc(dee)
+farg2 = c_loc(fyt)
+fresult = swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(farg1, farg2)
+swig_result = fresult
+end function
+
 function FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(dee, num_iters) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
diff --git a/src/sundomeigest/arnoldi/sundomeigest_arnoldi.c b/src/sundomeigest/arnoldi/sundomeigest_arnoldi.c
index 1f0cdd7c9b..034adcd867 100644
--- a/src/sundomeigest/arnoldi/sundomeigest_arnoldi.c
+++ b/src/sundomeigest/arnoldi/sundomeigest_arnoldi.c
@@ -113,6 +113,8 @@ SUNDomEigEstimator SUNDomEigEstimator_Arnoldi(N_Vector q, int kry_dim,
   DEE->ops->setrhs    = SUNDomEigEstimator_SetRhs_Arnoldi;
   DEE->ops->setrhslinearizationpoint =
     SUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi;
+  DEE->ops->setrhsatlinearizationpoint =
+    SUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi;
   DEE->ops->setnumpreprocessiters =
     SUNDomEigEstimator_SetNumPreprocessIters_Arnoldi;
   DEE->ops->setreltol         = SUNDomEigEstimator_SetRelTol_Arnoldi;
@@ -137,7 +139,6 @@ SUNDomEigEstimator SUNDomEigEstimator_Arnoldi(N_Vector q, int kry_dim,
   content->ATimes        = NULL;
   content->ATdata        = NULL;
   content->V             = NULL;
-  content->q             = NULL;
   content->rhs_linY      = NULL;
   content->rhs_linT      = ZERO;
   content->Fy            = NULL;
@@ -147,6 +148,7 @@ SUNDomEigEstimator SUNDomEigEstimator_Arnoldi(N_Vector q, int kry_dim,
   content->num_iters     = 0;
   content->num_ATimes    = 0;
   content->warmup_to_tol = SUNFALSE;
+  content->Fy_is_current = SUNFALSE;
   content->tol_warmup    = DEE_TOL_OF_WARMUPS_ARNOLDI_DEFAULT;
   content->rhsfn         = NULL;
   content->rhs_data      = NULL;
@@ -159,10 +161,6 @@ SUNDomEigEstimator SUNDomEigEstimator_Arnoldi(N_Vector q, int kry_dim,
   content->LAPACK_arr    = NULL;
   content->Hes           = NULL;
 
-  /* Allocate content */
-  content->q = N_VClone(q);
-  SUNCheckLastErrNull();
-
   content->V = N_VCloneVectorArray(kry_dim + 1, q);
   SUNCheckLastErrNull();
 
@@ -242,71 +240,88 @@ SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(
   return SUN_SUCCESS;
 }
 
-SUNErrCode SUNDomEigEstimator_Initialize_Arnoldi(SUNDomEigEstimator DEE)
+SUNErrCode SUNDomEigEstimator_SetRhsAtLinearizationPoint_Arnoldi(
+  SUNDomEigEstimator DEE, N_Vector Fyt)
 {
   SUNFunctionBegin(DEE->sunctx);
 
   SUNAssert(DEE, SUN_ERR_ARG_CORRUPT);
   SUNAssert(Arnoldi_CONTENT(DEE), SUN_ERR_ARG_CORRUPT);
+  SUNAssert(Fyt, SUN_ERR_ARG_CORRUPT);
 
-  if (Arnoldi_CONTENT(DEE)->kry_dim < 2)
-  {
-    Arnoldi_CONTENT(DEE)->kry_dim = DEE_KRYLOV_DIM_DEFAULT;
-  }
-  if (Arnoldi_CONTENT(DEE)->num_warmups < 0)
+  if (Arnoldi_CONTENT(DEE)->Fy == NULL)
   {
-    Arnoldi_CONTENT(DEE)->num_warmups = DEE_NUM_OF_WARMUPS_ARNOLDI_DEFAULT;
+    Arnoldi_CONTENT(DEE)->Fy = N_VClone(Fyt);
+    SUNCheckLastErr();
   }
 
+  N_VScale(ONE, Fyt, Arnoldi_CONTENT(DEE)->Fy);
+  SUNCheckLastErr();
+
+  Arnoldi_CONTENT(DEE)->Fy_is_current = SUNTRUE;
+
+  return SUN_SUCCESS;
+}
+
+SUNErrCode SUNDomEigEstimator_Initialize_Arnoldi(SUNDomEigEstimator DEE)
+{
+  SUNFunctionBegin(DEE->sunctx);
+
+  SUNAssert(DEE, SUN_ERR_ARG_CORRUPT);
+  SUNAssert(Arnoldi_CONTENT(DEE), SUN_ERR_ARG_CORRUPT);
   SUNAssert(Arnoldi_CONTENT(DEE)->ATimes, SUN_ERR_ARG_CORRUPT);
   SUNAssert(Arnoldi_CONTENT(DEE)->V, SUN_ERR_ARG_CORRUPT);
-  SUNAssert(Arnoldi_CONTENT(DEE)->q, SUN_ERR_ARG_CORRUPT);
+
+  int* kry_dim     = &(Arnoldi_CONTENT(DEE)->kry_dim);
+  int* num_warmups = &(Arnoldi_CONTENT(DEE)->num_warmups);
+
+  if (*kry_dim < 2) { *kry_dim = DEE_KRYLOV_DIM_DEFAULT; }
+  if (*num_warmups < 0) { *num_warmups = DEE_NUM_OF_WARMUPS_ARNOLDI_DEFAULT; }
 
   if (Arnoldi_CONTENT(DEE)->LAPACK_A == NULL)
   {
-    Arnoldi_CONTENT(DEE)->LAPACK_A = (sunrealtype*)malloc(
-      (Arnoldi_CONTENT(DEE)->kry_dim * Arnoldi_CONTENT(DEE)->kry_dim) *
-      sizeof(sunrealtype));
+    Arnoldi_CONTENT(DEE)->LAPACK_A =
+      (sunrealtype*)malloc((*kry_dim * *kry_dim) * sizeof(sunrealtype));
     SUNAssert(Arnoldi_CONTENT(DEE)->LAPACK_A, SUN_ERR_MALLOC_FAIL);
   }
   if (Arnoldi_CONTENT(DEE)->LAPACK_wr == NULL)
   {
-    Arnoldi_CONTENT(DEE)->LAPACK_wr =
-      malloc(Arnoldi_CONTENT(DEE)->kry_dim * sizeof(sunrealtype));
+    Arnoldi_CONTENT(DEE)->LAPACK_wr = malloc(*kry_dim * sizeof(sunrealtype));
     SUNAssert(Arnoldi_CONTENT(DEE)->LAPACK_wr, SUN_ERR_MALLOC_FAIL);
   }
   if (Arnoldi_CONTENT(DEE)->LAPACK_wi == NULL)
   {
-    Arnoldi_CONTENT(DEE)->LAPACK_wi =
-      malloc(Arnoldi_CONTENT(DEE)->kry_dim * sizeof(sunrealtype));
+    Arnoldi_CONTENT(DEE)->LAPACK_wi = malloc(*kry_dim * sizeof(sunrealtype));
     SUNAssert(Arnoldi_CONTENT(DEE)->LAPACK_wi, SUN_ERR_MALLOC_FAIL);
   }
 
   /* query the workspace size (call with lwork = -1) */
   char jobvl         = 'N';
   char jobvr         = 'N';
-  sunindextype N     = Arnoldi_CONTENT(DEE)->kry_dim;
-  sunindextype lda   = Arnoldi_CONTENT(DEE)->kry_dim;
+  sunindextype N     = *kry_dim;
+  sunindextype lda   = *kry_dim;
   sunindextype ldvl  = 1;
   sunindextype ldvr  = 1;
   sunindextype info  = 0;
   sunindextype lwork = -1;
   sunrealtype work   = SUN_RCONST(0.0);
+  sunrealtype* A     = Arnoldi_CONTENT(DEE)->LAPACK_A;
+  sunrealtype* wr    = Arnoldi_CONTENT(DEE)->LAPACK_wr;
+  sunrealtype* wi    = Arnoldi_CONTENT(DEE)->LAPACK_wi;
 
-  xgeev_f77(&jobvl, &jobvr, &N, Arnoldi_CONTENT(DEE)->LAPACK_A, &lda,
-            Arnoldi_CONTENT(DEE)->LAPACK_wr, Arnoldi_CONTENT(DEE)->LAPACK_wi,
-            NULL, &ldvl, NULL, &ldvr, &work, &lwork, &info);
+  xgeev_f77(&jobvl, &jobvr, &N, A, &lda, wr, wi, NULL, &ldvl, NULL, &ldvr,
+            &work, &lwork, &info);
 
   /* The workspace size is returned as the first entry of the work array */
   Arnoldi_CONTENT(DEE)->LAPACK_lwork = (sunindextype)work;
 
-  Arnoldi_CONTENT(DEE)->LAPACK_work = (sunrealtype*)malloc(
-    Arnoldi_CONTENT(DEE)->LAPACK_lwork * sizeof(sunrealtype));
+  Arnoldi_CONTENT(DEE)->LAPACK_work =
+    (sunrealtype*)malloc(((sunindextype)work) * sizeof(sunrealtype));
   SUNAssert(Arnoldi_CONTENT(DEE)->LAPACK_work, SUN_ERR_MALLOC_FAIL);
 
   /* LAPACK array */
   Arnoldi_CONTENT(DEE)->LAPACK_arr =
-    (sunrealtype**)malloc(Arnoldi_CONTENT(DEE)->kry_dim * sizeof(sunrealtype*));
+    (sunrealtype**)malloc(*kry_dim * sizeof(sunrealtype*));
   SUNAssert(Arnoldi_CONTENT(DEE)->LAPACK_arr, SUN_ERR_MALLOC_FAIL);
 
   for (int k = 0; k < Arnoldi_CONTENT(DEE)->kry_dim; k++)
@@ -317,14 +332,14 @@ SUNErrCode SUNDomEigEstimator_Initialize_Arnoldi(SUNDomEigEstimator DEE)
   }
 
   /* Hessenberg matrix Hes */
-  Arnoldi_CONTENT(DEE)->Hes = (sunrealtype**)malloc(
-    (Arnoldi_CONTENT(DEE)->kry_dim + 1) * sizeof(sunrealtype*));
+  Arnoldi_CONTENT(DEE)->Hes =
+    (sunrealtype**)malloc((*kry_dim + 1) * sizeof(sunrealtype*));
   SUNAssert(Arnoldi_CONTENT(DEE)->Hes, SUN_ERR_MALLOC_FAIL);
 
-  for (int k = 0; k <= Arnoldi_CONTENT(DEE)->kry_dim; k++)
+  for (int k = 0; k <= *kry_dim; k++)
   {
     Arnoldi_CONTENT(DEE)->Hes[k] =
-      (sunrealtype*)malloc(Arnoldi_CONTENT(DEE)->kry_dim * sizeof(sunrealtype));
+      (sunrealtype*)malloc(*kry_dim * sizeof(sunrealtype));
     SUNAssert(Arnoldi_CONTENT(DEE)->Hes[k], SUN_ERR_MALLOC_FAIL);
   }
 
@@ -407,41 +422,53 @@ SUNErrCode SUNDomEigEstimator_Estimate_Arnoldi(SUNDomEigEstimator DEE,
   SUNAssert(lambdaI, SUN_ERR_ARG_CORRUPT);
   SUNAssert(Arnoldi_CONTENT(DEE)->ATimes, SUN_ERR_ARG_CORRUPT);
   SUNAssert(Arnoldi_CONTENT(DEE)->V, SUN_ERR_ARG_CORRUPT);
-  SUNAssert(Arnoldi_CONTENT(DEE)->q, SUN_ERR_ARG_CORRUPT);
   SUNAssert(Arnoldi_CONTENT(DEE)->Hes, SUN_ERR_ARG_CORRUPT);
 
   int retval;
-  sunindextype n = Arnoldi_CONTENT(DEE)->kry_dim;
-  sunrealtype normq;
-  Arnoldi_CONTENT(DEE)->num_ATimes = 0;
-  Arnoldi_CONTENT(DEE)->num_iters  = 0;
+  int* num_warmups = &(Arnoldi_CONTENT(DEE)->num_warmups);
+  sunindextype n   = Arnoldi_CONTENT(DEE)->kry_dim;
+  sunrealtype normAv;
+  long int* num_ATimes = &(Arnoldi_CONTENT(DEE)->num_ATimes);
+  long int* num_iters  = &(Arnoldi_CONTENT(DEE)->num_iters);
+  *num_ATimes          = 0;
+  *num_iters           = 0;
+
+  N_Vector* V = Arnoldi_CONTENT(DEE)->V;
+  N_Vector Av = V[1];
+
+  sunrealtype** Hes = Arnoldi_CONTENT(DEE)->Hes;
+
+  SUNATimesFn ATimes = Arnoldi_CONTENT(DEE)->ATimes;
+  void* ATdata       = Arnoldi_CONTENT(DEE)->ATdata;
+
+  sunrealtype* A    = Arnoldi_CONTENT(DEE)->LAPACK_A;
+  sunrealtype* wr   = Arnoldi_CONTENT(DEE)->LAPACK_wr;
+  sunrealtype* wi   = Arnoldi_CONTENT(DEE)->LAPACK_wi;
+  sunrealtype** arr = Arnoldi_CONTENT(DEE)->LAPACK_arr;
 
   sunrealtype res;
   sunrealtype new_lambda = ZERO;
   sunrealtype old_lambda = ZERO;
 
   /* Set the initial q = A^{num_warmups}q/||A^{num_warmups}q|| */
-  for (int i = 0; i < Arnoldi_CONTENT(DEE)->num_warmups; i++)
+  for (int i = 0; i < *num_warmups; i++)
   {
-    retval = Arnoldi_CONTENT(DEE)->ATimes(Arnoldi_CONTENT(DEE)->ATdata,
-                                          Arnoldi_CONTENT(DEE)->V[0],
-                                          Arnoldi_CONTENT(DEE)->q);
-    Arnoldi_CONTENT(DEE)->num_ATimes++;
-    Arnoldi_CONTENT(DEE)->num_iters++;
+    retval = ATimes(ATdata, V[0], Av);
+    (*num_ATimes)++;
+    (*num_iters)++;
     if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; }
 
     if (Arnoldi_CONTENT(DEE)->warmup_to_tol)
     {
-      new_lambda = N_VDotProd(Arnoldi_CONTENT(DEE)->V[0],
-                              Arnoldi_CONTENT(DEE)->q); //Rayleigh quotient
+      new_lambda = N_VDotProd(V[0], Av); //Rayleigh quotient
       SUNCheckLastErr();
     }
 
-    normq = N_VDotProd(Arnoldi_CONTENT(DEE)->q, Arnoldi_CONTENT(DEE)->q);
+    normAv = N_VDotProd(Av, Av);
     SUNCheckLastErr();
 
-    normq = SUNRsqrt(normq);
-    N_VScale(ONE / normq, Arnoldi_CONTENT(DEE)->q, Arnoldi_CONTENT(DEE)->V[0]);
+    normAv = SUNRsqrt(normAv);
+    N_VScale(ONE / normAv, Av, V[0]);
     SUNCheckLastErr();
 
     if (Arnoldi_CONTENT(DEE)->warmup_to_tol)
@@ -458,20 +485,15 @@ SUNErrCode SUNDomEigEstimator_Estimate_Arnoldi(SUNDomEigEstimator DEE,
   for (int i = 0; i < n; i++)
   {
     /* Compute the next Krylov vector */
-    retval = Arnoldi_CONTENT(DEE)->ATimes(Arnoldi_CONTENT(DEE)->ATdata,
-                                          Arnoldi_CONTENT(DEE)->V[i],
-                                          Arnoldi_CONTENT(DEE)->V[i + 1]);
-    Arnoldi_CONTENT(DEE)->num_ATimes++;
-    Arnoldi_CONTENT(DEE)->num_iters++;
+    retval = ATimes(ATdata, V[i], V[i + 1]);
+    (*num_ATimes)++;
+    (*num_iters)++;
     if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; }
 
-    SUNCheckCall(SUNModifiedGS(Arnoldi_CONTENT(DEE)->V,
-                               Arnoldi_CONTENT(DEE)->Hes, i + 1, (int)n,
-                               &(Arnoldi_CONTENT(DEE)->Hes[i + 1][i])));
+    SUNCheckCall(SUNModifiedGS(V, Hes, i + 1, (int)n, &(Hes[i + 1][i])));
 
     /* Unitize the computed orthogonal vector */
-    N_VScale(SUN_RCONST(1.0) / Arnoldi_CONTENT(DEE)->Hes[i + 1][i],
-             Arnoldi_CONTENT(DEE)->V[i + 1], Arnoldi_CONTENT(DEE)->V[i + 1]);
+    N_VScale(SUN_RCONST(1.0) / Hes[i + 1][i], V[i + 1], V[i + 1]);
     SUNCheckLastErr();
   }
 
@@ -481,7 +503,7 @@ SUNErrCode SUNDomEigEstimator_Estimate_Arnoldi(SUNDomEigEstimator DEE,
   {
     for (int i = 0; i < n; i++)
     {
-      Arnoldi_CONTENT(DEE)->LAPACK_A[k] = Arnoldi_CONTENT(DEE)->Hes[i][j];
+      A[k] = Hes[i][j];
       k++;
     }
   }
@@ -503,34 +525,31 @@ SUNErrCode SUNDomEigEstimator_Estimate_Arnoldi(SUNDomEigEstimator DEE,
   sunindextype ldvr = n;
   sunindextype info;
   sunindextype lwork = Arnoldi_CONTENT(DEE)->LAPACK_lwork;
-  xgeev_f77(&jobvl, &jobvr, &n, Arnoldi_CONTENT(DEE)->LAPACK_A, &lda,
-            Arnoldi_CONTENT(DEE)->LAPACK_wr, Arnoldi_CONTENT(DEE)->LAPACK_wi,
-            NULL, &ldvl, NULL, &ldvr, Arnoldi_CONTENT(DEE)->LAPACK_work, &lwork,
-            &info);
+  xgeev_f77(&jobvl, &jobvr, &n, A, &lda, wr, wi, NULL, &ldvl, NULL, &ldvr,
+            Arnoldi_CONTENT(DEE)->LAPACK_work, &lwork, &info);
 
   if (info != 0) { return SUN_ERR_EXT_FAIL; }
 
   /* order the eigenvalues by their magnitude */
   for (int i = 0; i < n; i++)
   {
-    Arnoldi_CONTENT(DEE)->LAPACK_arr[i][0] = Arnoldi_CONTENT(DEE)->LAPACK_wr[i];
-    Arnoldi_CONTENT(DEE)->LAPACK_arr[i][1] = Arnoldi_CONTENT(DEE)->LAPACK_wi[i];
+    arr[i][0] = wr[i];
+    arr[i][1] = wi[i];
   }
 
   /* Sort the array using qsort */
-  qsort(Arnoldi_CONTENT(DEE)->LAPACK_arr, n,
-        sizeof(Arnoldi_CONTENT(DEE)->LAPACK_arr[0]), sundomeigest_Compare);
+  qsort(arr, n, sizeof(arr[0]), sundomeigest_Compare);
 
   /* Substitute the ordered eigenvalues back in LAPACK_w* */
   for (int i = 0; i < n; i++)
   {
-    Arnoldi_CONTENT(DEE)->LAPACK_wr[i] = Arnoldi_CONTENT(DEE)->LAPACK_arr[i][0];
-    Arnoldi_CONTENT(DEE)->LAPACK_wi[i] = Arnoldi_CONTENT(DEE)->LAPACK_arr[i][1];
+    wr[i] = arr[i][0];
+    wi[i] = arr[i][1];
   }
 
   /* Copy the dominant eigenvalue */
-  *lambdaR = Arnoldi_CONTENT(DEE)->LAPACK_wr[0];
-  *lambdaI = Arnoldi_CONTENT(DEE)->LAPACK_wi[0];
+  *lambdaR = wr[0];
+  *lambdaI = wi[0];
 
   return SUN_SUCCESS;
 }
@@ -607,11 +626,6 @@ SUNErrCode SUNDomEigEstimator_Destroy_Arnoldi(SUNDomEigEstimator* DEEptr)
   if (DEE->content)
   {
     /* delete items from within the content structure */
-    if (Arnoldi_CONTENT(DEE)->q)
-    {
-      N_VDestroy(Arnoldi_CONTENT(DEE)->q);
-      Arnoldi_CONTENT(DEE)->q = NULL;
-    }
     if (Arnoldi_CONTENT(DEE)->rhs_linY)
     {
       N_VDestroy(Arnoldi_CONTENT(DEE)->rhs_linY);
@@ -744,14 +758,26 @@ SUNErrCode dee_DQJtimes_Arnoldi(void* voidstarDEE, N_Vector v, N_Vector Jv)
   sunrealtype sig, siginv;
   int iter, retval;
 
-  N_Vector y    = Arnoldi_CONTENT(DEE)->rhs_linY;
-  N_Vector work = Arnoldi_CONTENT(DEE)->work;
-  N_Vector Fy   = Arnoldi_CONTENT(DEE)->Fy;
+  sunbooleantype* Fy_is_current = &(Arnoldi_CONTENT(DEE)->Fy_is_current);
+  N_Vector y                    = Arnoldi_CONTENT(DEE)->rhs_linY;
+  N_Vector work                 = Arnoldi_CONTENT(DEE)->work;
+  N_Vector Fy                   = Arnoldi_CONTENT(DEE)->Fy;
+
+  SUNRhsFn rhsfn = Arnoldi_CONTENT(DEE)->rhsfn;
+  void* rhs_data = Arnoldi_CONTENT(DEE)->rhs_data;
+
+  sunrealtype rhs_linT = Arnoldi_CONTENT(DEE)->rhs_linT;
+
+  long int* nfevals = &(Arnoldi_CONTENT(DEE)->nfevals);
 
-  retval = Arnoldi_CONTENT(DEE)->rhsfn(Arnoldi_CONTENT(DEE)->rhs_linT, y, Fy,
-                                       Arnoldi_CONTENT(DEE)->rhs_data);
-  Arnoldi_CONTENT(DEE)->nfevals++;
-  if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; }
+  if (!(*Fy_is_current))
+  {
+    /* If Fy is not current, compute it at the linearization point */
+    retval = rhsfn(rhs_linT, y, Fy, rhs_data);
+    (*nfevals)++;
+    if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; }
+    *Fy_is_current = SUNTRUE;
+  }
 
   /* Initialize perturbation */
   sunrealtype ydotv   = N_VDotProd(y, v);
@@ -766,9 +792,8 @@ SUNErrCode dee_DQJtimes_Arnoldi(void* voidstarDEE, N_Vector v, N_Vector Jv)
     N_VLinearSum(sig, v, ONE, y, work);
 
     /* Set Jv = f(tn, y+sig*v) */
-    retval = Arnoldi_CONTENT(DEE)->rhsfn(Arnoldi_CONTENT(DEE)->rhs_linT, work,
-                                         Jv, Arnoldi_CONTENT(DEE)->rhs_data);
-    Arnoldi_CONTENT(DEE)->nfevals++;
+    retval = rhsfn(rhs_linT, work, Jv, rhs_data);
+    (*nfevals)++;
     if (retval == 0) { break; }
     if (retval < 0) { return SUN_ERR_USER_FCN_FAIL; }
 
@@ -783,5 +808,6 @@ SUNErrCode dee_DQJtimes_Arnoldi(void* voidstarDEE, N_Vector v, N_Vector Jv)
   siginv = ONE / sig;
   N_VLinearSum(siginv, Jv, -siginv, Fy, Jv);
 
+  *Fy_is_current = SUNFALSE;
   return SUN_SUCCESS;
-}
\ No newline at end of file
+}
diff --git a/src/sundomeigest/power/fmod_int32/fsundomeigest_power_mod.c b/src/sundomeigest/power/fmod_int32/fsundomeigest_power_mod.c
index 324db92a27..6462b07889 100644
--- a/src/sundomeigest/power/fmod_int32/fsundomeigest_power_mod.c
+++ b/src/sundomeigest/power/fmod_int32/fsundomeigest_power_mod.c
@@ -372,49 +372,49 @@ SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Power__V_get(SwigClassWrappe
 }
 
 
-SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Power__q_set(SwigClassWrapper const *farg1, N_Vector farg2) {
+SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Power__Av_set(SwigClassWrapper const *farg1, N_Vector farg2) {
   struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
   N_Vector arg2 = (N_Vector) 0 ;
   
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::q", return );
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::Av", return );
   arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
   arg2 = (N_Vector)(farg2);
-  if (arg1) (arg1)->q = arg2;
+  if (arg1) (arg1)->Av = arg2;
 }
 
 
-SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Power__q_get(SwigClassWrapper const *farg1) {
+SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Power__Av_get(SwigClassWrapper const *farg1) {
   N_Vector fresult ;
   struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
   N_Vector result;
   
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::q", return 0);
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::Av", return 0);
   arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
-  result = (N_Vector) ((arg1)->q);
+  result = (N_Vector) ((arg1)->Av);
   fresult = result;
   return fresult;
 }
 
 
-SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Power__q_prev_set(SwigClassWrapper const *farg1, N_Vector farg2) {
+SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Power__v_prev_set(SwigClassWrapper const *farg1, N_Vector farg2) {
   struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
   N_Vector arg2 = (N_Vector) 0 ;
   
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::q_prev", return );
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::v_prev", return );
   arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
   arg2 = (N_Vector)(farg2);
-  if (arg1) (arg1)->q_prev = arg2;
+  if (arg1) (arg1)->v_prev = arg2;
 }
 
 
-SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Power__q_prev_get(SwigClassWrapper const *farg1) {
+SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Power__v_prev_get(SwigClassWrapper const *farg1) {
   N_Vector fresult ;
   struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
   N_Vector result;
   
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::q_prev", return 0);
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::v_prev", return 0);
   arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
-  result = (N_Vector) ((arg1)->q_prev);
+  result = (N_Vector) ((arg1)->v_prev);
   fresult = result;
   return fresult;
 }
@@ -756,6 +756,30 @@ SWIGEXPORT int _wrap_SUNDomEigEstimatorContent_Power__is_complex_get(SwigClassWr
 }
 
 
+SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Power__Fy_is_current_set(SwigClassWrapper const *farg1, int const *farg2) {
+  struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
+  int arg2 ;
+  
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::Fy_is_current", return );
+  arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
+  arg2 = (int)(*farg2);
+  if (arg1) (arg1)->Fy_is_current = arg2;
+}
+
+
+SWIGEXPORT int _wrap_SUNDomEigEstimatorContent_Power__Fy_is_current_get(SwigClassWrapper const *farg1) {
+  int fresult ;
+  struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
+  int result;
+  
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::Fy_is_current", return 0);
+  arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
+  result = (int) ((arg1)->Fy_is_current);
+  fresult = (int)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT SwigClassWrapper _wrap_new_SUNDomEigEstimatorContent_Power_() {
   SwigClassWrapper fresult ;
   struct SUNDomEigEstimatorContent_Power_ *result = 0 ;
@@ -909,6 +933,20 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint_Power(SUNDomEi
 }
 
 
+SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(SUNDomEigEstimator farg1, N_Vector farg2) {
+  int fresult ;
+  SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
+  N_Vector arg2 = (N_Vector) 0 ;
+  SUNErrCode result;
+  
+  arg1 = (SUNDomEigEstimator)(farg1);
+  arg2 = (N_Vector)(farg2);
+  result = (SUNErrCode)SUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(arg1,arg2);
+  fresult = (SUNErrCode)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetIsReal_Power(SUNDomEigEstimator farg1, int const *farg2) {
   int fresult ;
   SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
diff --git a/src/sundomeigest/power/fmod_int32/fsundomeigest_power_mod.f90 b/src/sundomeigest/power/fmod_int32/fsundomeigest_power_mod.f90
index 41b78edd03..5c38720a8f 100644
--- a/src/sundomeigest/power/fmod_int32/fsundomeigest_power_mod.f90
+++ b/src/sundomeigest/power/fmod_int32/fsundomeigest_power_mod.f90
@@ -46,10 +46,10 @@ module fsundomeigest_power_mod
   procedure :: get_ATdata => swigf_SUNDomEigEstimatorContent_Power__ATdata_get
   procedure :: set_V => swigf_SUNDomEigEstimatorContent_Power__V_set
   procedure :: get_V => swigf_SUNDomEigEstimatorContent_Power__V_get
-  procedure :: set_q => swigf_SUNDomEigEstimatorContent_Power__q_set
-  procedure :: get_q => swigf_SUNDomEigEstimatorContent_Power__q_get
-  procedure :: set_q_prev => swigf_SUNDomEigEstimatorContent_Power__q_prev_set
-  procedure :: get_q_prev => swigf_SUNDomEigEstimatorContent_Power__q_prev_get
+  procedure :: set_Av => swigf_SUNDomEigEstimatorContent_Power__Av_set
+  procedure :: get_Av => swigf_SUNDomEigEstimatorContent_Power__Av_get
+  procedure :: set_v_prev => swigf_SUNDomEigEstimatorContent_Power__v_prev_set
+  procedure :: get_v_prev => swigf_SUNDomEigEstimatorContent_Power__v_prev_get
   procedure :: set_rhs_linY => swigf_SUNDomEigEstimatorContent_Power__rhs_linY_set
   procedure :: get_rhs_linY => swigf_SUNDomEigEstimatorContent_Power__rhs_linY_get
   procedure :: set_Fy => swigf_SUNDomEigEstimatorContent_Power__Fy_set
@@ -78,6 +78,8 @@ module fsundomeigest_power_mod
   procedure :: get_nfevals => swigf_SUNDomEigEstimatorContent_Power__nfevals_get
   procedure :: set_is_complex => swigf_SUNDomEigEstimatorContent_Power__is_complex_set
   procedure :: get_is_complex => swigf_SUNDomEigEstimatorContent_Power__is_complex_get
+  procedure :: set_Fy_is_current => swigf_SUNDomEigEstimatorContent_Power__Fy_is_current_set
+  procedure :: get_Fy_is_current => swigf_SUNDomEigEstimatorContent_Power__Fy_is_current_get
   procedure :: release => swigf_release_SUNDomEigEstimatorContent_Power_
   procedure, private :: swigf_SUNDomEigEstimatorContent_Power__op_assign__
   generic :: assignment(=) => swigf_SUNDomEigEstimatorContent_Power__op_assign__
@@ -93,6 +95,7 @@ module fsundomeigest_power_mod
  public :: FSUNDomEigEstimator_SetRelTol_Power
  public :: FSUNDomEigEstimator_SetInitialGuess_Power
  public :: FSUNDomEigEstimator_SetRhsLinearizationPoint_Power
+ public :: FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power
  public :: FSUNDomEigEstimator_SetIsReal_Power
  public :: FSUNDomEigEstimator_Initialize_Power
  public :: FSUNDomEigEstimator_Estimate_Power
@@ -156,16 +159,16 @@ function swigc_SUNDomEigEstimatorContent_Power__V_get(farg1) &
 type(C_PTR) :: fresult
 end function
 
-subroutine swigc_SUNDomEigEstimatorContent_Power__q_set(farg1, farg2) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__q_set")
+subroutine swigc_SUNDomEigEstimatorContent_Power__Av_set(farg1, farg2) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__Av_set")
 use, intrinsic :: ISO_C_BINDING
 import :: swigclasswrapper
 type(SwigClassWrapper) :: farg1
 type(C_PTR), value :: farg2
 end subroutine
 
-function swigc_SUNDomEigEstimatorContent_Power__q_get(farg1) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__q_get") &
+function swigc_SUNDomEigEstimatorContent_Power__Av_get(farg1) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__Av_get") &
 result(fresult)
 use, intrinsic :: ISO_C_BINDING
 import :: swigclasswrapper
@@ -173,16 +176,16 @@ function swigc_SUNDomEigEstimatorContent_Power__q_get(farg1) &
 type(C_PTR) :: fresult
 end function
 
-subroutine swigc_SUNDomEigEstimatorContent_Power__q_prev_set(farg1, farg2) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__q_prev_set")
+subroutine swigc_SUNDomEigEstimatorContent_Power__v_prev_set(farg1, farg2) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__v_prev_set")
 use, intrinsic :: ISO_C_BINDING
 import :: swigclasswrapper
 type(SwigClassWrapper) :: farg1
 type(C_PTR), value :: farg2
 end subroutine
 
-function swigc_SUNDomEigEstimatorContent_Power__q_prev_get(farg1) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__q_prev_get") &
+function swigc_SUNDomEigEstimatorContent_Power__v_prev_get(farg1) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__v_prev_get") &
 result(fresult)
 use, intrinsic :: ISO_C_BINDING
 import :: swigclasswrapper
@@ -428,6 +431,23 @@ function swigc_SUNDomEigEstimatorContent_Power__is_complex_get(farg1) &
 integer(C_INT) :: fresult
 end function
 
+subroutine swigc_SUNDomEigEstimatorContent_Power__Fy_is_current_set(farg1, farg2) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__Fy_is_current_set")
+use, intrinsic :: ISO_C_BINDING
+import :: swigclasswrapper
+type(SwigClassWrapper) :: farg1
+integer(C_INT), intent(in) :: farg2
+end subroutine
+
+function swigc_SUNDomEigEstimatorContent_Power__Fy_is_current_get(farg1) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__Fy_is_current_get") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+import :: swigclasswrapper
+type(SwigClassWrapper) :: farg1
+integer(C_INT) :: fresult
+end function
+
 function swigc_new_SUNDomEigEstimatorContent_Power_() &
 bind(C, name="_wrap_new_SUNDomEigEstimatorContent_Power_") &
 result(fresult)
@@ -528,6 +548,15 @@ function swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint_Power(farg1, farg2,
 integer(C_INT) :: fresult
 end function
 
+function swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(farg1, farg2) &
+bind(C, name="_wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+type(C_PTR), value :: farg1
+type(C_PTR), value :: farg2
+integer(C_INT) :: fresult
+end function
+
 function swigc_FSUNDomEigEstimator_SetIsReal_Power(farg1, farg2) &
 bind(C, name="_wrap_FSUNDomEigEstimator_SetIsReal_Power") &
 result(fresult)
@@ -688,19 +717,19 @@ function swigf_SUNDomEigEstimatorContent_Power__V_get(self) &
 call c_f_pointer(fresult, swig_result)
 end function
 
-subroutine swigf_SUNDomEigEstimatorContent_Power__q_set(self, q)
+subroutine swigf_SUNDomEigEstimatorContent_Power__Av_set(self, av)
 use, intrinsic :: ISO_C_BINDING
 class(SUNDomEigEstimatorContent_Power_), intent(in) :: self
-type(N_Vector), target, intent(inout) :: q
+type(N_Vector), target, intent(inout) :: av
 type(SwigClassWrapper) :: farg1 
 type(C_PTR) :: farg2 
 
 farg1 = self%swigdata
-farg2 = c_loc(q)
-call swigc_SUNDomEigEstimatorContent_Power__q_set(farg1, farg2)
+farg2 = c_loc(av)
+call swigc_SUNDomEigEstimatorContent_Power__Av_set(farg1, farg2)
 end subroutine
 
-function swigf_SUNDomEigEstimatorContent_Power__q_get(self) &
+function swigf_SUNDomEigEstimatorContent_Power__Av_get(self) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
 type(N_Vector), pointer :: swig_result
@@ -709,23 +738,23 @@ function swigf_SUNDomEigEstimatorContent_Power__q_get(self) &
 type(SwigClassWrapper) :: farg1 
 
 farg1 = self%swigdata
-fresult = swigc_SUNDomEigEstimatorContent_Power__q_get(farg1)
+fresult = swigc_SUNDomEigEstimatorContent_Power__Av_get(farg1)
 call c_f_pointer(fresult, swig_result)
 end function
 
-subroutine swigf_SUNDomEigEstimatorContent_Power__q_prev_set(self, q_prev)
+subroutine swigf_SUNDomEigEstimatorContent_Power__v_prev_set(self, v_prev)
 use, intrinsic :: ISO_C_BINDING
 class(SUNDomEigEstimatorContent_Power_), intent(in) :: self
-type(N_Vector), target, intent(inout) :: q_prev
+type(N_Vector), target, intent(inout) :: v_prev
 type(SwigClassWrapper) :: farg1 
 type(C_PTR) :: farg2 
 
 farg1 = self%swigdata
-farg2 = c_loc(q_prev)
-call swigc_SUNDomEigEstimatorContent_Power__q_prev_set(farg1, farg2)
+farg2 = c_loc(v_prev)
+call swigc_SUNDomEigEstimatorContent_Power__v_prev_set(farg1, farg2)
 end subroutine
 
-function swigf_SUNDomEigEstimatorContent_Power__q_prev_get(self) &
+function swigf_SUNDomEigEstimatorContent_Power__v_prev_get(self) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
 type(N_Vector), pointer :: swig_result
@@ -734,7 +763,7 @@ function swigf_SUNDomEigEstimatorContent_Power__q_prev_get(self) &
 type(SwigClassWrapper) :: farg1 
 
 farg1 = self%swigdata
-fresult = swigc_SUNDomEigEstimatorContent_Power__q_prev_get(farg1)
+fresult = swigc_SUNDomEigEstimatorContent_Power__v_prev_get(farg1)
 call c_f_pointer(fresult, swig_result)
 end function
 
@@ -1088,6 +1117,31 @@ function swigf_SUNDomEigEstimatorContent_Power__is_complex_get(self) &
 swig_result = fresult
 end function
 
+subroutine swigf_SUNDomEigEstimatorContent_Power__Fy_is_current_set(self, fy_is_current)
+use, intrinsic :: ISO_C_BINDING
+class(SUNDomEigEstimatorContent_Power_), intent(in) :: self
+integer(C_INT), intent(in) :: fy_is_current
+type(SwigClassWrapper) :: farg1 
+integer(C_INT) :: farg2 
+
+farg1 = self%swigdata
+farg2 = fy_is_current
+call swigc_SUNDomEigEstimatorContent_Power__Fy_is_current_set(farg1, farg2)
+end subroutine
+
+function swigf_SUNDomEigEstimatorContent_Power__Fy_is_current_get(self) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT) :: swig_result
+class(SUNDomEigEstimatorContent_Power_), intent(in) :: self
+integer(C_INT) :: fresult 
+type(SwigClassWrapper) :: farg1 
+
+farg1 = self%swigdata
+fresult = swigc_SUNDomEigEstimatorContent_Power__Fy_is_current_get(farg1)
+swig_result = fresult
+end function
+
 function swigf_create_SUNDomEigEstimatorContent_Power_() &
 result(self)
 use, intrinsic :: ISO_C_BINDING
@@ -1268,6 +1322,22 @@ function FSUNDomEigEstimator_SetRhsLinearizationPoint_Power(dee, t, v) &
 swig_result = fresult
 end function
 
+function FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(dee, fyt) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT) :: swig_result
+type(SUNDomEigEstimator), target, intent(inout) :: dee
+type(N_Vector), target, intent(inout) :: fyt
+integer(C_INT) :: fresult 
+type(C_PTR) :: farg1 
+type(C_PTR) :: farg2 
+
+farg1 = c_loc(dee)
+farg2 = c_loc(fyt)
+fresult = swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(farg1, farg2)
+swig_result = fresult
+end function
+
 function FSUNDomEigEstimator_SetIsReal_Power(dee, real) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
diff --git a/src/sundomeigest/power/fmod_int64/fsundomeigest_power_mod.c b/src/sundomeigest/power/fmod_int64/fsundomeigest_power_mod.c
index 324db92a27..6462b07889 100644
--- a/src/sundomeigest/power/fmod_int64/fsundomeigest_power_mod.c
+++ b/src/sundomeigest/power/fmod_int64/fsundomeigest_power_mod.c
@@ -372,49 +372,49 @@ SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Power__V_get(SwigClassWrappe
 }
 
 
-SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Power__q_set(SwigClassWrapper const *farg1, N_Vector farg2) {
+SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Power__Av_set(SwigClassWrapper const *farg1, N_Vector farg2) {
   struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
   N_Vector arg2 = (N_Vector) 0 ;
   
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::q", return );
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::Av", return );
   arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
   arg2 = (N_Vector)(farg2);
-  if (arg1) (arg1)->q = arg2;
+  if (arg1) (arg1)->Av = arg2;
 }
 
 
-SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Power__q_get(SwigClassWrapper const *farg1) {
+SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Power__Av_get(SwigClassWrapper const *farg1) {
   N_Vector fresult ;
   struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
   N_Vector result;
   
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::q", return 0);
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::Av", return 0);
   arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
-  result = (N_Vector) ((arg1)->q);
+  result = (N_Vector) ((arg1)->Av);
   fresult = result;
   return fresult;
 }
 
 
-SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Power__q_prev_set(SwigClassWrapper const *farg1, N_Vector farg2) {
+SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Power__v_prev_set(SwigClassWrapper const *farg1, N_Vector farg2) {
   struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
   N_Vector arg2 = (N_Vector) 0 ;
   
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::q_prev", return );
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::v_prev", return );
   arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
   arg2 = (N_Vector)(farg2);
-  if (arg1) (arg1)->q_prev = arg2;
+  if (arg1) (arg1)->v_prev = arg2;
 }
 
 
-SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Power__q_prev_get(SwigClassWrapper const *farg1) {
+SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Power__v_prev_get(SwigClassWrapper const *farg1) {
   N_Vector fresult ;
   struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
   N_Vector result;
   
-  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::q_prev", return 0);
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::v_prev", return 0);
   arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
-  result = (N_Vector) ((arg1)->q_prev);
+  result = (N_Vector) ((arg1)->v_prev);
   fresult = result;
   return fresult;
 }
@@ -756,6 +756,30 @@ SWIGEXPORT int _wrap_SUNDomEigEstimatorContent_Power__is_complex_get(SwigClassWr
 }
 
 
+SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Power__Fy_is_current_set(SwigClassWrapper const *farg1, int const *farg2) {
+  struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
+  int arg2 ;
+  
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::Fy_is_current", return );
+  arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
+  arg2 = (int)(*farg2);
+  if (arg1) (arg1)->Fy_is_current = arg2;
+}
+
+
+SWIGEXPORT int _wrap_SUNDomEigEstimatorContent_Power__Fy_is_current_get(SwigClassWrapper const *farg1) {
+  int fresult ;
+  struct SUNDomEigEstimatorContent_Power_ *arg1 = (struct SUNDomEigEstimatorContent_Power_ *) 0 ;
+  int result;
+  
+  SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Power_ *", "SUNDomEigEstimatorContent_Power_", "SUNDomEigEstimatorContent_Power_::Fy_is_current", return 0);
+  arg1 = (struct SUNDomEigEstimatorContent_Power_ *)(farg1->cptr);
+  result = (int) ((arg1)->Fy_is_current);
+  fresult = (int)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT SwigClassWrapper _wrap_new_SUNDomEigEstimatorContent_Power_() {
   SwigClassWrapper fresult ;
   struct SUNDomEigEstimatorContent_Power_ *result = 0 ;
@@ -909,6 +933,20 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint_Power(SUNDomEi
 }
 
 
+SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(SUNDomEigEstimator farg1, N_Vector farg2) {
+  int fresult ;
+  SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
+  N_Vector arg2 = (N_Vector) 0 ;
+  SUNErrCode result;
+  
+  arg1 = (SUNDomEigEstimator)(farg1);
+  arg2 = (N_Vector)(farg2);
+  result = (SUNErrCode)SUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(arg1,arg2);
+  fresult = (SUNErrCode)(result);
+  return fresult;
+}
+
+
 SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetIsReal_Power(SUNDomEigEstimator farg1, int const *farg2) {
   int fresult ;
   SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ;
diff --git a/src/sundomeigest/power/fmod_int64/fsundomeigest_power_mod.f90 b/src/sundomeigest/power/fmod_int64/fsundomeigest_power_mod.f90
index 41b78edd03..5c38720a8f 100644
--- a/src/sundomeigest/power/fmod_int64/fsundomeigest_power_mod.f90
+++ b/src/sundomeigest/power/fmod_int64/fsundomeigest_power_mod.f90
@@ -46,10 +46,10 @@ module fsundomeigest_power_mod
   procedure :: get_ATdata => swigf_SUNDomEigEstimatorContent_Power__ATdata_get
   procedure :: set_V => swigf_SUNDomEigEstimatorContent_Power__V_set
   procedure :: get_V => swigf_SUNDomEigEstimatorContent_Power__V_get
-  procedure :: set_q => swigf_SUNDomEigEstimatorContent_Power__q_set
-  procedure :: get_q => swigf_SUNDomEigEstimatorContent_Power__q_get
-  procedure :: set_q_prev => swigf_SUNDomEigEstimatorContent_Power__q_prev_set
-  procedure :: get_q_prev => swigf_SUNDomEigEstimatorContent_Power__q_prev_get
+  procedure :: set_Av => swigf_SUNDomEigEstimatorContent_Power__Av_set
+  procedure :: get_Av => swigf_SUNDomEigEstimatorContent_Power__Av_get
+  procedure :: set_v_prev => swigf_SUNDomEigEstimatorContent_Power__v_prev_set
+  procedure :: get_v_prev => swigf_SUNDomEigEstimatorContent_Power__v_prev_get
   procedure :: set_rhs_linY => swigf_SUNDomEigEstimatorContent_Power__rhs_linY_set
   procedure :: get_rhs_linY => swigf_SUNDomEigEstimatorContent_Power__rhs_linY_get
   procedure :: set_Fy => swigf_SUNDomEigEstimatorContent_Power__Fy_set
@@ -78,6 +78,8 @@ module fsundomeigest_power_mod
   procedure :: get_nfevals => swigf_SUNDomEigEstimatorContent_Power__nfevals_get
   procedure :: set_is_complex => swigf_SUNDomEigEstimatorContent_Power__is_complex_set
   procedure :: get_is_complex => swigf_SUNDomEigEstimatorContent_Power__is_complex_get
+  procedure :: set_Fy_is_current => swigf_SUNDomEigEstimatorContent_Power__Fy_is_current_set
+  procedure :: get_Fy_is_current => swigf_SUNDomEigEstimatorContent_Power__Fy_is_current_get
   procedure :: release => swigf_release_SUNDomEigEstimatorContent_Power_
   procedure, private :: swigf_SUNDomEigEstimatorContent_Power__op_assign__
   generic :: assignment(=) => swigf_SUNDomEigEstimatorContent_Power__op_assign__
@@ -93,6 +95,7 @@ module fsundomeigest_power_mod
  public :: FSUNDomEigEstimator_SetRelTol_Power
  public :: FSUNDomEigEstimator_SetInitialGuess_Power
  public :: FSUNDomEigEstimator_SetRhsLinearizationPoint_Power
+ public :: FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power
  public :: FSUNDomEigEstimator_SetIsReal_Power
  public :: FSUNDomEigEstimator_Initialize_Power
  public :: FSUNDomEigEstimator_Estimate_Power
@@ -156,16 +159,16 @@ function swigc_SUNDomEigEstimatorContent_Power__V_get(farg1) &
 type(C_PTR) :: fresult
 end function
 
-subroutine swigc_SUNDomEigEstimatorContent_Power__q_set(farg1, farg2) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__q_set")
+subroutine swigc_SUNDomEigEstimatorContent_Power__Av_set(farg1, farg2) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__Av_set")
 use, intrinsic :: ISO_C_BINDING
 import :: swigclasswrapper
 type(SwigClassWrapper) :: farg1
 type(C_PTR), value :: farg2
 end subroutine
 
-function swigc_SUNDomEigEstimatorContent_Power__q_get(farg1) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__q_get") &
+function swigc_SUNDomEigEstimatorContent_Power__Av_get(farg1) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__Av_get") &
 result(fresult)
 use, intrinsic :: ISO_C_BINDING
 import :: swigclasswrapper
@@ -173,16 +176,16 @@ function swigc_SUNDomEigEstimatorContent_Power__q_get(farg1) &
 type(C_PTR) :: fresult
 end function
 
-subroutine swigc_SUNDomEigEstimatorContent_Power__q_prev_set(farg1, farg2) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__q_prev_set")
+subroutine swigc_SUNDomEigEstimatorContent_Power__v_prev_set(farg1, farg2) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__v_prev_set")
 use, intrinsic :: ISO_C_BINDING
 import :: swigclasswrapper
 type(SwigClassWrapper) :: farg1
 type(C_PTR), value :: farg2
 end subroutine
 
-function swigc_SUNDomEigEstimatorContent_Power__q_prev_get(farg1) &
-bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__q_prev_get") &
+function swigc_SUNDomEigEstimatorContent_Power__v_prev_get(farg1) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__v_prev_get") &
 result(fresult)
 use, intrinsic :: ISO_C_BINDING
 import :: swigclasswrapper
@@ -428,6 +431,23 @@ function swigc_SUNDomEigEstimatorContent_Power__is_complex_get(farg1) &
 integer(C_INT) :: fresult
 end function
 
+subroutine swigc_SUNDomEigEstimatorContent_Power__Fy_is_current_set(farg1, farg2) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__Fy_is_current_set")
+use, intrinsic :: ISO_C_BINDING
+import :: swigclasswrapper
+type(SwigClassWrapper) :: farg1
+integer(C_INT), intent(in) :: farg2
+end subroutine
+
+function swigc_SUNDomEigEstimatorContent_Power__Fy_is_current_get(farg1) &
+bind(C, name="_wrap_SUNDomEigEstimatorContent_Power__Fy_is_current_get") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+import :: swigclasswrapper
+type(SwigClassWrapper) :: farg1
+integer(C_INT) :: fresult
+end function
+
 function swigc_new_SUNDomEigEstimatorContent_Power_() &
 bind(C, name="_wrap_new_SUNDomEigEstimatorContent_Power_") &
 result(fresult)
@@ -528,6 +548,15 @@ function swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint_Power(farg1, farg2,
 integer(C_INT) :: fresult
 end function
 
+function swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(farg1, farg2) &
+bind(C, name="_wrap_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+type(C_PTR), value :: farg1
+type(C_PTR), value :: farg2
+integer(C_INT) :: fresult
+end function
+
 function swigc_FSUNDomEigEstimator_SetIsReal_Power(farg1, farg2) &
 bind(C, name="_wrap_FSUNDomEigEstimator_SetIsReal_Power") &
 result(fresult)
@@ -688,19 +717,19 @@ function swigf_SUNDomEigEstimatorContent_Power__V_get(self) &
 call c_f_pointer(fresult, swig_result)
 end function
 
-subroutine swigf_SUNDomEigEstimatorContent_Power__q_set(self, q)
+subroutine swigf_SUNDomEigEstimatorContent_Power__Av_set(self, av)
 use, intrinsic :: ISO_C_BINDING
 class(SUNDomEigEstimatorContent_Power_), intent(in) :: self
-type(N_Vector), target, intent(inout) :: q
+type(N_Vector), target, intent(inout) :: av
 type(SwigClassWrapper) :: farg1 
 type(C_PTR) :: farg2 
 
 farg1 = self%swigdata
-farg2 = c_loc(q)
-call swigc_SUNDomEigEstimatorContent_Power__q_set(farg1, farg2)
+farg2 = c_loc(av)
+call swigc_SUNDomEigEstimatorContent_Power__Av_set(farg1, farg2)
 end subroutine
 
-function swigf_SUNDomEigEstimatorContent_Power__q_get(self) &
+function swigf_SUNDomEigEstimatorContent_Power__Av_get(self) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
 type(N_Vector), pointer :: swig_result
@@ -709,23 +738,23 @@ function swigf_SUNDomEigEstimatorContent_Power__q_get(self) &
 type(SwigClassWrapper) :: farg1 
 
 farg1 = self%swigdata
-fresult = swigc_SUNDomEigEstimatorContent_Power__q_get(farg1)
+fresult = swigc_SUNDomEigEstimatorContent_Power__Av_get(farg1)
 call c_f_pointer(fresult, swig_result)
 end function
 
-subroutine swigf_SUNDomEigEstimatorContent_Power__q_prev_set(self, q_prev)
+subroutine swigf_SUNDomEigEstimatorContent_Power__v_prev_set(self, v_prev)
 use, intrinsic :: ISO_C_BINDING
 class(SUNDomEigEstimatorContent_Power_), intent(in) :: self
-type(N_Vector), target, intent(inout) :: q_prev
+type(N_Vector), target, intent(inout) :: v_prev
 type(SwigClassWrapper) :: farg1 
 type(C_PTR) :: farg2 
 
 farg1 = self%swigdata
-farg2 = c_loc(q_prev)
-call swigc_SUNDomEigEstimatorContent_Power__q_prev_set(farg1, farg2)
+farg2 = c_loc(v_prev)
+call swigc_SUNDomEigEstimatorContent_Power__v_prev_set(farg1, farg2)
 end subroutine
 
-function swigf_SUNDomEigEstimatorContent_Power__q_prev_get(self) &
+function swigf_SUNDomEigEstimatorContent_Power__v_prev_get(self) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
 type(N_Vector), pointer :: swig_result
@@ -734,7 +763,7 @@ function swigf_SUNDomEigEstimatorContent_Power__q_prev_get(self) &
 type(SwigClassWrapper) :: farg1 
 
 farg1 = self%swigdata
-fresult = swigc_SUNDomEigEstimatorContent_Power__q_prev_get(farg1)
+fresult = swigc_SUNDomEigEstimatorContent_Power__v_prev_get(farg1)
 call c_f_pointer(fresult, swig_result)
 end function
 
@@ -1088,6 +1117,31 @@ function swigf_SUNDomEigEstimatorContent_Power__is_complex_get(self) &
 swig_result = fresult
 end function
 
+subroutine swigf_SUNDomEigEstimatorContent_Power__Fy_is_current_set(self, fy_is_current)
+use, intrinsic :: ISO_C_BINDING
+class(SUNDomEigEstimatorContent_Power_), intent(in) :: self
+integer(C_INT), intent(in) :: fy_is_current
+type(SwigClassWrapper) :: farg1 
+integer(C_INT) :: farg2 
+
+farg1 = self%swigdata
+farg2 = fy_is_current
+call swigc_SUNDomEigEstimatorContent_Power__Fy_is_current_set(farg1, farg2)
+end subroutine
+
+function swigf_SUNDomEigEstimatorContent_Power__Fy_is_current_get(self) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT) :: swig_result
+class(SUNDomEigEstimatorContent_Power_), intent(in) :: self
+integer(C_INT) :: fresult 
+type(SwigClassWrapper) :: farg1 
+
+farg1 = self%swigdata
+fresult = swigc_SUNDomEigEstimatorContent_Power__Fy_is_current_get(farg1)
+swig_result = fresult
+end function
+
 function swigf_create_SUNDomEigEstimatorContent_Power_() &
 result(self)
 use, intrinsic :: ISO_C_BINDING
@@ -1268,6 +1322,22 @@ function FSUNDomEigEstimator_SetRhsLinearizationPoint_Power(dee, t, v) &
 swig_result = fresult
 end function
 
+function FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(dee, fyt) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+integer(C_INT) :: swig_result
+type(SUNDomEigEstimator), target, intent(inout) :: dee
+type(N_Vector), target, intent(inout) :: fyt
+integer(C_INT) :: fresult 
+type(C_PTR) :: farg1 
+type(C_PTR) :: farg2 
+
+farg1 = c_loc(dee)
+farg2 = c_loc(fyt)
+fresult = swigc_FSUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(farg1, farg2)
+swig_result = fresult
+end function
+
 function FSUNDomEigEstimator_SetIsReal_Power(dee, real) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
diff --git a/src/sundomeigest/power/sundomeigest_power.c b/src/sundomeigest/power/sundomeigest_power.c
index 7b8b5eb905..3062645378 100644
--- a/src/sundomeigest/power/sundomeigest_power.c
+++ b/src/sundomeigest/power/sundomeigest_power.c
@@ -111,6 +111,8 @@ SUNDomEigEstimator SUNDomEigEstimator_Power(N_Vector q, long int max_iters,
   DEE->ops->setrhs    = SUNDomEigEstimator_SetRhs_Power;
   DEE->ops->setrhslinearizationpoint =
     SUNDomEigEstimator_SetRhsLinearizationPoint_Power;
+  DEE->ops->setrhsatlinearizationpoint =
+    SUNDomEigEstimator_SetRhsAtLinearizationPoint_Power;
   DEE->ops->setmaxiters = SUNDomEigEstimator_SetMaxIters_Power;
   DEE->ops->setnumpreprocessiters = SUNDomEigEstimator_SetNumPreprocessIters_Power;
   DEE->ops->setreltol         = SUNDomEigEstimator_SetRelTol_Power;
@@ -133,28 +135,32 @@ SUNDomEigEstimator SUNDomEigEstimator_Power(N_Vector q, long int max_iters,
   DEE->content = content;
 
   /* Fill content */
-  content->ATimes      = NULL;
-  content->ATdata      = NULL;
-  content->V           = NULL;
-  content->q           = NULL;
-  content->q_prev      = NULL;
-  content->rhs_linY    = NULL;
-  content->rhs_linT    = ZERO;
-  content->Fy          = NULL;
-  content->work        = NULL;
-  content->is_complex  = SUNTRUE;
-  content->max_iters   = max_iters;
-  content->num_warmups = DEE_NUM_OF_WARMUPS_PI_DEFAULT;
-  content->rel_tol     = rel_tol;
-  content->res         = ZERO;
-  content->rhsfn       = NULL;
-  content->rhs_data    = NULL;
-  content->nfevals     = 0;
-  content->num_iters   = 0;
-  content->num_ATimes  = 0;
+  content->ATimes        = NULL;
+  content->ATdata        = NULL;
+  content->V             = NULL;
+  content->Av            = NULL;
+  content->v_prev        = NULL;
+  content->rhs_linY      = NULL;
+  content->rhs_linT      = ZERO;
+  content->Fy            = NULL;
+  content->work          = NULL;
+  content->is_complex    = SUNTRUE;
+  content->Fy_is_current = SUNFALSE;
+  content->max_iters     = max_iters;
+  content->num_warmups   = DEE_NUM_OF_WARMUPS_PI_DEFAULT;
+  content->rel_tol       = rel_tol;
+  content->res           = ZERO;
+  content->rhsfn         = NULL;
+  content->rhs_data      = NULL;
+  content->nfevals       = 0;
+  content->num_iters     = 0;
+  content->num_ATimes    = 0;
 
   /* Allocate content */
-  content->q = N_VClone(q);
+  content->Av = N_VClone(q);
+  SUNCheckLastErrNull();
+
+  content->v_prev = N_VClone(q);
   SUNCheckLastErrNull();
 
   content->V = N_VClone(q);
@@ -237,6 +243,30 @@ SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint_Power(
   return SUN_SUCCESS;
 }
 
+SUNErrCode SUNDomEigEstimator_SetRhsAtLinearizationPoint_Power(
+  SUNDomEigEstimator DEE, N_Vector Fyt)
+{
+  SUNFunctionBegin(DEE->sunctx);
+
+  SUNAssert(DEE, SUN_ERR_ARG_CORRUPT);
+  SUNAssert(PI_CONTENT(DEE), SUN_ERR_ARG_CORRUPT);
+  SUNAssert(Fyt, SUN_ERR_ARG_CORRUPT);
+
+  if (PI_CONTENT(DEE)->Fy == NULL)
+  {
+    PI_CONTENT(DEE)->Fy = N_VClone(Fyt);
+    SUNCheckLastErr();
+  }
+
+  N_VScale(ONE, Fyt, PI_CONTENT(DEE)->Fy);
+  SUNCheckLastErr();
+
+  /* Set the flag to indicate that Fy is current */
+  PI_CONTENT(DEE)->Fy_is_current = SUNTRUE;
+
+  return SUN_SUCCESS;
+}
+
 SUNErrCode SUNDomEigEstimator_SetIsReal_Power(SUNDomEigEstimator DEE,
                                               sunbooleantype real)
 {
@@ -248,13 +278,20 @@ SUNErrCode SUNDomEigEstimator_SetIsReal_Power(SUNDomEigEstimator DEE,
   /* set the complex flag to the opposite of the real flag */
   PI_CONTENT(DEE)->is_complex = !real;
 
-  /* q_prev is allocated in SUNDomEigEstimator_Initialize_Power, which is expected to be 
-  called after this routine. If the user calls this routine after initialization, we need 
-  to free q_prev here. */
-  if (!(PI_CONTENT(DEE)->is_complex) && PI_CONTENT(DEE)->q_prev)
+  /* v_prev is allocated in the constructor. If the user calls this routine later, 
+  we need to free v_prev here. */
+  if (real && PI_CONTENT(DEE)->v_prev != NULL)
+  {
+    N_VDestroy(PI_CONTENT(DEE)->v_prev);
+    PI_CONTENT(DEE)->v_prev = NULL;
+  }
+
+  /* If the user requests complex eigenvalues, we need to allocate v_prev. */
+  if (!real && PI_CONTENT(DEE)->v_prev == NULL)
   {
-    N_VDestroy(PI_CONTENT(DEE)->q_prev);
-    PI_CONTENT(DEE)->q_prev = NULL;
+    /* allocate v_prev vector */
+    PI_CONTENT(DEE)->v_prev = N_VClone(PI_CONTENT(DEE)->Av);
+    SUNCheckLastErr();
   }
 
   return SUN_SUCCESS;
@@ -283,11 +320,11 @@ SUNErrCode SUNDomEigEstimator_Initialize_Power(SUNDomEigEstimator DEE)
 
   SUNAssert(PI_CONTENT(DEE)->ATimes, SUN_ERR_ARG_CORRUPT);
   SUNAssert(PI_CONTENT(DEE)->V, SUN_ERR_ARG_CORRUPT);
-  SUNAssert(PI_CONTENT(DEE)->q, SUN_ERR_ARG_CORRUPT);
+  SUNAssert(PI_CONTENT(DEE)->Av, SUN_ERR_ARG_CORRUPT);
 
   if (PI_CONTENT(DEE)->is_complex)
   {
-    SUNAssert(PI_CONTENT(DEE)->q_prev == NULL, SUN_ERR_ARG_CORRUPT);
+    SUNAssert(PI_CONTENT(DEE)->v_prev, SUN_ERR_ARG_CORRUPT);
   }
 
   return SUN_SUCCESS;
@@ -377,70 +414,80 @@ SUNErrCode SUNDomEigEstimator_Estimate_Power(SUNDomEigEstimator DEE,
   SUNAssert(lambdaI, SUN_ERR_ARG_CORRUPT);
   SUNAssert(PI_CONTENT(DEE)->ATimes, SUN_ERR_ARG_CORRUPT);
   SUNAssert(PI_CONTENT(DEE)->V, SUN_ERR_ARG_CORRUPT);
-  SUNAssert(PI_CONTENT(DEE)->q, SUN_ERR_ARG_CORRUPT);
   SUNAssert((PI_CONTENT(DEE)->max_iters >= 0), SUN_ERR_ARG_CORRUPT);
-  if (PI_CONTENT(DEE)->is_complex && (PI_CONTENT(DEE)->q_prev == NULL))
-  {
-    /* allocate q_prev vector */
-    PI_CONTENT(DEE)->q_prev = N_VClone(PI_CONTENT(DEE)->q);
-    SUNCheckLastErr();
-  }
+
+  sunbooleantype is_complex = PI_CONTENT(DEE)->is_complex;
+
+  if (is_complex) { SUNAssert(PI_CONTENT(DEE)->v_prev, SUN_ERR_ARG_CORRUPT); }
+
+  N_Vector Av     = PI_CONTENT(DEE)->Av;
+  N_Vector V      = PI_CONTENT(DEE)->V;
+  N_Vector v_prev = PI_CONTENT(DEE)->v_prev;
+
+  int num_warmups = PI_CONTENT(DEE)->num_warmups;
+
+  long int* num_ATimes = &(PI_CONTENT(DEE)->num_ATimes);
+  long int* num_iters  = &(PI_CONTENT(DEE)->num_iters);
+  long int max_iters   = PI_CONTENT(DEE)->max_iters;
+
+  sunrealtype rel_tol = PI_CONTENT(DEE)->rel_tol;
+  sunrealtype* res    = &(PI_CONTENT(DEE)->res);
+
+  void* ATdata = PI_CONTENT(DEE)->ATdata;
+
+  SUNATimesFn ATimes = PI_CONTENT(DEE)->ATimes;
 
   sunrealtype newlambdaR = ZERO;
   sunrealtype oldlambdaR = ZERO;
 
   int retval;
   sunbooleantype converged;
-  sunrealtype normq;
-  PI_CONTENT(DEE)->num_ATimes = 0;
-  PI_CONTENT(DEE)->num_iters  = 0;
+  sunrealtype normw;
+  *num_ATimes = 0;
+  *num_iters  = 0;
 
   /* Set the initial q = A^{num_warmups}q/||A^{num_warmups}q|| */
-  for (int i = 0; i < PI_CONTENT(DEE)->num_warmups; i++)
+  for (int i = 0; i < num_warmups; i++)
   {
-    retval = PI_CONTENT(DEE)->ATimes(PI_CONTENT(DEE)->ATdata,
-                                     PI_CONTENT(DEE)->V, PI_CONTENT(DEE)->q);
-    PI_CONTENT(DEE)->num_ATimes++;
-    PI_CONTENT(DEE)->num_iters++;
+    retval = ATimes(ATdata, V, Av);
+    (*num_ATimes)++;
+    (*num_iters)++;
     if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; }
 
-    normq = N_VDotProd(PI_CONTENT(DEE)->q, PI_CONTENT(DEE)->q);
+    normw = N_VDotProd(Av, Av);
     SUNCheckLastErr();
 
-    normq = SUNRsqrt(normq);
-    N_VScale(ONE / normq, PI_CONTENT(DEE)->q, PI_CONTENT(DEE)->V);
+    normw = SUNRsqrt(normw);
+    N_VScale(ONE / normw, Av, V);
     SUNCheckLastErr();
   }
 
-  for (int k = 0; k < PI_CONTENT(DEE)->max_iters; k++)
+  for (int k = 0; k < max_iters; k++)
   {
-    if (PI_CONTENT(DEE)->is_complex)
+    if (is_complex)
     {
-      N_VScale(ONE, PI_CONTENT(DEE)->V, PI_CONTENT(DEE)->q_prev);
+      N_VScale(ONE, V, v_prev);
       SUNCheckLastErr();
     }
 
-    retval = PI_CONTENT(DEE)->ATimes(PI_CONTENT(DEE)->ATdata,
-                                     PI_CONTENT(DEE)->V, PI_CONTENT(DEE)->q);
-    PI_CONTENT(DEE)->num_ATimes++;
-    PI_CONTENT(DEE)->num_iters++;
+    retval = ATimes(ATdata, V, Av);
+    (*num_ATimes)++;
+    (*num_iters)++;
     if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; }
 
-    newlambdaR = N_VDotProd(PI_CONTENT(DEE)->V,
-                            PI_CONTENT(DEE)->q); //Rayleigh quotient
+    newlambdaR = N_VDotProd(V, Av); //Rayleigh quotient
     SUNCheckLastErr();
 
-    PI_CONTENT(DEE)->res = SUNRabs(newlambdaR - oldlambdaR);
-    converged =
-      (PI_CONTENT(DEE)->res <= PI_CONTENT(DEE)->rel_tol * SUNRabs(newlambdaR));
+    *res      = SUNRabs(newlambdaR - oldlambdaR);
+    converged = (*res <= rel_tol * SUNRabs(newlambdaR));
 
-    if (converged && !PI_CONTENT(DEE)->is_complex) { break; }
+    if (converged && !is_complex) { break; }
 
-    normq = N_VDotProd(PI_CONTENT(DEE)->q, PI_CONTENT(DEE)->q);
+    normw = N_VDotProd(Av, Av);
     SUNCheckLastErr();
 
-    normq = SUNRsqrt(normq);
-    N_VScale(ONE / normq, PI_CONTENT(DEE)->q, PI_CONTENT(DEE)->V);
+    normw = SUNRsqrt(normw);
+    N_VScale(ONE / normw, Av, V);
     SUNCheckLastErr();
 
     if (converged) { break; }
@@ -448,12 +495,11 @@ SUNErrCode SUNDomEigEstimator_Estimate_Power(SUNDomEigEstimator DEE,
     oldlambdaR = newlambdaR;
   }
 
-  if (PI_CONTENT(DEE)->is_complex)
+  if (is_complex)
   {
-    retval = sundomeigestimator_complex_dom_eigs_from_PI(DEE, newlambdaR, normq,
-                                                         PI_CONTENT(DEE)->q_prev,
-                                                         PI_CONTENT(DEE)->V,
-                                                         lambdaR, lambdaI);
+    retval = sundomeigestimator_complex_dom_eigs_from_PI(DEE, newlambdaR, normw,
+                                                         v_prev, V, lambdaR,
+                                                         lambdaI);
     if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; }
   }
   else
@@ -555,11 +601,17 @@ SUNErrCode sundomeigestimator_complex_dom_eigs_from_PI(
   SUNFunctionBegin(DEE->sunctx);
 
   int retval;
+  N_Vector Av          = PI_CONTENT(DEE)->Av;
+  SUNATimesFn ATimes   = PI_CONTENT(DEE)->ATimes;
+  void* ATdata         = PI_CONTENT(DEE)->ATdata;
+  sunrealtype rel_tol  = PI_CONTENT(DEE)->rel_tol;
+  long int* num_ATimes = &(PI_CONTENT(DEE)->num_ATimes);
+
   sunrealtype cos_qs, gram_det, det_G_inv, h11, h12, h22, p11, p12, p21, p22;
   /* The threshold for identifying real or complex DEE is experimentally 
-  determined based on the relative tolerance PI_CONTENT(DEE)->rel_tol */
+  determined based on the relative tolerance rel_tol */
   sunrealtype gram_det_tol = SUNMAX(SUN_RCONST(10.0) * SUN_UNIT_ROUNDOFF,
-                                    SUN_RCONST(10.0) * PI_CONTENT(DEE)->rel_tol);
+                                    SUN_RCONST(10.0) * rel_tol);
   cos_qs                   = N_VDotProd(v_prev, v);
   SUNCheckLastErr();
 
@@ -588,14 +640,13 @@ SUNErrCode sundomeigestimator_complex_dom_eigs_from_PI(
 
     h11 = lambdaR;
 
-    retval = PI_CONTENT(DEE)->ATimes(PI_CONTENT(DEE)->ATdata, v,
-                                     PI_CONTENT(DEE)->q);
-    PI_CONTENT(DEE)->num_ATimes++;
+    retval = ATimes(ATdata, v, Av);
+    (*num_ATimes)++;
     if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; }
 
-    h12 = N_VDotProd(v_prev, PI_CONTENT(DEE)->q);
+    h12 = N_VDotProd(v_prev, Av);
     SUNCheckLastErr();
-    h22 = N_VDotProd(v, PI_CONTENT(DEE)->q);
+    h22 = N_VDotProd(v, Av);
     SUNCheckLastErr();
 
     p11 = det_G_inv * (h11 - cos_qs * h21);
@@ -640,15 +691,15 @@ SUNErrCode SUNDomEigEstimator_Destroy_Power(SUNDomEigEstimator* DEEptr)
   if (DEE->content)
   {
     /* delete items from within the content structure */
-    if (PI_CONTENT(DEE)->q)
+    if (PI_CONTENT(DEE)->Av)
     {
-      N_VDestroy(PI_CONTENT(DEE)->q);
-      PI_CONTENT(DEE)->q = NULL;
+      N_VDestroy(PI_CONTENT(DEE)->Av);
+      PI_CONTENT(DEE)->Av = NULL;
     }
-    if (PI_CONTENT(DEE)->q_prev)
+    if (PI_CONTENT(DEE)->v_prev)
     {
-      N_VDestroy(PI_CONTENT(DEE)->q_prev);
-      PI_CONTENT(DEE)->q_prev = NULL;
+      N_VDestroy(PI_CONTENT(DEE)->v_prev);
+      PI_CONTENT(DEE)->v_prev = NULL;
     }
     if (PI_CONTENT(DEE)->V)
     {
@@ -726,14 +777,26 @@ SUNErrCode dee_DQJtimes_Power(void* voidstarDEE, N_Vector v, N_Vector Jv)
   sunrealtype sig, siginv;
   int iter, retval;
 
-  N_Vector y    = PI_CONTENT(DEE)->rhs_linY;
-  N_Vector work = PI_CONTENT(DEE)->work;
-  N_Vector Fy   = PI_CONTENT(DEE)->Fy;
+  sunbooleantype* Fy_is_current = &(PI_CONTENT(DEE)->Fy_is_current);
+  N_Vector y                    = PI_CONTENT(DEE)->rhs_linY;
+  N_Vector work                 = PI_CONTENT(DEE)->work;
+  N_Vector Fy                   = PI_CONTENT(DEE)->Fy;
+
+  SUNRhsFn rhsfn = PI_CONTENT(DEE)->rhsfn;
+  void* rhs_data = PI_CONTENT(DEE)->rhs_data;
 
-  retval = PI_CONTENT(DEE)->rhsfn(PI_CONTENT(DEE)->rhs_linT, y, Fy,
-                                  PI_CONTENT(DEE)->rhs_data);
-  PI_CONTENT(DEE)->nfevals++;
-  if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; }
+  sunrealtype rhs_linT = PI_CONTENT(DEE)->rhs_linT;
+
+  long int* nfevals = &(PI_CONTENT(DEE)->nfevals);
+
+  if (!(*Fy_is_current))
+  {
+    /* If Fy is not current, compute it at the linearization point */
+    retval = rhsfn(rhs_linT, y, Fy, rhs_data);
+    (*nfevals)++;
+    if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; }
+    *Fy_is_current = SUNTRUE;
+  }
 
   /* Initialize perturbation */
   sunrealtype ydotv   = N_VDotProd(y, v);
@@ -748,9 +811,8 @@ SUNErrCode dee_DQJtimes_Power(void* voidstarDEE, N_Vector v, N_Vector Jv)
     N_VLinearSum(sig, v, ONE, y, work);
 
     /* Set Jv = f(tn, y+sig*v) */
-    retval = PI_CONTENT(DEE)->rhsfn(PI_CONTENT(DEE)->rhs_linT, work, Jv,
-                                    PI_CONTENT(DEE)->rhs_data);
-    PI_CONTENT(DEE)->nfevals++;
+    retval = rhsfn(rhs_linT, work, Jv, rhs_data);
+    (*nfevals)++;
     if (retval == 0) { break; }
     if (retval < 0) { return SUN_ERR_USER_FCN_FAIL; }
 
@@ -765,5 +827,6 @@ SUNErrCode dee_DQJtimes_Power(void* voidstarDEE, N_Vector v, N_Vector Jv)
   siginv = ONE / sig;
   N_VLinearSum(siginv, Jv, -siginv, Fy, Jv);
 
+  *Fy_is_current = SUNFALSE;
   return SUN_SUCCESS;
 }