Feature/dom eig

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

@@ -174,20 +174,65 @@ 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);
@@ -201,7 +246,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:: If LSRKStepSetDomEigFrequency routine is not called, then the default ``nsteps`` is set to :math:`25` as recommended in :cite:p:`VSH:04`.
    Calling this function with ``nsteps < 0`` resets the default value while ``nsteps = 0`` refers to constant dominant eigenvalue.
@@ -218,7 +263,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:: If LSRKStepSetMaxNumStages routine is not called, then the default ``stage_max_limit`` is
    set to :math:`200`. Calling this function with ``stage_max_limit < 2`` resets the default value.
@@ -239,12 +284,28 @@ 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:: If LSRKStepSetDomEigSafetyFactor routine is not called, then the default ``dom_eig_safety`` is
    set to :math:`1.01`. Calling this function with ``dom_eig_safety < 1`` resets the default value.
 
 
+.. 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.
@@ -259,7 +320,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:: If LSRKStepSetNumSSPStages routine is not called, then the default ``num_of_stages`` is
@@ -298,6 +359,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}
+}