Feature: SUNDomEigEstimator upgrades

.github/workflows/ci.yml

@@ -96,7 +96,8 @@ jobs:
           # Check for changes to trigger python interface checks
           # Updates in include or bindings/sundials4py directories
           if git diff --name-only "$BASE" "$HEAD" | grep -qE '^(include/|bindings/sundials4py/)'; then
-            echo "python=true" >> "$GITHUB_OUTPUT"
+            #echo "python=true" >> "$GITHUB_OUTPUT"
+            echo "python=false" >> "$GITHUB_OUTPUT"
           else
             echo "python=false" >> "$GITHUB_OUTPUT"
           fi

doc/arkode/guide/source/Usage/LSRKStep/User_callable.rst

@@ -315,6 +315,43 @@ Allowable Method Families
       This routine will be called by :c:func:`ARKodeSetOptions`
       when using the key "arkid.dom_eig_safety_factor".
 
+.. c:function:: int LSRKStepSetUseAnalyticStabRegion(void* arkode_mem, sunbooleantype analytic_stab_region);
+
+   Specifies whether to use the analytic stability region for determining the number of stages in STS methods.
+
+   **Arguments:**
+      * *arkode_mem* -- pointer to the LSRKStep memory block.
+      * *analytic_stab_region* -- Use the analytic stability region if ``SUNTRUE``; use the inscribed ellipse stability region if ``SUNFALSE``.
+
+   **Return value:**
+      * *ARK_SUCCESS* if successful
+      * *ARK_MEM_NULL* if ``arkode_mem`` was ``NULL``.
+
+   .. note::
+
+      If :c:func:`LSRKStepSetUseAnalyticStabRegion` is not called, then the default
+      ``analytic_stab_region`` is set to ``SUNFALSE``.  This routine will be called by
+      :c:func:`ARKodeSetOptions` when using the key "arkid.use_analytic_stab_region".
+
+      :c:func:`LSRKStepSetUseAnalyticStabRegion` sets whether to use the ellipse or the exact stability region for 
+      stability checks. The stability region check for RKC and RKL methods is performed with the dominant
+      eigenvalue and the current step size to ensure stability. While this is sufficient for a stability
+      region of disk, the stability region of RKC and RKL methods is not a disk but rather a more complicated 
+      region that can be approximated by an inscribed ellipse. By default, the ellipse is used for stability 
+      checks, which is a conservative approximation of the stability region that possibly reduces the step 
+      sizes to ensure stability. Setting use_analytic_stab_region to SUNTRUE allows the use of the exact 
+      stability region, which can potentially allow for larger step sizes but possibly cause stability failures 
+      for the second dominant eigenvalue since it might be outside of the stability region even if the dominant 
+      eigenvalue is inside the stability region. Using ellipse for stability checks can be beneficial when two 
+      dominant eigenvalues are close to the stability boundary. Nevertheless, unless the full spectrum is used 
+      for stability checks, there is always a risk of stability failures one way or another. Thus, users 
+      have the option to choose their preference for a more conservative or more aggressive 
+      approach to stability. This input is only used for RKC and RKL methods.
+
+      .. figure:: ../../../../../arkode/guide/source/figs/arkode/STS2_region_s10.png
+         :alt: Stability region of RKL method with 10 stages
+         :align: center
+         :width: 50%
 
 .. c:function:: int LSRKStepSetNumDomEigEstInitPreprocessIters(void* arkode_mem, int num_iters);
 

doc/shared/sundomeigest/SUNDomEigEst_API.rst

@@ -259,6 +259,26 @@ instead of supplying a dummy routine.
 
    .. note::
 
+      The relative tolerance is used as a stopping criterion for the Power Iteration method. Specifically, 
+      it defines the acceptable relative change between successive dominant eigenvalue estimates.
+
+      When used in combination with Arnoldi Iteration, this routine sets the stopping criterion for 
+      the preprocessing Power Iteration phase. In this workflow, Power Iteration is first applied to 
+      perform an initial convergence check for the dominant eigenvalue estimate. The convergence
+      criterion is based on the change in magnitude between successive estimates and is defined as
+
+      .. math::
+
+         \left|\,|\lambda_{k}| - |\lambda_{k-1}|\,\right|
+         \le \texttt{rel_tol} \cdot |\lambda_{k}|.
+
+      The implementation performs this inexpensive preprocessing check using only the magnitude of 
+      the eigenvalue estimates. Arnoldi Iteration is executed only after this convergence criterion is 
+      satisfied. This approach ensures that Arnoldi is performed only once, rather than repeatedly 
+      running Arnoldi and checking convergence based on its complex-valued results. By relying on 
+      the cheaper magnitude-based preprocessing step, the routine avoids multiple Arnoldi runs that would
+      yield only marginal improvements in accuracy while incurring significantly higher computational cost.
+
       This routine will be called by :c:func:`SUNDomEigEstimator_SetOptions`
       when using the key "Did.rel_tol".
 

doc/shared/sundomeigest/SUNDomEigEst_Introduction.rst

@@ -38,15 +38,54 @@ Moreover, advanced users can provide a customized :c:type:`SUNDomEigEstimator`
 implementation to any SUNDIALS package, particularly in cases where they
 provide their own :c:type:`N_Vector`.
 
-While Krylov-based estimators preset the number of Krylov subspace
-dimensions, resulting in a tolerance-free estimation, SUNDIALS requires
-that iterative estimators stop when the residual meets a prescribed
-tolerance, :math:`\tau`,
+The default Power iteration implementation estimates complex-valued dominant
+eigenvalues. After the iteration phase, a postprocessing step is performed
+using the two most recent iterate vectors (approximations of the dominant
+eigenvector). These vectors are used to construct a 2×2 projection of the
+original matrix.
+
+If the two iterates are (numerically) linearly dependent, this indicates
+convergence to a one-dimensional invariant subspace, consistent with a
+real-valued dominant eigenvalue. In this case, the dominant eigenvalue
+estimate is taken as the Rayleigh quotient of the final iterate.
+
+If the iterates are not linearly dependent, they span a two-dimensional
+subspace. A 2×2 projection of the original matrix onto this subspace is
+constructed, and the eigenvalues of this projected matrix are used as the
+dominant eigenvalue estimates. This allows the method to capture complex
+conjugate dominant eigenvalue pairs.
+
+Convergence is determined using a magnitude-based relative tolerance
+criterion. Given successive eigenvalue estimates :math:`\lambda_k` and
+:math:`\lambda_{k-1}`, convergence is achieved when
 
 .. math::
-  :name: pi_rel_tol
+   :name: pi_rel_tol
+
+   \frac{\left||\lambda_k| - |\lambda_{k-1}|\right|}{\left|\lambda_k \right|} < \tau.
+
+where :math:`\tau` is a prescribed tolerance.
+
+An option is also provided to estimate only a real-valued dominant
+eigenvalue. In this mode, the 2×2 projection step is skipped and the
+Rayleigh quotient of the final iterate is returned directly.
+
+Krylov-based estimators use a fixed number of Krylov subspace dimensions in
+order to guarantee a bounded memory footprint, which is essential for
+large-scale problems. Unlike the Power iteration, these implementations do
+not perform tolerance-based convergence checks at every Arnoldi step since
+repeating an Arnoldi iteration due to failed convergence would be
+computationally expensive.
+
+Instead, the magnitude-based convergence criterion defined in
+:ref:`relative tolerance <pi_rel_tol>` is used as a preliminary screening
+mechanism before invoking the Krylov-based estimator. While this approach
+is slightly less robust than explicitly monitoring both the real and
+imaginary components of the eigenvalue residual during Arnoldi iteration,
+it significantly reduces computational cost. This trade-off is particularly
+advantageous for large-scale problems, where each Arnoldi cycle may be
+expensive.
 
-  \frac{\left|\lambda_k - \lambda_{k-1}\right|}{\left|\lambda_k \right|} < \tau.
 
 For users interested in providing their own :c:func:`SUNDomEigEstimator`, the
 following section presents the :c:type:`SUNDomEigEstimator` class and its implementation

doc/shared/sundomeigest/SUNDomEigEst_Power.rst

@@ -109,6 +109,27 @@ routines:
       :c:func:`SUNDomEigEstimator_SetInitialGuess`.
 
 
+.. c:function:: SUNErrCode SUNDomEigEstimator_SetDEEisReal_Power(SUNDomEigEstimator DEE, sunbooleantype real);
+
+   This routine sets the ``real`` flag in the Power iteration dominant eigenvalue estimator object.
+
+   :param DEE: the dominant eigenvalue estimator object.
+   :param real: the value to set the ``real`` flag to.
+
+   :returns: If successful, ``SUN_SUCCESS`` otherwise an appropriate error code.
+
+   .. note::
+
+      If ``real`` is ``SUNTRUE``, then the convergence criterion is based on the relative change 
+      in the magnitude of successive eigenvalue estimates.  If ``real`` is ``SUNFALSE``, then a 
+      postprocessing step is performed to compute the complex-valued dominant eigenvalue estimate, 
+      and the convergence criterion is based on the relative change in the magnitude of successive 
+      eigenvalue estimates. :c:func:`SUNDomEigEstimator_SetRelTol` can be used to set 
+      the relative tolerance for this convergence criterion.
+
+      The default value is ``SUNFALSE``.
+
+
 .. _SUNDomEigEst.Power.Description:
 
 SUNDomEigEstimator_Power Description
@@ -195,6 +216,8 @@ eigenvalue estimator operations listed in :numref:`SUNDomEigEst.API`:
 
 * ``SUNDomEigEstimator_SetRelTol_Power``
 
+* ``SUNDomEigEstimator_SetDEEisReal_Power``
+
 * ``SUNDomEigEstimator_Initialize_Power``
 
 * ``SUNDomEigEstimator_Estimate_Power``

examples/arkode/C_serial/ark_analytic_lsrk.out

@@ -19,22 +19,22 @@ Analytical ODE test problem:
    ---------------------
 
 Final Statistics:
-Current time                  = 10.0040469322476
-Steps                         = 1120
-Step attempts                 = 1124
+Current time                  = 10.0013387372969
+Steps                         = 1122
+Step attempts                 = 1125
 Stability limited steps       = 232
-Accuracy limited steps        = 1124
-Error test fails              = 4
+Accuracy limited steps        = 1125
+Error test fails              = 3
 NLS step fails                = 0
 Inequality constraint fails   = 0
 Initial step size             = 1.93010111094261e-10
 Last step size                = 0.01791
-Current step size             = 0.0370361796227445
-RHS fn evals                  = 141974
+Current step size             = 0.0393007765090647
+RHS fn evals                  = 141846
 Number of dom_eig updates     = 1
 Max. num. of stages used      = 199
 Max. num. of stages allowed   = 200
 Max. spectral radius          = 1010000
 Min. spectral radius          = 1010000
 
-ACCURACY at the final time   = 1.52101e-13
+ACCURACY at the final time   = 1.04139e-13

examples/arkode/C_serial/ark_analytic_lsrk_domeigest.c

@@ -237,6 +237,10 @@ int main(int argc, char* argv[])
   flag = ARKodeSetOptions(arkode_mem, NULL, NULL, argc, argv);
   if (check_flag(&flag, "ARKodeSetOptions", 1)) { return 1; }
 
+  /* Set real type dominant eigenvalue */
+  flag = SUNDomEigEstimator_SetDEEisReal_Power(DEE, SUNTRUE);
+  if (check_flag(&flag, "SUNDomEigEstimator_SetDEEisReal_Power", 1)) { return 1; }
+
   /* Open output stream for results, output comment line */
   UFID = fopen("solution.txt", "w");
   fprintf(UFID, "# t u\n");

examples/arkode/C_serial/ark_analytic_lsrk_domeigest.out

@@ -23,24 +23,24 @@ The stiffness of the problem is directly proportional to
    ---------------------
 
 Final Statistics:
-Current time                  = 10.0225297852833
-Steps                         = 1228
-Step attempts                 = 1293
-Stability limited steps       = 100
-Accuracy limited steps        = 1293
-Error test fails              = 65
+Current time                  = 10.0071779045287
+Steps                         = 1223
+Step attempts                 = 1280
+Stability limited steps       = 53
+Accuracy limited steps        = 1280
+Error test fails              = 57
 NLS step fails                = 0
 Inequality constraint fails   = 0
 Initial step size             = 1.93010111094261e-10
-Last step size                = 0.0231428574002256
-Current step size             = 0.0355362892083706
-RHS fn evals                  = 134017
-Number of dom_eig updates     = 85
-Number of fe calls for DEE    = 304
-Number of iterations for DEE  = 180
+Last step size                = 0.0232720439260913
+Current step size             = 0.026140590649044
+RHS fn evals                  = 132154
+Number of dom_eig updates     = 81
+Number of fe calls for DEE    = 278
+Number of iterations for DEE  = 172
 Max. num. of stages used      = 199
 Max. num. of stages allowed   = 200
-Max. spectral radius          = 1010100.99886782
+Max. spectral radius          = 1010100.56057129
 Min. spectral radius          = 1009899.00000001
 
-ACCURACY at the final time   = 2.41362e-13
+ACCURACY at the final time   = 5.97078e-13

examples/arkode/C_serial/ark_analytic_lsrk_domeigest_arkid.dom_eig_est_init_preprocess_iters_1_sundomeigestimator.max_iters_1.out

@@ -23,24 +23,24 @@ The stiffness of the problem is directly proportional to
    ---------------------
 
 Final Statistics:
-Current time                  = 10.0078124292941
-Steps                         = 1220
-Step attempts                 = 1279
-Stability limited steps       = 109
-Accuracy limited steps        = 1279
-Error test fails              = 59
+Current time                  = 10.0012842870746
+Steps                         = 1225
+Step attempts                 = 1282
+Stability limited steps       = 44
+Accuracy limited steps        = 1282
+Error test fails              = 57
 NLS step fails                = 0
 Inequality constraint fails   = 0
 Initial step size             = 1.93010111094261e-10
-Last step size                = 0.0231428888979688
-Current step size             = 0.0371525856132128
-RHS fn evals                  = 132879
+Last step size                = 0.0232720481291628
+Current step size             = 0.0257097908695329
+RHS fn evals                  = 132323
 Number of dom_eig updates     = 81
-Number of fe calls for DEE    = 147
+Number of fe calls for DEE    = 144
 Number of iterations for DEE  = 91
 Max. num. of stages used      = 199
 Max. num. of stages allowed   = 200
-Max. spectral radius          = 1010099.62410851
+Max. spectral radius          = 1010100.37814108
 Min. spectral radius          = 1009899
 
-ACCURACY at the final time   = 1.98952e-13
+ACCURACY at the final time   = 7.18314e-13

examples/arkode/C_serial/ark_analytic_lsrk_varjac.out

@@ -23,22 +23,22 @@ The stiffness of the problem is directly proportional to
    ---------------------
 
 Final Statistics:
-Current time                  = 10.0176850761672
-Steps                         = 1218
-Step attempts                 = 1274
-Stability limited steps       = 102
-Accuracy limited steps        = 1274
-Error test fails              = 56
+Current time                  = 10.0024706965917
+Steps                         = 1225
+Step attempts                 = 1283
+Stability limited steps       = 47
+Accuracy limited steps        = 1283
+Error test fails              = 58
 NLS step fails                = 0
 Inequality constraint fails   = 0
 Initial step size             = 1.93010111094261e-10
-Last step size                = 0.0231428715572232
-Current step size             = 0.035071188930862
-RHS fn evals                  = 133006
-Number of dom_eig updates     = 80
+Last step size                = 0.0232720622166514
+Current step size             = 0.0258318595494996
+RHS fn evals                  = 132509
+Number of dom_eig updates     = 81
 Max. num. of stages used      = 199
 Max. num. of stages allowed   = 200
-Max. spectral radius          = 1010100.38096708
+Max. spectral radius          = 1010099.76668785
 Min. spectral radius          = 1009899
 
-ACCURACY at the final time   = 2.42473e-13
+ACCURACY at the final time   = 6.7879e-13