Feature: Dominant Eigenvalue Estimator (#710)

Add dominant eigenvalue estimator class

---------

Co-authored-by: Daniel R. Reynolds <reynolds@smu.edu>
Co-authored-by: Daniel R. Reynolds <dreynolds@umbc.edu>
Co-authored-by: Balos, Cody, J <balos1@llnl.gov>
Co-authored-by: David J. Gardner <gardner48@llnl.gov>

Author: 4 people

CHANGELOG.md

@@ -1,9 +1,19 @@
 # SUNDIALS Changelog
 
+
 ## Changes to SUNDIALS in release X.Y.Z
 
 ### Major Features
 
+Added the `SUNDomEigEstimator` interface for estimating the dominant value of a system. 
+Two implementations are provided: Power Iteration and Arnoldi Iteration. The latter
+method requires building with LAPACK support enabled.
+
+Added the function `LSRKStepSetDomEigEstimator` in LSRKStep to attach a
+`SUNDomEigEstimator`, when using Runge-Kutta-Chebyshev or Runge-Kutta-Legendre
+methods, as an alternative to supplying a user-defined function to compute the dominant
+eigenvalue. 
+
 ### New Features and Enhancements
 
 The functions `KINSetMAA` and `KINSetOrthAA` have been updated to allow for
@@ -36,6 +46,7 @@ erroneous forcing terms.
 
 ### Deprecation Notices
 
+
 ## Changes to SUNDIALS in release 7.4.0
 
 ### New Features and Enhancements

CMakeLists.txt

@@ -102,6 +102,9 @@ set(sunlinsollib_SOVERSION "5")
 set(sunnonlinsollib_VERSION "4.4.0")
 set(sunnonlinsollib_SOVERSION "4")
 
+set(sundomeigestlib_VERSION "1.0.0")
+set(sundomeigestlib_SOVERSION "1")
+
 set(sundialslib_VERSION
     "${PACKAGE_VERSION_MAJOR}.${PACKAGE_VERSION_MINOR}.${PACKAGE_VERSION_PATCH}"
 )

cmake/SUNDIALSConfig.cmake.in

@@ -72,8 +72,8 @@ if("@ENABLE_ADIAK@" AND NOT TARGET adiak::adiak)
   find_dependency(adiak PATHS "@adiak_DIR@")
 endif()
 
-if("@ENABLE_CUDA@" AND NOT (TARGET CUDA::cudart AND TARGET CUDA::cublas
-   AND TARGET CUDA::cusparse AND TARGET CUDA::cusolver))
+if("@ENABLE_CUDA@" AND NOT (TARGET CUDA::cudart AND TARGET CUDA::curand
+   AND TARGET CUDA::cublas AND TARGET CUDA::cusparse AND TARGET CUDA::cusolver))
   find_dependency(CUDAToolkit)
 endif()
 

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

@@ -182,23 +182,70 @@ Allowable Method Families
 
 .. c:function:: int LSRKStepSetDomEigFn(void* arkode_mem, ARKDomEigFn dom_eig);
 
-   Specifies the dominant eigenvalue approximation routine to
+   Specifies the user-supplied dominant eigenvalue approximation routine to
    be used for determining the number of stages that will be used by either the
-   RKC or RKL methods.
+   RKC or RKL methods.  ``dom_eig = NULL`` input creates internal SUNDIALS Dominant
+   Eigenvalue Estimator (DEE) to serve this goal.
 
    **Arguments:**
       * *arkode_mem* -- pointer to the LSRKStep memory block.
       * *dom_eig* -- name of user-supplied dominant eigenvalue approximation function (of type :c:func:`ARKDomEigFn()`).
 
    **Return value:**
       * *ARK_SUCCESS* if successful
-      * *ARKLS_MEM_NULL* if ``arkode_mem`` was ``NULL``.
-      * *ARK_ILL_INPUT* ``dom_eig = NULL`` and LSRKStep does not currently estimate this internally.
+      * *ARK_MEM_NULL* if ``arkode_mem`` was ``NULL``.
+      * *ARK_DEE_FAIL* if ``dom_eig = NULL`` and DEE failed.
 
    .. note::
 
       This function is currently required when either the RKC or RKL methods are used.
 
+   .. note:: *ARK_DEE_FAIL* return should also produce error messages due to DEE error(s).  These errors
+      are handled by :c:type:`SUNErrCode`.
+
+      Either this function or the DEE creator function, :c:func:`LSRKStepSetDomEigEstimator` is required
+      when either the RKC or RKL methods are used.
+
+      :c:func:`LSRKStepSetDomEigFn` expects a user-provided spectral radius function pointer of type
+      :c:func:`ARKDomEigFn()`.
+
+      Note that although a DEE creation routine requires :c:func:`SUNDomEigEst_SetATimes` with a valid
+      matrix-vector product function pointer, creating a DEE with :c:func:`LSRKStepSetDomEigFn`
+      uses an internal Jacobian-vector product estimation that is passed with the *arkode_mem* pointer.
+      Similarly, it estimates the eigenvalue as needed internally without requiring a call to
+      :c:func:`SUNDomEig_Estimate`.
+
+
+.. c:function:: int LSRKStepSetDomEigEstimator(void* arkode_mem, SUNDomEigEstimator DEE);
+
+   Specifies the user-supplied dominant eigenvalue estimator (DEE) to be used for
+   determining the number of stages that will be used by either the RKC or RKL methods.
+   This function is an alternative to :c:func:`LSRKStepSetDomEigFn` and is used when a DEE
+   has already been created. TODO: Add extra information about how we implement this function 
+   after deciding on it by the group.
+
+   **Arguments:**
+      * *arkode_mem* -- pointer to the LSRKStep memory block.
+      * *DEE* -- pointer to the SUNDIALS Dominant Eigenvalue Estimator (of type :c:type:`SUNDomEigEstimator`).
+
+   **Return value:**
+      * *ARK_SUCCESS* if successful
+      * *ARK_MEM_NULL* if ``arkode_mem`` was ``NULL``.
+      * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. DEE itself or some of the required options is ``NULL``)
+      * *ARK_DEE_FAIL* if DEE failed.
+
+   .. note:: *ARK_DEE_FAIL* return should also produce error messages due to DEE error(s).  These errors
+      are handled by :c:type:`SUNErrCode`.
+
+      Either this function or the user-supplied dominant eigenvalue estimator creator function,
+      :c:func:`LSRKStepSetDomEigFn` is required when either the RKC or RKL methods are used.
+
+      Note that although a DEE creation routine requires :c:func:`SUNDomEigEst_SetATimes` with a valid
+      matrix-vector product function pointer, setting a DEE with :c:func:`LSRKStepSetDomEigEstimator`
+      uses an internal Jacobian-vector product estimation that is passed with the *arkode_mem* pointer.
+      Similarly, it estimates the eigenvalue as needed internally without requiring a call to
+      :c:func:`SUNDomEig_Estimate`.
+
 
 .. c:function:: int LSRKStepSetDomEigFrequency(void* arkode_mem, long int nsteps);
 
@@ -211,7 +258,7 @@ Allowable Method Families
 
    **Return value:**
       * *ARK_SUCCESS* if successful
-      * *ARKLS_MEM_NULL* if ``arkode_mem`` was ``NULL``.
+      * *ARK_MEM_NULL* if ``arkode_mem`` was ``NULL``.
 
    .. note::
 
@@ -235,7 +282,7 @@ Allowable Method Families
 
    **Return value:**
       * *ARK_SUCCESS* if successful
-      * *ARKLS_MEM_NULL* if ``arkode_mem`` was ``NULL``.
+      * *ARK_MEM_NULL* if ``arkode_mem`` was ``NULL``.
 
    .. note::
 
@@ -264,7 +311,7 @@ Allowable Method Families
 
    **Return value:**
       * *ARK_SUCCESS* if successful
-      * *ARKLS_MEM_NULL* if ``arkode_mem`` was ``NULL``.
+      * *ARK_MEM_NULL* if ``arkode_mem`` was ``NULL``.
 
    .. note::
 
@@ -276,6 +323,22 @@ Allowable Method Families
       when using the key "arkid.dom_eig_safety_factor".
 
 
+.. c:function:: int LSRKStepSetNumSucceedingWarmups(void* arkode_mem, int num_succ_warmups);
+
+   Specifies the number of the preprocessing warmups before each estimate call succeeding the very first estimate call.
+
+   **Arguments:**
+      * *arkode_mem* -- pointer to the LSRKStep memory block.
+      * *num_succ_warmups* -- the number of succeeding warmups.
+
+   **Return value:**
+      * *ARK_SUCCESS* if successful
+      * *ARK_MEM_NULL* if ``arkode_mem`` was ``NULL``.
+
+.. note:: If LSRKStepSetNumSucceedingWarmups routine is not called, then the default ``num_succ_warmups`` is set to :math:`0`.
+   Calling this function with ``num_succ_warmups < 0`` resets the default.
+
+
 .. c:function:: int LSRKStepSetNumSSPStages(void* arkode_mem, int num_of_stages);
 
    Sets the number of stages, ``s`` in ``SSP(s, p)`` methods. This input is only utilized by SSPRK methods.
@@ -290,7 +353,7 @@ Allowable Method Families
 
    **Return value:**
       * *ARK_SUCCESS* if successful
-      * *ARKLS_MEM_NULL* if ``arkode_mem`` was ``NULL``.
+      * *ARK_MEM_NULL* if ``arkode_mem`` was ``NULL``.
       * *ARK_ILL_INPUT* if an argument had an illegal value (e.g. SSP method is not declared)
 
    .. note::
@@ -334,6 +397,26 @@ Optional output functions
    **Return value:**
       * *ARK_SUCCESS* if successful
       * *ARK_MEM_NULL* if the LSRKStep memory was ``NULL``
+      * *ARK_ILL_INPUT* if ``stage_max`` is illegal
+
+
+.. c:function:: int LSRKStepGetNumDomEigEstRhsEvals(void* arkode_mem, long int* nfeDQ);
+
+   Returns the number of RHS function evaluations used in the difference quotient
+   Jacobian approximations (so far).
+
+   **Arguments:**
+      * *arkode_mem* -- pointer to the LSRKStep memory block.
+      * *nfeDQ* -- number of rhs calls.
+
+   **Return value:**
+      * *ARK_SUCCESS* if successful
+      * *ARK_MEM_NULL* if the LSRKStep memory was ``NULL``
+      * *ARK_ILL_INPUT* if ``nfeDQ`` is illegal
+
+.. note:: The internal SUNDIALS dominant eigenvalue estimator (DEE) uses this approximation.
+   Therefore, it returns 0 if DEE the was not created.
+
 
 .. _ARKODE.Usage.LSRKStep.Reinitialization:
 

doc/arkode/guide/source/index.rst

@@ -64,6 +64,7 @@ with support by the `US Department of Energy <http://www.doe.gov>`_,
    sunmatrix/index.rst
    sunlinsol/index.rst
    sunnonlinsol/index.rst
+   sundomeigest/index.rst
    sunadaptcontroller/index.rst
    sunstepper/index.rst
    sunadjoint/index.rst

doc/arkode/guide/source/sundomeigest/SUNDomEigEst_links.rst

@@ -0,0 +1,17 @@
+.. ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2025, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. include:: ../../../../shared/sundomeigest/SUNDomEigEst_Introduction.rst
+.. include:: ../../../../shared/sundomeigest/SUNDomEigEst_API.rst
+.. include:: ../../../../shared/sundomeigest/SUNDomEigEst_Power.rst
+.. include:: ../../../../shared/sundomeigest/SUNDomEigEst_Arnoldi.rst
+.. include:: ../../../../shared/sundomeigest/SUNDomEigEst_Examples.rst

doc/arkode/guide/source/sundomeigest/index.rst

@@ -0,0 +1,24 @@
+.. ----------------------------------------------------------------
+   Mustafa Aggul @ SMU
+   ----------------------------------------------------------------
+   SUNDIALS Copyright Start
+   Copyright (c) 2002-2025, Lawrence Livermore National Security
+   and Southern Methodist University.
+   All rights reserved.
+
+   See the top-level LICENSE and NOTICE files for details.
+
+   SPDX-License-Identifier: BSD-3-Clause
+   SUNDIALS Copyright End
+   ----------------------------------------------------------------
+
+.. _SUNDomEigEst:
+
+##############################
+Dominant Eigenvalue Estimators
+##############################
+
+.. toctree::
+   :maxdepth: 1
+
+   SUNDomEigEst_links.rst

doc/shared/sundials.bib

@@ -2645,3 +2645,31 @@ @article{rackauckas2020universal
   journal={arXiv preprint arXiv:2001.04385},
   year={2020}
 }
+
+%
+% Arnoldi Iteration
+%
+
+@article{arnoldi51,
+  title={The principle of minimized iterations in the solution of the matrix eigenvalue problem},
+  author={Arnoldi, W. E.},
+  journal={Quarterly of Applied Mathematics},
+  volume={9},
+  number={1},
+  pages={17--29},
+  year={1951}
+}
+
+%
+% Power Iteration
+%
+
+@article{vonmises29,
+  title={Praktische Verfahren der Gleichungsauflösung},
+  author={von Mises, Richard and Pollaczek-Geiringer, Hilda},
+  journal={Zeitschrift für Angewandte Mathematik und Mechanik},
+  volume={9},
+  pages={152--164},
+  doi = {10.1002/zamm.19290090206},
+  year={1929}
+}