SUNDomEigEstimator Cleanup and remove extra vector memory allocation

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

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,

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

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

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
   ===============================================================*/

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

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;
   }
 

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 ;

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

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 ;

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