diff --git a/bindings/sundials4py/arkode/arkode_lsrkstep_generated.hpp b/bindings/sundials4py/arkode/arkode_lsrkstep_generated.hpp index 252f93ff92..cae989c9fe 100644 --- a/bindings/sundials4py/arkode/arkode_lsrkstep_generated.hpp +++ b/bindings/sundials4py/arkode/arkode_lsrkstep_generated.hpp @@ -42,6 +42,9 @@ m.def("LSRKStepSetMaxNumStages", LSRKStepSetMaxNumStages, nb::arg("arkode_mem"), m.def("LSRKStepSetDomEigSafetyFactor", LSRKStepSetDomEigSafetyFactor, nb::arg("arkode_mem"), nb::arg("dom_eig_safety")); +m.def("LSRKStepSetUseAnalyticStabRegion", LSRKStepSetUseAnalyticStabRegion, + nb::arg("arkode_mem"), nb::arg("analytic_stab_region")); + m.def("LSRKStepSetNumDomEigEstInitPreprocessIters", LSRKStepSetNumDomEigEstInitPreprocessIters, nb::arg("arkode_mem"), nb::arg("num_iters")); diff --git a/bindings/sundials4py/sundials/generate.yaml b/bindings/sundials4py/sundials/generate.yaml index ce6a745aac..281ae1692d 100644 --- a/bindings/sundials4py/sundials/generate.yaml +++ b/bindings/sundials4py/sundials/generate.yaml @@ -63,6 +63,7 @@ modules: - ../../include/sundials/sundials_domeigestimator.h fn_exclude_by_name__regex: - "SUNDomEigEstimator_SetATimes" # nanobind cannot bind to functions which take a function pointer, so we do something custom + - "SUNDomEigEstimator_SetRhs" # nanobind cannot bind to functions which take a function pointer, so we do something custom sundials_errors: path: sundials/sundials_errors_generated.hpp headers: diff --git a/bindings/sundials4py/sundials/sundials_domeigestimator.cpp b/bindings/sundials4py/sundials/sundials_domeigestimator.cpp index 3e77254954..ee0ac40326 100644 --- a/bindings/sundials4py/sundials/sundials_domeigestimator.cpp +++ b/bindings/sundials4py/sundials/sundials_domeigestimator.cpp @@ -83,6 +83,26 @@ void bind_sundomeigestimator(nb::module_& m) else { return SUNDomEigEstimator_SetATimes(dee, fntable, nullptr); } }, nb::arg("DEE"), nb::arg("ATimes").none()); + + m.def( + "SUNDomEigEstimator_SetRhs", + [](SUNDomEigEstimator DEE, + std::function> RHSfn) -> SUNErrCode + { + if (!DEE->python) { DEE->python = new SUNDomEigEstimatorFunctionTable; } + + auto fntable = static_cast(DEE->python); + + fntable->deerhs = nb::cast(RHSfn); + + if (RHSfn) + { + return SUNDomEigEstimator_SetRhs(DEE, fntable, + sundomeigestimator_setrhs_wrapper); + } + else { return SUNDomEigEstimator_SetRhs(DEE, fntable, nullptr); } + }, + nb::arg("DEE"), nb::arg("RHSfn").none()); } } // namespace sundials4py diff --git a/bindings/sundials4py/sundials/sundials_domeigestimator_usersupplied.hpp b/bindings/sundials4py/sundials/sundials_domeigestimator_usersupplied.hpp index d21b3f2c4c..4939125dee 100644 --- a/bindings/sundials4py/sundials/sundials_domeigestimator_usersupplied.hpp +++ b/bindings/sundials4py/sundials/sundials_domeigestimator_usersupplied.hpp @@ -33,6 +33,7 @@ using namespace sundials::experimental; struct SUNDomEigEstimatorFunctionTable { nb::object atimes; + nb::object deerhs; }; template @@ -43,4 +44,12 @@ SUNErrCode sundomeigestimator_atimes_wrapper(Args... args) 3>(&SUNDomEigEstimatorFunctionTable::atimes, std::forward(args)...); } +template +SUNErrCode sundomeigestimator_setrhs_wrapper(Args... args) +{ + return sundials4py::user_supplied_fn_caller< + std::remove_pointer_t, SUNDomEigEstimatorFunctionTable, + 1>(&SUNDomEigEstimatorFunctionTable::deerhs, std::forward(args)...); +} + #endif \ No newline at end of file diff --git a/bindings/sundials4py/sundomeigest/generate.yaml b/bindings/sundials4py/sundomeigest/generate.yaml index e82948a83a..98d6675357 100644 --- a/bindings/sundials4py/sundomeigest/generate.yaml +++ b/bindings/sundials4py/sundomeigest/generate.yaml @@ -26,11 +26,13 @@ modules: fn_exclude_by_name__regex: # Don't interface the implementation specific overrides of the generic routines - "^SUNDomEigEstimator_SetATimes_.*" + - "^SUNDomEigEstimator_SetRhs_.*" - "^SUNDomEigEstimator_SetOptions_.*" - "^SUNDomEigEstimator_SetMaxIters_.*" - "^SUNDomEigEstimator_SetNumPreprocessIters_.*" - "^SUNDomEigEstimator_SetRelTol_.*" - "^SUNDomEigEstimator_SetInitialGuess_.*" + - "^SUNDomEigEstimator_SetRhsLinearizationPoint_.*" - "^SUNDomEigEstimator_Initialize_.*" - "^SUNDomEigEstimator_Estimate_.*" - "^SUNDomEigEstimator_GetRes_.*" diff --git a/bindings/sundials4py/sundomeigest/sundomeigest_power_generated.hpp b/bindings/sundials4py/sundomeigest/sundomeigest_power_generated.hpp index 996531f676..27b613dd62 100644 --- a/bindings/sundials4py/sundomeigest/sundomeigest_power_generated.hpp +++ b/bindings/sundials4py/sundomeigest/sundomeigest_power_generated.hpp @@ -32,6 +32,9 @@ 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_SetIsReal_Power", SUNDomEigEstimator_SetIsReal_Power, + nb::arg("DEE"), nb::arg("real")); // #ifdef __cplusplus // // #endif diff --git a/doc/arkode/guide/source/Usage/LSRKStep/User_callable.rst b/doc/arkode/guide/source/Usage/LSRKStep/User_callable.rst index a006ea62a4..141651f44d 100644 --- a/doc/arkode/guide/source/Usage/LSRKStep/User_callable.rst +++ b/doc/arkode/guide/source/Usage/LSRKStep/User_callable.rst @@ -315,6 +315,47 @@ 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:: + + This input is only used for RKC and RKL methods. + + 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 exact stability region or an + ellipse that is fully inscribed in the stability region for determining stability in RKC and RKL + methods. Whichever region is selected, LSRKStep will ensure that the complex number + :math:`z=h\lambda`, where :math:`h` is the current time step size and :math:`\lambda` is the + estimated dominant eigenvalue, is in this region. Since the ellipse is smaller than the analytical + stability region it provides a more conservative estimate, which may be appropriate for problems + wherein sub-dominant eigenvalues may also limit stability of the method. Thus, the ellipse may + result in a smaller time step size, more internal stages, or both, when compared to using the + analytical stability region. We note that in either case, since LSRKStep does not examine the + full eigenvalue spectrum, there is always some risk of instability. For applications where + this is a concern, we recommend calling :c:func:`LSRKStepSetDomEigSafetyFactor` with a `dom_eig_safety` + value significantly larger than 1. + + If :c:func:`LSRKStepSetUseAnalyticStabRegion` is called during integration, the change will take effect + at the next step attempt. Both analytic and ellipse stability regions of RKC and RKL methods with 10 stages + are shown in the figure below. + + .. figure:: ../../../../../shared/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); diff --git a/doc/shared/Changelog.rst b/doc/shared/Changelog.rst index fc6a4bb95e..177b3357fa 100644 --- a/doc/shared/Changelog.rst +++ b/doc/shared/Changelog.rst @@ -29,6 +29,12 @@ Changelog Changes to SUNDIALS in release X.Y.Z ==================================== +**Bug Fixes** + +Fixed a minor bug where the number of required stages for STS methods +in the LSRKStep module was incorrectly computed using the spectral +radius instead of the real part of the Jacobian eigenvalues. + .. include:: RecentChanges_link.rst .. _Changelog.7.7.0: diff --git a/doc/shared/RecentChanges.rst b/doc/shared/RecentChanges.rst index 22e47b0213..137c809a7f 100644 --- a/doc/shared/RecentChanges.rst +++ b/doc/shared/RecentChanges.rst @@ -7,4 +7,8 @@ **Bug Fixes** +Fixed a minor bug where the number of required stages for STS methods +in the LSRKStep module was incorrectly computed using the spectral +radius instead of the real part of the Jacobian eigenvalues. + **Deprecation Notices** diff --git a/doc/shared/figs/arkode/STS2_region_s10.png b/doc/shared/figs/arkode/STS2_region_s10.png new file mode 100644 index 0000000000..47115cacbf Binary files /dev/null and b/doc/shared/figs/arkode/STS2_region_s10.png differ diff --git a/doc/shared/sundomeigest/SUNDomEigEst_API.rst b/doc/shared/sundomeigest/SUNDomEigEst_API.rst index fb9cfb030f..7465d41886 100644 --- a/doc/shared/sundomeigest/SUNDomEigEst_API.rst +++ b/doc/shared/sundomeigest/SUNDomEigEst_API.rst @@ -209,6 +209,46 @@ instead of supplying a dummy routine. A :c:type:`SUNErrCode`. +.. c:function:: SUNErrCode SUNDomEigEstimator_SetRhs(SUNDomEigEstimator DEE, void* rhs_data, SUNRhsFn RHSfn) + + For applications that do not provide a :c:type:`SUNATimesFn` function to :c:func:`SUNDomEigEstimator_SetATimes`, + the action of this matrix-vector product may be approximated internally. If the matrix corresponds with + the Jacobian of a vector-valued function, :math:`A = \frac{\partial f_{RHS}}{\partial y}(t,y)`, then the + function :math:`f_{RHS}` may optionally be input via this routine, and the Jacobian-vector products will + be approximated as + + .. math:: + Av \approx \frac{1}{\sigma}\left[ f_{RHS}(t,y+\sigma v) - f_{RHS}(t,y)\right] + + The linearization point :math:`(t,y)` should be separately supplied by calling :c:func:`SUNDomEigEstimator_SetRhsLinearizationPoint`. + + **Arguments:** + + * *DEE* -- a SUNDomEigEstimator object. + * *rhs_data* -- pointer to structure for ``RHSfn``. + * *RHSfn* -- function pointer to perform right-hand side evaluations. This is typically the same as the problem-defining function supplied to CVODE or ARKODE. + + **Return value:** + + A :c:type:`SUNErrCode`. + + +.. c:function:: SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint(SUNDomEigEstimator DEE, sunrealtype t, N_Vector y) + + This *optional* function sets the linearization point for the right-hand side function when using + :c:func:`SUNDomEigEstimator_SetRhs`. + + **Arguments:** + + * *DEE* -- a SUNDomEigEstimator object. + * *t* -- the time at which the linearization point is specified. + * *y* -- the linearization point for the right-hand side function. + + **Return value:** + + A :c:type:`SUNErrCode`. + + .. c:function:: SUNErrCode SUNDomEigEstimator_SetNumPreprocessIters(SUNDomEigEstimator DEE, int num_iters) This *optional* routine sets the number of preprocessing matrix-vector @@ -244,6 +284,7 @@ instead of supplying a dummy routine. This routine will be called by :c:func:`SUNDomEigEstimator_SetOptions` when using the key "Did.num_preprocess_iters". + .. c:function:: SUNErrCode SUNDomEigEstimator_SetRelTol(SUNDomEigEstimator DEE, sunrealtype rel_tol) This *optional* routine sets the estimator's :ref:`relative tolerance `. @@ -259,6 +300,52 @@ instead of supplying a dummy routine. .. note:: + The relative tolerance is used as a stopping criterion for either the Power Iteration method, or for the + preprocessing phase of the Arnoldi iteration. Specifically, it defines the acceptable relative change + between successive dominant eigenvalue estimates. It also serves as a threshold for determining + whether the dominant eigenvalue estimated by the Power Iteration is real- or complex-valued. + + When used to check convergence, we declare the iteration converged when the change in magnitude between successive estimates satisfies + + .. math:: + + \left|\,|\lambda_{k}| - |\lambda_{k-1}|\,\right| + \le \mathtt{rel\_tol} \cdot |\lambda_{k}|. + + When used to assess preprocessing iterations for the Arnoldi method, once the above test is satisfied, + the Arnoldi iteration commences, ensuring that Arnoldi is performed only once, as opposed to restarting + Arnoldi repeatedly. + + When ``rel_tol`` is used as a threshold to assess whether the Power Iteration dominant eigenvalue is + real versus complex, we define + + .. math:: + \mathtt{gram\_det\_tol} = 10 \cdot \max\left(\varepsilon,\; \mathtt{rel\_tol}\right) + + where :math:`\varepsilon` represents the machine precision. The ``gram_det_tol`` value is used to + assess the numerical rank of the 2×2 Gram matrix formed by the two most recent iterates + in the Power Iteration method -- if the determinant of this Gram matrix is less than or equal to + ``gram_det_tol``, the iterates are considered linearly dependent, and the dominant eigenvalue is + treated as real. For this use case, to ensure that small imaginary parts of the dominant eigenvalue are + not ignored, ``rel_tol`` should not be chosen too large. In practice, the smallest reliably detectable + imaginary part is proportional to the chosen relative tolerance, i.e., + + .. math:: + |\beta| \gtrsim \mathcal{O}(\mathtt{rel\_tol}). + + Therefore, to resolve an expected imaginary part of magnitude :math:`|\beta|`, it is recommended to choose + + .. math:: + \mathtt{rel\_tol} \ll |\beta|. + + Choosing a smaller relative tolerance improves the ability to detect weakly complex eigenvalues, + but may increase computational cost. + + Acceptable inputs to this routine include :math:`0 < \mathtt{rel\_tol} < 1 - \varepsilon`. For Power + Iteration, values outside this range result in the default value of ``rel_tol = 0.005``. For + Arnoldi, ``rel_tol < 0`` disables preprocessing, while + :math:`\mathtt{rel\_tol} \ge 1-\varepsilon` result in the default value of ``rel_tol = 0.005``. + This routine will be called by :c:func:`SUNDomEigEstimator_SetOptions` when using the key "Did.rel_tol". @@ -388,14 +475,25 @@ dominant eigenvalue estimator. *All routines are optional.* .. _SUNDomEigEst.SUNSuppliedFn: -Functions provided by SUNDIALS packages +User or SUNDIALS package provided functions --------------------------------------------- To interface with SUNDomEigEst modules, the SUNDIALS packages supply a -:c:type:`SUNATimesFn` function for evaluating the matrix-vector product. This -package-provided routine translates between the user-supplied ODE or DAE systems -and the generic dominant eigenvalue estimator API. The function types for these -routines are defined in the header file ``sundials/sundials_iterative.h``. +:c:type:`SUNATimesFn` function for evaluating the matrix-vector +product. This package-provided routine translates between the user-supplied ODE or DAE +systems and the generic dominant eigenvalue estimator API. The function type :c:type:`SUNATimesFn` is defined in the header file ``sundials/sundials_iterative.h``. + +Users who wish to use a SUNDomEigEst module in "standalone" mode, however, must provide either a :c:type:`SUNATimesFn` or a :c:type:`SUNRhsFn`, as described below. + + +.. c:type:: int (*SUNRhsFn)(sunrealtype t, N_Vector y, N_Vector ydot, void* rhs_data) + + Used to compute the right-hand side of an ODE system. This function is used + when the dominant eigenvalue estimator is configured to perform a discrete + Jacobian-vector product using quotient approximations of the Jacobian. The parameter + *rhs_data* is a pointer to any information about RHS which the function needs in order + to do its job. The time parameter :math:`t` and the vector :math:`y` should be left + unchanged. .. _SUNDomEigEst.Generic: diff --git a/doc/shared/sundomeigest/SUNDomEigEst_Arnoldi.rst b/doc/shared/sundomeigest/SUNDomEigEst_Arnoldi.rst index c12869c44c..76c16f2246 100644 --- a/doc/shared/sundomeigest/SUNDomEigEst_Arnoldi.rst +++ b/doc/shared/sundomeigest/SUNDomEigEst_Arnoldi.rst @@ -116,15 +116,24 @@ The SUNDomEigEstimator_Arnoldi module defines the *content* field of a void* ATdata; N_Vector* V; N_Vector q; + N_Vector rhs_linY; + N_Vector Fy; + N_Vector work; int kry_dim; int num_warmups; long int num_iters; + sunbooleantype warmup_to_tol; + sunrealtype tol_preprocess; + sunrealtype rhs_linT; long int num_ATimes; + SUNRhsFn rhsfn; + void* rhs_data; + long int nfevals; sunrealtype* LAPACK_A; sunrealtype* LAPACK_wr; sunrealtype* LAPACK_wi; sunrealtype* LAPACK_work; - snuindextype LAPACK_lwork; + sunindextype LAPACK_lwork; sunrealtype** LAPACK_arr; sunrealtype** Hes; }; @@ -137,7 +146,7 @@ information: * ``ATData`` - pointer to structure for ``ATimes``, -* ``V, q`` - vectors used for workspace by the Arnoldi algorithm. +* ``V, q, Fy, work`` - vectors used for workspace. * ``kry_dim`` - dimension of Krylov subspaces (default is 3), @@ -146,6 +155,21 @@ information: * ``num_iters`` - number of iterations (preprocessing and estimation) in the last :c:func:`SUNDomEigEstimator_Estimate` call, +* ``warmup_to_tol`` - enable warmup iterations (default is ``SUNFALSE``) + +* ``tol_preprocess`` - tolerance for preprocessing iterations (default is 0.005; + only used if ``warmup_to_tol`` is ``SUNTRUE``), + +* ``rhs_linY`` - state vector for linearization point, + +* ``rhs_linT`` - time value for linearization point, + +* ``rhsfn`` - user provided RHS function, + +* ``rhs_data`` - pointer to the data structure for ``rhsfn``, + +* ``nfevals`` - number of RHS evaluations, + * ``num_ATimes`` - number of calls to the ``ATimes`` function, * ``LAPACK_A, LAPACK_wr, LAPACK_wi, LAPACK_work`` - ``sunrealtype`` used for workspace by LAPACK, diff --git a/doc/shared/sundomeigest/SUNDomEigEst_Introduction.rst b/doc/shared/sundomeigest/SUNDomEigEst_Introduction.rst index 8db9db05dc..b5f9fd31a3 100644 --- a/doc/shared/sundomeigest/SUNDomEigEst_Introduction.rst +++ b/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 ` 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 diff --git a/doc/shared/sundomeigest/SUNDomEigEst_Power.rst b/doc/shared/sundomeigest/SUNDomEigEst_Power.rst index af79b90622..94d6bcca0b 100644 --- a/doc/shared/sundomeigest/SUNDomEigEst_Power.rst +++ b/doc/shared/sundomeigest/SUNDomEigEst_Power.rst @@ -109,6 +109,28 @@ routines: :c:func:`SUNDomEigEstimator_SetInitialGuess`. +.. c:function:: SUNErrCode SUNDomEigEstimator_SetIsReal_Power(SUNDomEigEstimator DEE, sunbooleantype real) + + This routine informs the Power iteration that the dominant eigenvalue is + real-valued, so that the complex projection described in Section + :numref:`SUNDomEigEst.Introduction` can be omitted. + + :param DEE: the dominant eigenvalue estimator object. + :param real: flag indicating that the dominant eigenvalue is real-valued. + + :returns: ``SUN_SUCCESS`` if successful, otherwise an appropriate error code. + + .. note:: + + No matter the value of ``real``, the convergence criterion is based on the relative change in the + magnitude of successive eigenvalue estimates (with tolerance set using + :c:func:`SUNDomEigEstimator_SetRelTol`). If ``real`` is ``SUNTRUE``, then the final Power + Iteration estimate is returned. Otherwise, a postprocessing step is performed to compute + the complex-valued dominant eigenvalue estimate. + + The default value is ``SUNFALSE``. + + .. _SUNDomEigEst.Power.Description: SUNDomEigEstimator_Power Description @@ -125,12 +147,21 @@ The SUNDomEigEstimator_Power module defines the *content* field of a void* ATdata; N_Vector* V; N_Vector q; + N_Vector q_prev; + N_Vector rhs_linY; + N_Vector Fy; + N_Vector work; int num_warmups; long int max_iters; long int num_iters; long int num_ATimes; + sunrealtype rhs_linT; sunrealtype rel_tol; sunrealtype res; + SUNRhsFn rhsfn; + void* rhs_data; + long int nfevals; + sunbooleantype complex; }; @@ -141,7 +172,7 @@ information: * ``ATData`` - pointer to structure for ``ATimes``, -* ``V, q`` - ``N_Vector`` used for workspace by the PI algorithm. +* ``V, q, q_prev, Fy, work`` - ``N_Vector`` used for workspace. * ``num_warmups`` - number of preprocessing iterations (default is 100), @@ -152,11 +183,23 @@ information: * ``num_ATimes`` - number of calls to the ``ATimes`` function, +* ``rhs_linY`` - state vector for linearization point, + +* ``rhs_linT`` - time value for linearization point, + * ``rel_tol`` - relative tolerance for the convergence criteria (default is 0.005), * ``res`` - the residual from the last :c:func:`SUNDomEigEstimator_Estimate` call. +* ``rhsfn`` - user provided RHS function, + +* ``rhs_data`` - pointer to the data structure for ``rhsfn``, + +* ``nfevals`` - number of RHS evaluations, + +* ``complex`` - flag indicating whether the dominant eigenvalue is + complex-valued (default is ``SUNTRUE``). This estimator is constructed to perform the following operations: @@ -195,6 +238,8 @@ eigenvalue estimator operations listed in :numref:`SUNDomEigEst.API`: * ``SUNDomEigEstimator_SetRelTol_Power`` +* ``SUNDomEigEstimator_SetIsReal_Power`` + * ``SUNDomEigEstimator_Initialize_Power`` * ``SUNDomEigEstimator_Estimate_Power`` diff --git a/examples/arkode/C_serial/CMakeLists.txt b/examples/arkode/C_serial/CMakeLists.txt index 1f4c7c4010..1279e06889 100644 --- a/examples/arkode/C_serial/CMakeLists.txt +++ b/examples/arkode/C_serial/CMakeLists.txt @@ -42,6 +42,8 @@ set(ARKODE_examples "ark_analytic_ssprk.c\;\;develop" "ark_brusselator_1D_mri.c\;\;develop" "ark_brusselator_fp.c\;\;exclude-single" + "ark_brusselator_lsrk_domeigest.c\;\;develop" + "ark_brusselator_lsrk_externaldomeigest.c\;\;develop" "ark_brusselator_mri.c\;\;develop" "ark_brusselator.c\;\;develop" "ark_brusselator1D_imexmri.c\;0 0.001\;exclude-single" diff --git a/examples/arkode/C_serial/ark_analytic_lsrk.out b/examples/arkode/C_serial/ark_analytic_lsrk.out index d8f00ed4ad..2441d6b625 100644 --- a/examples/arkode/C_serial/ark_analytic_lsrk.out +++ b/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 \ No newline at end of file diff --git a/examples/arkode/C_serial/ark_analytic_lsrk_domeigest.c b/examples/arkode/C_serial/ark_analytic_lsrk_domeigest.c index 5c35c8802a..6327634e11 100644 --- a/examples/arkode/C_serial/ark_analytic_lsrk_domeigest.c +++ b/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_SetIsReal_Power(DEE, SUNTRUE); + if (check_flag(&flag, "SUNDomEigEstimator_SetIsReal_Power", 1)) { return 1; } + /* Open output stream for results, output comment line */ UFID = fopen("solution.txt", "w"); fprintf(UFID, "# t u\n"); diff --git a/examples/arkode/C_serial/ark_analytic_lsrk_domeigest.out b/examples/arkode/C_serial/ark_analytic_lsrk_domeigest.out index 6518d38b4b..7bdb6371a6 100644 --- a/examples/arkode/C_serial/ark_analytic_lsrk_domeigest.out +++ b/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 \ No newline at end of file +ACCURACY at the final time = 5.97078e-13 \ No newline at end of file diff --git a/examples/arkode/C_serial/ark_analytic_lsrk_domeigest_arkid.dom_eig_est_init_preprocess_iters_1_sundomeigestimator.max_iters_1.out b/examples/arkode/C_serial/ark_analytic_lsrk_domeigest_arkid.dom_eig_est_init_preprocess_iters_1_sundomeigestimator.max_iters_1.out index c0fa551a34..17f974ce3a 100644 --- a/examples/arkode/C_serial/ark_analytic_lsrk_domeigest_arkid.dom_eig_est_init_preprocess_iters_1_sundomeigestimator.max_iters_1.out +++ b/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 diff --git a/examples/arkode/C_serial/ark_analytic_lsrk_varjac.out b/examples/arkode/C_serial/ark_analytic_lsrk_varjac.out index 5deb40defe..06ef355cc0 100644 --- a/examples/arkode/C_serial/ark_analytic_lsrk_varjac.out +++ b/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 \ No newline at end of file diff --git a/examples/arkode/C_serial/ark_brusselator_lsrk_domeigest.c b/examples/arkode/C_serial/ark_brusselator_lsrk_domeigest.c new file mode 100644 index 0000000000..26ab2d1ae0 --- /dev/null +++ b/examples/arkode/C_serial/ark_brusselator_lsrk_domeigest.c @@ -0,0 +1,373 @@ +/*----------------------------------------------------------------- + * Programmer(s): Mustafa Aggul @ SMU + * Based on + * ark_analytic_lsrk_domeigest.c by Mustafa Aggul @ SMU and + * ark_brusselator.c by Daniel R. Reynolds @ UMBC + *--------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2025-2026, Lawrence Livermore National Security, + * University of Maryland Baltimore County, and the SUNDIALS contributors. + * Copyright (c) 2013-2025, Lawrence Livermore National Security + * and Southern Methodist University. + * Copyright (c) 2002-2013, Lawrence Livermore National Security. + * All rights reserved. + * + * See the top-level LICENSE and NOTICE files for details. + * + * SPDX-License-Identifier: BSD-3-Clause + * SUNDIALS Copyright End + *--------------------------------------------------------------- + * Example problem: + * + * The following test simulates a brusselator problem from chemical + * kinetics. This is an ODE system with 3 components, Y = [u,v,w], + * satisfying the equations, + * du/dt = a - (w+1)*u + v*u^2 + * dv/dt = w*u - v*u^2 + * dw/dt = (b-w)/ep - w*u + * for t in the interval [0.0, 10.0], with initial conditions + * Y0 = [u0,v0,w0]. + * + * We have 3 different testing scenarios: + * + * Test 1: u0=3.9, v0=1.1, w0=2.8, a=1.2, b=2.5, ep=1.0e-5 + * Here, all three components exhibit a rapid transient change + * during the first 0.2 time units, followed by a slow and + * smooth evolution. + * + * Test 2: u0=1.2, v0=3.1, w0=3, a=1, b=3.5, ep=5.0e-6 + * Here, w experiences a fast initial transient, jumping 0.5 + * within a few steps. All values proceed smoothly until + * around t=6.5, when both u and v undergo a sharp transition, + * with u increaseing from around 0.5 to 5 and v decreasing + * from around 6 to 1 in less than 0.5 time units. After this + * transition, both u and v continue to evolve somewhat + * rapidly for another 1.4 time units, and finish off smoothly. + * + * Test 3: u0=3, v0=3, w0=3.5, a=0.5, b=3, ep=5.0e-4 + * Here, all components undergo very rapid initial transients + * during the first 0.3 time units, and all then proceed very + * smoothly for the remainder of the simulation. + * + * This file is hard-coded to use test 2. + * + * This program solves the problem with the LSRK method using internal + * SUNDIALS dominant eigenvalue estimation (DEE) module. + * + * 100 outputs are printed at equal intervals, and run statistics + * are printed at the end. + *-----------------------------------------------------------------*/ + +/* Header files */ +#include /* prototypes for LSRKStep fcts., consts */ +#include +#include /* serial N_Vector types, fcts., macros */ +#include +#include /* def. of SUNRsqrt, etc. */ +#include /* definition of type sunrealtype */ +#include /* access to Power Iteration module */ + +#if defined(SUNDIALS_EXTENDED_PRECISION) +#define GSYM "Lg" +#define ESYM "Le" +#define FSYM "Lf" +#else +#define GSYM "g" +#define ESYM "e" +#define FSYM "f" +#endif + +/* User-supplied Functions Called by the Solver */ +static int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); + +/* Private function to check function return values */ +static int check_flag(void* flagvalue, const char* funcname, int opt); + +/* Main Program */ +int main(int argc, char* argv[]) +{ + /* general problem parameters */ + sunrealtype T0 = SUN_RCONST(0.0); /* initial time */ + sunrealtype Tf = SUN_RCONST(10.0); /* final time */ + sunrealtype dTout = SUN_RCONST(1.0); /* time between outputs */ + sunindextype NEQ = 3; /* number of dependent vars. */ + int Nt = (int)ceil(Tf / dTout); /* number of output times */ + int test = 2; /* test problem to run */ + sunrealtype a, b, ep, u0, v0, w0; + +#if defined(SUNDIALS_DOUBLE_PRECISION) + sunrealtype reltol = SUN_RCONST(1.0e-6); /* tolerances */ + sunrealtype abstol = SUN_RCONST(1.0e-10); +#elif defined(SUNDIALS_SINGLE_PRECISION) + sunrealtype reltol = SUN_RCONST(1.0e-4); /* tolerances */ + sunrealtype abstol = SUN_RCONST(1.0e-8); +#elif defined(SUNDIALS_EXTENDED_PRECISION) + sunrealtype reltol = SUN_RCONST(1.0e-6); /* tolerances */ + sunrealtype abstol = SUN_RCONST(1.0e-10); +#endif + + /* general problem variables */ + int flag; /* reusable error-checking flag */ + N_Vector y = NULL; /* empty vector for storing solution */ + void* arkode_mem = NULL; /* empty ARKode memory structure */ + sunrealtype rdata[3]; + sunrealtype t, tout; + int iout; + + /* Dominant Eigenvalue Estimator (DEE) pointers and variables */ + SUNDomEigEstimator DEE = NULL; /* domeig estimator object */ + sunindextype max_iters = 100; /* max number of power iterations (PI)*/ + sunindextype numwarmup = 10; /* number of preprocessing warmups */ + sunrealtype rel_tol = SUN_RCONST(5.0e-3); /* relative error for PI*/ + N_Vector q = NULL; /* random initial eigenvector */ + + /* Create the SUNDIALS context object for this simulation */ + SUNContext ctx; + flag = SUNContext_Create(SUN_COMM_NULL, &ctx); + if (check_flag(&flag, "SUNContext_Create", 1)) { return 1; } + + /* set up the test problem according to the desired test */ + if (test == 1) + { + u0 = SUN_RCONST(3.9); + v0 = SUN_RCONST(1.1); + w0 = SUN_RCONST(2.8); + a = SUN_RCONST(1.2); + b = SUN_RCONST(2.5); + ep = SUN_RCONST(1.0e-5); + } + else if (test == 3) + { + u0 = SUN_RCONST(3.0); + v0 = SUN_RCONST(3.0); + w0 = SUN_RCONST(3.5); + a = SUN_RCONST(0.5); + b = SUN_RCONST(3.0); + ep = SUN_RCONST(5.0e-4); + } + else + { + u0 = SUN_RCONST(1.2); + v0 = SUN_RCONST(3.1); + w0 = SUN_RCONST(3.0); + a = SUN_RCONST(1.0); + b = SUN_RCONST(3.5); + ep = SUN_RCONST(5.0e-6); + } + + /* Initial problem output */ + printf("\nBrusselator ODE test problem:\n"); + printf(" initial conditions: u0 = %" GSYM ", v0 = %" GSYM + ", w0 = %" GSYM "\n", + u0, v0, w0); + printf(" problem parameters: a = %" GSYM ", b = %" GSYM ", ep = %" GSYM + "\n", + a, b, ep); + printf(" reltol = %.1" ESYM ", abstol = %.1" ESYM "\n\n", reltol, abstol); + + /* Initialize data structures */ + rdata[0] = a; /* set user data */ + rdata[1] = b; + rdata[2] = ep; + y = N_VNew_Serial(NEQ, ctx); /* Create serial vector for solution */ + if (check_flag((void*)y, "N_VNew_Serial", 0)) { return 1; } + NV_Ith_S(y, 0) = u0; /* Set initial conditions */ + NV_Ith_S(y, 1) = v0; + NV_Ith_S(y, 2) = w0; + + /* Call LSRKStepCreateSTS to initialize the ARK timestepper module and + specify the right-hand side function in y'=f(t,y), the initial time + T0, and the initial dependent variable vector y. */ + arkode_mem = LSRKStepCreateSTS(f, T0, y, ctx); + if (check_flag((void*)arkode_mem, "LSRKStepCreateSTS", 0)) { return 1; } + + /* Set routines */ + flag = ARKodeSetUserData(arkode_mem, + (void*)rdata); /* Pass rdata to user functions */ + if (check_flag(&flag, "ARKodeSetUserData", 1)) { return 1; } + + flag = ARKodeSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */ + if (check_flag(&flag, "ARKodeSStolerances", 1)) { return 1; } + + flag = ARKodeSetInterpolantType(arkode_mem, + ARK_INTERP_LAGRANGE); /* Specify stiff interpolant */ + if (check_flag(&flag, "ARKodeSetInterpolantType", 1)) { return 1; } + + /* Set the initial random eigenvector for the DEE */ + q = N_VClone(y); + if (check_flag(q, "N_VClone", 0)) { return 1; } + + sunrealtype* qd = N_VGetArrayPointer(q); + for (int i = 0; i < NEQ; i++) + { + qd[i] = (sunrealtype)rand() / (sunrealtype)RAND_MAX; + } + + /* Create power iteration dominant eigenvalue estimator (DEE) */ + DEE = SUNDomEigEstimator_Power(q, max_iters, rel_tol, ctx); + if (check_flag(DEE, "SUNDomEigEstimator_Power", 0)) { return 1; } + + /* After the DEE creation, random q vector is no longer needed. + It is used only to initialize the DEE */ + N_VDestroy(q); + + /* Attach the DEE to the LSRKStep module. + There is no need to set Atimes or initialize since LSRKStep provides + a default Atimes, and initialized the DEE, after it is attached. */ + flag = LSRKStepSetDomEigEstimator(arkode_mem, DEE); + if (check_flag(&flag, "LSRKStepSetDomEigEstimator", 1)) { return 1; } + + /* Set the number of preprocessing warmups. The warmup + is used to compute a "better" initial eigenvector and so an initial + eigenvalue. The warmup is performed only once by the LSRKStep module + internally unless LSRKStepSetNumDomEigEstPreprocessIters is called to set + a new number of succeeding warmups that would be executed before + every dominant eigenvalue estimation call */ + flag = LSRKStepSetNumDomEigEstInitPreprocessIters(arkode_mem, numwarmup); + if (check_flag(&flag, "LSRKStepSetNumDomEigEstInitPreprocessIters", 1)) + { + return 1; + } + + /* Specify after how many successful steps dom_eig is recomputed */ + flag = LSRKStepSetDomEigFrequency(arkode_mem, 0); + if (check_flag(&flag, "LSRKStepSetDomEigFrequency", 1)) { return 1; } + + /* Specify max number of stages allowed */ + flag = LSRKStepSetMaxNumStages(arkode_mem, 200); + if (check_flag(&flag, "LSRKStepSetMaxNumStages", 1)) { return 1; } + + /* Specify max number of steps allowed */ + flag = ARKodeSetMaxNumSteps(arkode_mem, 2000); + if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) { return 1; } + + /* Specify safety factor for user provided dom_eig */ + flag = LSRKStepSetDomEigSafetyFactor(arkode_mem, SUN_RCONST(1.01)); + if (check_flag(&flag, "LSRKStepSetDomEigSafetyFactor", 1)) { return 1; } + + /* Specify the Runge--Kutta--Chebyshev LSRK method by name */ + flag = LSRKStepSetSTSMethodByName(arkode_mem, "ARKODE_LSRK_RKL_2"); + if (check_flag(&flag, "LSRKStepSetSTSMethodByName", 1)) { return 1; } + + /* Override any current settings with command-line options */ + flag = SUNDomEigEstimator_SetOptions(DEE, NULL, NULL, argc, argv); + if (check_flag(&flag, "SUNDomEigEstimator_SetOptions", 1)) { return 1; } + + /* Override any current settings with command-line options */ + flag = ARKodeSetOptions(arkode_mem, NULL, NULL, argc, argv); + if (check_flag(&flag, "ARKodeSetOptions", 1)) { return 1; } + + /* Main time-stepping loop: calls ARKodeEvolve to perform the integration, then + prints results. Stops when the final time has been reached */ + t = T0; + tout = T0 + dTout; + printf(" t u v w\n"); + printf(" -------------------------------------------\n"); + printf(" %10.6" FSYM " %10.6" FSYM " %10.6" FSYM " %10.6" FSYM "\n", t, + NV_Ith_S(y, 0), NV_Ith_S(y, 1), NV_Ith_S(y, 2)); + + for (iout = 0; iout < Nt; iout++) + { + flag = ARKodeEvolve(arkode_mem, tout, y, &t, ARK_NORMAL); /* call integrator */ + if (check_flag(&flag, "ARKodeEvolve", 1)) { break; } + printf(" %10.6" FSYM " %10.6" FSYM " %10.6" FSYM " %10.6" FSYM + "\n", /* access/print solution */ + t, NV_Ith_S(y, 0), NV_Ith_S(y, 1), NV_Ith_S(y, 2)); + if (flag >= 0) + { /* successful solve: update time */ + tout += dTout; + tout = (tout > Tf) ? Tf : tout; + } + else + { /* unsuccessful solve: break */ + fprintf(stderr, "Solver failure, stopping integration\n"); + break; + } + } + printf(" -------------------------------------------\n"); + + /* Print final statistics */ + printf("\nFinal Statistics:\n"); + flag = ARKodePrintAllStats(arkode_mem, stdout, SUN_OUTPUTFORMAT_TABLE); + if (check_flag(&flag, "ARKodePrintAllStats", 1)) { return 1; } + + /* Clean up and return with successful completion */ + N_VDestroy(y); /* Free y vector */ + ARKodeFree(&arkode_mem); /* Free integrator memory */ + SUNDomEigEstimator_Destroy(&DEE); /* Free DEE object */ + SUNContext_Free(&ctx); /* Free context */ + + return flag; +} + +/*------------------------------- + * Functions called by the solver + *-------------------------------*/ + +/* f routine to compute the ODE RHS function f(t,y). */ +static int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) +{ + sunrealtype* rdata = (sunrealtype*)user_data; /* cast user_data to sunrealtype */ + sunrealtype a = rdata[0]; /* access data entries */ + sunrealtype b = rdata[1]; + sunrealtype ep = rdata[2]; + sunrealtype u = NV_Ith_S(y, 0); /* access solution values */ + sunrealtype v = NV_Ith_S(y, 1); + sunrealtype w = NV_Ith_S(y, 2); + + /* fill in the RHS function */ + NV_Ith_S(ydot, 0) = a - (w + 1.0) * u + v * u * u; + NV_Ith_S(ydot, 1) = w * u - v * u * u; + NV_Ith_S(ydot, 2) = (b - w) / ep - w * u; + + return 0; /* Return with success */ +} + +/*------------------------------- + * Private helper functions + *-------------------------------*/ + +/* Check function return value... + opt == 0 means SUNDIALS function allocates memory so check if + returned NULL pointer + opt == 1 means SUNDIALS function returns a flag so check if + flag >= 0 + opt == 2 means function allocates memory so check if returned + NULL pointer +*/ +static int check_flag(void* flagvalue, const char* funcname, int opt) +{ + int* errflag; + + /* Check if SUNDIALS function returned NULL pointer - no memory allocated */ + if (opt == 0 && flagvalue == NULL) + { + fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n", + funcname); + return 1; + } + + /* Check if flag < 0 */ + else if (opt == 1) + { + errflag = (int*)flagvalue; + if (*errflag < 0) + { + fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n", + funcname, *errflag); + return 1; + } + } + + /* Check if function returned NULL pointer - no memory allocated */ + else if (opt == 2 && flagvalue == NULL) + { + fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n", + funcname); + return 1; + } + + return 0; +} + +/*---- end of file ----*/ diff --git a/examples/arkode/C_serial/ark_brusselator_lsrk_domeigest.out b/examples/arkode/C_serial/ark_brusselator_lsrk_domeigest.out new file mode 100644 index 0000000000..83cb9efb40 --- /dev/null +++ b/examples/arkode/C_serial/ark_brusselator_lsrk_domeigest.out @@ -0,0 +1,41 @@ + +Brusselator ODE test problem: + initial conditions: u0 = 1.2, v0 = 3.1, w0 = 3 + problem parameters: a = 1, b = 3.5, ep = 5e-06 + reltol = 1.0e-06, abstol = 1.0e-10 + + t u v w + ------------------------------------------- + 0.000000 1.200000 3.100000 3.000000 + 1.000000 1.103852 3.013128 3.499981 + 2.000000 0.687945 3.521426 3.499988 + 3.000000 0.409445 4.277957 3.499993 + 4.000000 0.367889 4.942066 3.499994 + 5.000000 0.413862 5.510684 3.499993 + 6.000000 0.589251 5.855724 3.499990 + 7.000000 4.756506 0.735406 3.499917 + 8.000000 1.813465 1.575773 3.499968 + 9.000000 0.527922 2.807379 3.499991 + 10.000000 0.305608 3.657399 3.499995 + ------------------------------------------- + +Final Statistics: +Current time = 10.0124653670818 +Steps = 752 +Step attempts = 762 +Stability limited steps = 0 +Accuracy limited steps = 762 +Error test fails = 10 +NLS step fails = 0 +Inequality constraint fails = 0 +Initial step size = 1.13978962117636e-08 +Last step size = 0.020921379243018 +Current step size = 0.0212806907861296 +RHS fn evals = 47265 +Number of dom_eig updates = 1 +Number of fe calls for DEE = 12 +Number of iterations for DEE = 12 +Max. num. of stages used = 195 +Max. num. of stages allowed = 200 +Max. spectral radius = 202001.212020988 +Min. spectral radius = 202001.212020988 \ No newline at end of file diff --git a/examples/arkode/C_serial/ark_brusselator_lsrk_externaldomeigest.c b/examples/arkode/C_serial/ark_brusselator_lsrk_externaldomeigest.c new file mode 100644 index 0000000000..4c7a19b4d9 --- /dev/null +++ b/examples/arkode/C_serial/ark_brusselator_lsrk_externaldomeigest.c @@ -0,0 +1,387 @@ +/*----------------------------------------------------------------- + * Programmer(s): Mustafa Aggul @ SMU + * Based on + * ark_analytic_lsrk_domeigest.c by Mustafa Aggul @ SMU and + * ark_brusselator.c by Daniel R. Reynolds @ UMBC + *--------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2025-2026, Lawrence Livermore National Security, + * University of Maryland Baltimore County, and the SUNDIALS contributors. + * Copyright (c) 2013-2025, Lawrence Livermore National Security + * and Southern Methodist University. + * Copyright (c) 2002-2013, Lawrence Livermore National Security. + * All rights reserved. + * + * See the top-level LICENSE and NOTICE files for details. + * + * SPDX-License-Identifier: BSD-3-Clause + * SUNDIALS Copyright End + *--------------------------------------------------------------- + * Example problem: + * + * The following test simulates the same problem as + * ark_brusselator_lsrk_domeigest.c. + * + * The only difference with that example is that this one + * attaches a "user provided" dominant eigenvalue estimate, + * and where that estimate is provided by a user-constructed + * SUNDomEigEstimator object (see the dom_eig function below). + * + * While this example and ark_brusselator_lsrk_domeigest.c should provide + * the same results, the purpose of this example is to show + * users how they may manually construct SUNDomEigEstimator + * objects to query the dominant eigenvalues of a desired + * function. In particular, we note that there is no + * requirement for the SUNDomEigEstimator to be used purely + * for super-time-stepping methods in LSRKStep, so users + * interested in other purposes may focus on the contents + * of the dom_eig function below as a template. + *-----------------------------------------------------------------*/ + +/* Header files */ +#include /* prototypes for LSRKStep fcts., consts */ +#include +#include /* serial N_Vector types, fcts., macros */ +#include +#include /* def. of SUNRsqrt, etc. */ +#include /* definition of type sunrealtype */ +#include /* access to Power Iteration module */ + +#if defined(SUNDIALS_EXTENDED_PRECISION) +#define GSYM "Lg" +#define ESYM "Le" +#define FSYM "Lf" +#else +#define GSYM "g" +#define ESYM "e" +#define FSYM "f" +#endif + +/* User-supplied Functions Called by the Solver */ +static int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data); + +/* User-supplied Dominated Eigenvalue Called by the Solver */ +static int dom_eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR, + sunrealtype* lambdaI, void* user_data, N_Vector temp1, + N_Vector temp2, N_Vector temp3); + +/* Private function to check function return values */ +static int check_flag(void* flagvalue, const char* funcname, int opt); + +/* user data structure */ +typedef struct +{ + SUNContext ctx; + sunrealtype rdata[3]; + SUNDomEigEstimator DEE; + sunrealtype rel_tol; + sunindextype max_iters; +} UserData; + +/* Main Program */ +int main(int argc, char* argv[]) +{ + /* general problem parameters */ + sunrealtype T0 = SUN_RCONST(0.0); /* initial time */ + sunrealtype Tf = SUN_RCONST(10.0); /* final time */ + sunrealtype dTout = SUN_RCONST(1.0); /* time between outputs */ + sunindextype NEQ = 3; /* number of dependent vars. */ + int Nt = (int)ceil(Tf / dTout); /* number of output times */ + int test = 2; /* test problem to run */ + sunrealtype a, b, ep, u0, v0, w0; + +#if defined(SUNDIALS_DOUBLE_PRECISION) + sunrealtype reltol = SUN_RCONST(1.0e-6); /* tolerances */ + sunrealtype abstol = SUN_RCONST(1.0e-10); +#elif defined(SUNDIALS_SINGLE_PRECISION) + sunrealtype reltol = SUN_RCONST(1.0e-4); /* tolerances */ + sunrealtype abstol = SUN_RCONST(1.0e-8); +#elif defined(SUNDIALS_EXTENDED_PRECISION) + sunrealtype reltol = SUN_RCONST(1.0e-6); /* tolerances */ + sunrealtype abstol = SUN_RCONST(1.0e-10); +#endif + + /* general problem variables */ + int flag; /* reusable error-checking flag */ + N_Vector y = NULL; /* empty vector for storing solution */ + void* arkode_mem = NULL; /* empty ARKode memory structure */ + UserData ProbData; /* problem data structure */ + sunrealtype t, tout; + int iout; + + /* Create the SUNDIALS context object for this simulation */ + SUNContext ctx; + flag = SUNContext_Create(SUN_COMM_NULL, &ctx); + if (check_flag(&flag, "SUNContext_Create", 1)) { return 1; } + + /* set up the test problem according to the desired test */ + if (test == 1) + { + u0 = SUN_RCONST(3.9); + v0 = SUN_RCONST(1.1); + w0 = SUN_RCONST(2.8); + a = SUN_RCONST(1.2); + b = SUN_RCONST(2.5); + ep = SUN_RCONST(1.0e-5); + } + else if (test == 3) + { + u0 = SUN_RCONST(3.0); + v0 = SUN_RCONST(3.0); + w0 = SUN_RCONST(3.5); + a = SUN_RCONST(0.5); + b = SUN_RCONST(3.0); + ep = SUN_RCONST(5.0e-4); + } + else + { + u0 = SUN_RCONST(1.2); + v0 = SUN_RCONST(3.1); + w0 = SUN_RCONST(3.0); + a = SUN_RCONST(1.0); + b = SUN_RCONST(3.5); + ep = SUN_RCONST(5.0e-6); + } + + /* Initial problem output */ + printf("\nBrusselator ODE test problem:\n"); + printf(" initial conditions: u0 = %" GSYM ", v0 = %" GSYM + ", w0 = %" GSYM "\n", + u0, v0, w0); + printf(" problem parameters: a = %" GSYM ", b = %" GSYM ", ep = %" GSYM + "\n", + a, b, ep); + printf(" reltol = %.1" ESYM ", abstol = %.1" ESYM "\n\n", reltol, abstol); + + /* Initialize data structures */ + ProbData.ctx = ctx; + ProbData.rdata[0] = a; /* set user data */ + ProbData.rdata[1] = b; + ProbData.rdata[2] = ep; + ProbData.DEE = NULL; + ProbData.rel_tol = SUN_RCONST(5.0e-3); + ProbData.max_iters = 100; + y = N_VNew_Serial(NEQ, ctx); /* Create serial vector for solution */ + if (check_flag((void*)y, "N_VNew_Serial", 0)) { return 1; } + NV_Ith_S(y, 0) = u0; /* Set initial conditions */ + NV_Ith_S(y, 1) = v0; + NV_Ith_S(y, 2) = w0; + + /* Call LSRKStepCreateSTS to initialize the ARK timestepper module and + specify the right-hand side function in y'=f(t,y), the initial time + T0, and the initial dependent variable vector y. */ + arkode_mem = LSRKStepCreateSTS(f, T0, y, ctx); + if (check_flag((void*)arkode_mem, "LSRKStepCreateSTS", 0)) { return 1; } + + /* Set routines */ + flag = ARKodeSetUserData(arkode_mem, + (void*)&ProbData); /* Pass rdata to user functions */ + if (check_flag(&flag, "ARKodeSetUserData", 1)) { return 1; } + + flag = ARKodeSStolerances(arkode_mem, reltol, abstol); /* Specify tolerances */ + if (check_flag(&flag, "ARKodeSStolerances", 1)) { return 1; } + + flag = ARKodeSetInterpolantType(arkode_mem, + ARK_INTERP_LAGRANGE); /* Specify stiff interpolant */ + if (check_flag(&flag, "ARKodeSetInterpolantType", 1)) { return 1; } + + /* Specify user provided spectral radius */ + flag = LSRKStepSetDomEigFn(arkode_mem, dom_eig); + if (check_flag(&flag, "LSRKStepSetDomEigFn", 1)) { return 1; } + + /* Specify after how many successful steps dom_eig is recomputed */ + flag = LSRKStepSetDomEigFrequency(arkode_mem, 0); + if (check_flag(&flag, "LSRKStepSetDomEigFrequency", 1)) { return 1; } + + /* Specify max number of stages allowed */ + flag = LSRKStepSetMaxNumStages(arkode_mem, 200); + if (check_flag(&flag, "LSRKStepSetMaxNumStages", 1)) { return 1; } + + /* Specify max number of steps allowed */ + flag = ARKodeSetMaxNumSteps(arkode_mem, 2000); + if (check_flag(&flag, "ARKodeSetMaxNumSteps", 1)) { return 1; } + + /* Specify safety factor for user provided dom_eig */ + flag = LSRKStepSetDomEigSafetyFactor(arkode_mem, SUN_RCONST(1.01)); + if (check_flag(&flag, "LSRKStepSetDomEigSafetyFactor", 1)) { return 1; } + + /* Specify the Runge--Kutta--Chebyshev LSRK method by name */ + flag = LSRKStepSetSTSMethodByName(arkode_mem, "ARKODE_LSRK_RKL_2"); + if (check_flag(&flag, "LSRKStepSetSTSMethodByName", 1)) { return 1; } + + /* Override any current settings with command-line options */ + flag = ARKodeSetOptions(arkode_mem, NULL, NULL, argc, argv); + if (check_flag(&flag, "ARKodeSetOptions", 1)) { return 1; } + + /* Main time-stepping loop: calls ARKodeEvolve to perform the integration, then + prints results. Stops when the final time has been reached */ + t = T0; + tout = T0 + dTout; + printf(" t u v w\n"); + printf(" -------------------------------------------\n"); + printf(" %10.6" FSYM " %10.6" FSYM " %10.6" FSYM " %10.6" FSYM "\n", t, + NV_Ith_S(y, 0), NV_Ith_S(y, 1), NV_Ith_S(y, 2)); + + for (iout = 0; iout < Nt; iout++) + { + flag = ARKodeEvolve(arkode_mem, tout, y, &t, ARK_NORMAL); /* call integrator */ + if (check_flag(&flag, "ARKodeEvolve", 1)) { break; } + printf(" %10.6" FSYM " %10.6" FSYM " %10.6" FSYM " %10.6" FSYM + "\n", /* access/print solution */ + t, NV_Ith_S(y, 0), NV_Ith_S(y, 1), NV_Ith_S(y, 2)); + if (flag >= 0) + { /* successful solve: update time */ + tout += dTout; + tout = (tout > Tf) ? Tf : tout; + } + else + { /* unsuccessful solve: break */ + fprintf(stderr, "Solver failure, stopping integration\n"); + break; + } + } + printf(" -------------------------------------------\n"); + + /* Print final statistics */ + printf("\nFinal Statistics:\n"); + flag = ARKodePrintAllStats(arkode_mem, stdout, SUN_OUTPUTFORMAT_TABLE); + if (check_flag(&flag, "ARKodePrintAllStats", 1)) { return 1; } + + /* Clean up and return with successful completion */ + N_VDestroy(y); /* Free y vector */ + ARKodeFree(&arkode_mem); /* Free integrator memory */ + SUNDomEigEstimator_Destroy(&ProbData.DEE); /* Free DEE object */ + SUNContext_Free(&ctx); /* Free context */ + + return flag; +} + +/*------------------------------- + * Functions called by the solver + *-------------------------------*/ + +/* f routine to compute the ODE RHS function f(t,y). */ +static int f(sunrealtype t, N_Vector y, N_Vector ydot, void* user_data) +{ + UserData* data = (UserData*)user_data; /* cast user_data to UserData */ + sunrealtype* rdata = data->rdata; /* access rdata from UserData */ + sunrealtype a = rdata[0]; /* access data entries */ + sunrealtype b = rdata[1]; + sunrealtype ep = rdata[2]; + sunrealtype u = NV_Ith_S(y, 0); /* access solution values */ + sunrealtype v = NV_Ith_S(y, 1); + sunrealtype w = NV_Ith_S(y, 2); + + /* fill in the RHS function */ + NV_Ith_S(ydot, 0) = a - (w + 1.0) * u + v * u * u; + NV_Ith_S(ydot, 1) = w * u - v * u * u; + NV_Ith_S(ydot, 2) = (b - w) / ep - w * u; + + return 0; /* Return with success */ +} + +/* dom_eig routine to estimate the dominated eigenvalue */ +static int dom_eig(sunrealtype t, N_Vector y, N_Vector fn, sunrealtype* lambdaR, + sunrealtype* lambdaI, void* user_data, N_Vector temp1, + N_Vector temp2, N_Vector temp3) +{ + int flag; + UserData* data = (UserData*)user_data; /* cast user_data to UserData */ + + SUNContext ctx = data->ctx; /* access context from UserData */ + SUNDomEigEstimator DEE = data->DEE; /* access DEE from UserData */ + + /* DEE is initialized to NULL, so on the first dom_eig call we need + to create and initialize this object */ + if (DEE == NULL) + { + /* Create random initial vector for power iteration */ + sunrealtype* qd = N_VGetArrayPointer(temp1); + sunindextype NEQ = N_VGetLength(temp1); + for (int i = 0; i < NEQ; i++) + { + qd[i] = (sunrealtype)rand() / (sunrealtype)RAND_MAX; + } + + /* Create power iteration dominant eigenvalue estimator (DEE) */ + DEE = SUNDomEigEstimator_Power(temp1, data->max_iters, data->rel_tol, ctx); + if (check_flag(DEE, "SUNDomEigEstimator_Power", 0)) { return 1; } + data->DEE = DEE; + + /* Set the ODE right-hand side function at t for the Jacobian-vector products */ + flag = SUNDomEigEstimator_SetRhs(DEE, user_data, f); + if (check_flag(&flag, "SUNDomEigEstimator_SetRhs", 1)) { return 1; } + + /* Set the linearization vector for the Jacobian-vector products */ + flag = SUNDomEigEstimator_SetRhsLinearizationPoint(DEE, t, y); + if (check_flag(&flag, "SUNDomEigEstimator_SetRhsLinearizationPoint", 1)) + { + return 1; + } + + flag = SUNDomEigEstimator_Initialize(DEE); + if (check_flag(&flag, "SUNDomEigEstimator_Initialize", 1)) { return 1; } + } + + /* Update the linearization vector and time for the Jacobian-vector products */ + flag = SUNDomEigEstimator_SetRhsLinearizationPoint(DEE, t, y); + if (check_flag(&flag, "SUNDomEigEstimator_SetRhsLinearizationPoint", 1)) + { + return 1; + } + + /* Estimate the dominant eigenvalue with power iteration */ + flag = SUNDomEigEstimator_Estimate(DEE, lambdaR, lambdaI); + if (check_flag(&flag, "SUNDomEigEstimator_Estimate", 1)) { return 1; } + + return 0; /* return with success */ +} + +/*------------------------------- + * Private helper functions + *-------------------------------*/ + +/* Check function return value... + opt == 0 means SUNDIALS function allocates memory so check if + returned NULL pointer + opt == 1 means SUNDIALS function returns a flag so check if + flag >= 0 + opt == 2 means function allocates memory so check if returned + NULL pointer +*/ +static int check_flag(void* flagvalue, const char* funcname, int opt) +{ + int* errflag; + + /* Check if SUNDIALS function returned NULL pointer - no memory allocated */ + if (opt == 0 && flagvalue == NULL) + { + fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed - returned NULL pointer\n\n", + funcname); + return 1; + } + + /* Check if flag < 0 */ + else if (opt == 1) + { + errflag = (int*)flagvalue; + if (*errflag < 0) + { + fprintf(stderr, "\nSUNDIALS_ERROR: %s() failed with flag = %d\n\n", + funcname, *errflag); + return 1; + } + } + + /* Check if function returned NULL pointer - no memory allocated */ + else if (opt == 2 && flagvalue == NULL) + { + fprintf(stderr, "\nMEMORY_ERROR: %s() failed - returned NULL pointer\n\n", + funcname); + return 1; + } + + return 0; +} + +/*---- end of file ----*/ diff --git a/examples/arkode/C_serial/ark_brusselator_lsrk_externaldomeigest.out b/examples/arkode/C_serial/ark_brusselator_lsrk_externaldomeigest.out new file mode 100644 index 0000000000..ccc4759935 --- /dev/null +++ b/examples/arkode/C_serial/ark_brusselator_lsrk_externaldomeigest.out @@ -0,0 +1,39 @@ + +Brusselator ODE test problem: + initial conditions: u0 = 1.2, v0 = 3.1, w0 = 3 + problem parameters: a = 1, b = 3.5, ep = 5e-06 + reltol = 1.0e-06, abstol = 1.0e-10 + + t u v w + ------------------------------------------- + 0.000000 1.200000 3.100000 3.000000 + 1.000000 1.103852 3.013128 3.499981 + 2.000000 0.687945 3.521426 3.499988 + 3.000000 0.409445 4.277957 3.499993 + 4.000000 0.367889 4.942066 3.499994 + 5.000000 0.413862 5.510684 3.499993 + 6.000000 0.589251 5.855724 3.499990 + 7.000000 4.756506 0.735406 3.499917 + 8.000000 1.813465 1.575773 3.499968 + 9.000000 0.527922 2.807379 3.499991 + 10.000000 0.305608 3.657399 3.499995 + ------------------------------------------- + +Final Statistics: +Current time = 10.0124653670818 +Steps = 752 +Step attempts = 762 +Stability limited steps = 0 +Accuracy limited steps = 762 +Error test fails = 10 +NLS step fails = 0 +Inequality constraint fails = 0 +Initial step size = 1.13978962117636e-08 +Last step size = 0.020921379243018 +Current step size = 0.0212806907861296 +RHS fn evals = 47265 +Number of dom_eig updates = 1 +Max. num. of stages used = 195 +Max. num. of stages allowed = 200 +Max. spectral radius = 202001.210790713 +Min. spectral radius = 202001.210790713 \ No newline at end of file diff --git a/examples/arkode/C_serial/ark_heat1D_adapt.out b/examples/arkode/C_serial/ark_heat1D_adapt.out index 486aa1afa7..22318e1d33 100644 --- a/examples/arkode/C_serial/ark_heat1D_adapt.out +++ b/examples/arkode/C_serial/ark_heat1D_adapt.out @@ -1,4 +1,4 @@ - + 1D adaptive Heat PDE test problem: diffusion coefficient: k = 0.5 initial N = 21 @@ -72,7 +72,7 @@ 65 2.014255531625573e-02 2.014255531625573e-02 3.250909359584335e-02 61 5 105 ---------------------------------------------------------------------------------------- Final solver statistics: - Total number of time steps = 65 - Total nonlinear iterations = 433 - Total linear iterations = 8275 + Total number of time steps = 60 + Total nonlinear iterations = 370 + Total linear iterations = 6972 diff --git a/include/arkode/arkode_lsrkstep.h b/include/arkode/arkode_lsrkstep.h index 050f9d9905..09412bdceb 100644 --- a/include/arkode/arkode_lsrkstep.h +++ b/include/arkode/arkode_lsrkstep.h @@ -93,6 +93,9 @@ SUNDIALS_EXPORT int LSRKStepSetMaxNumStages(void* arkode_mem, SUNDIALS_EXPORT int LSRKStepSetDomEigSafetyFactor(void* arkode_mem, sunrealtype dom_eig_safety); +SUNDIALS_EXPORT int LSRKStepSetUseAnalyticStabRegion( + void* arkode_mem, sunbooleantype analytic_stab_region); + SUNDIALS_EXPORT int LSRKStepSetNumDomEigEstInitPreprocessIters(void* arkode_mem, int num_iters); diff --git a/include/sundials/sundials_domeigestimator.h b/include/sundials/sundials_domeigestimator.h index 0723678d21..fd27517e35 100644 --- a/include/sundials/sundials_domeigestimator.h +++ b/include/sundials/sundials_domeigestimator.h @@ -32,6 +32,9 @@ extern "C" { #endif +typedef int (*SUNRhsFn)(sunrealtype t, N_Vector y, N_Vector ydot, + void* user_data); + /* ----------------------------------------------------------------- * Generic definition of SUNDomEigEstimator (DEE) * ----------------------------------------------------------------- */ @@ -46,6 +49,9 @@ typedef _SUNDIALS_STRUCT_ SUNDomEigEstimator_* SUNDomEigEstimator; struct SUNDomEigEstimator_Ops_ { SUNErrCode (*setatimes)(SUNDomEigEstimator, void*, SUNATimesFn); + SUNErrCode (*setrhs)(SUNDomEigEstimator, void*, SUNRhsFn); + SUNErrCode (*setrhslinearizationpoint)(SUNDomEigEstimator, sunrealtype, + N_Vector); SUNErrCode (*setoptions)(SUNDomEigEstimator DEE, const char* Did, const char* file_name, int argc, char* argv[]); SUNErrCode (*setmaxiters)(SUNDomEigEstimator, long int); @@ -86,6 +92,14 @@ SUNDIALS_EXPORT SUNErrCode SUNDomEigEstimator_SetATimes(SUNDomEigEstimator DEE, void* A_data, SUNATimesFn ATimes); +SUNDIALS_EXPORT +SUNErrCode SUNDomEigEstimator_SetRhs(SUNDomEigEstimator DEE, void* rhs_data, + SUNRhsFn RHSfn); + +SUNDIALS_EXPORT +SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint(SUNDomEigEstimator DEE, + sunrealtype t, N_Vector v); + 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 2a451924dc..846aad3588 100644 --- a/include/sundomeigest/sundomeigest_arnoldi.h +++ b/include/sundomeigest/sundomeigest_arnoldi.h @@ -32,6 +32,10 @@ extern "C" { #endif +#ifndef MAX_DQITERS +#define MAX_DQITERS 3 +#endif + /* ----------------------------------------------------- * Arnoldi Iteration Implementation of SUNDomEigEstimator * ----------------------------------------------------- */ @@ -43,14 +47,21 @@ struct SUNDomEigEstimatorContent_Arnoldi_ /* Krylov subspace vectors */ N_Vector* V; - N_Vector q; + N_Vector q, 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 */ + 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 */ + sunrealtype tol_preprocess; /* Tolerance for preprocessing iterations */ + sunrealtype rhs_linT; /* Time value for linearization point */ long int num_ATimes; /* Number of ATimes calls */ + SUNRhsFn rhsfn; /* User provided RHS function */ + void* rhs_data; /* RHS function data */ + long int nfevals; /* Number of RHS evaluations */ + sunrealtype* LAPACK_A; /* The vector which holds rows of the Hessenberg matrix in the given order */ sunrealtype* LAPACK_wr; /* Real parts of eigenvalues */ sunrealtype* LAPACK_wi; /* Imaginary parts of eigenvalues */ @@ -75,10 +86,22 @@ SUNDIALS_EXPORT SUNErrCode SUNDomEigEstimator_SetATimes_Arnoldi(SUNDomEigEstimator DEE, void* A_data, SUNATimesFn ATimes); +SUNDIALS_EXPORT +SUNErrCode SUNDomEigEstimator_SetRhs_Arnoldi(SUNDomEigEstimator DEE, + void* rhs_data, SUNRhsFn RHSfn); + +SUNDIALS_EXPORT +SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi( + SUNDomEigEstimator DEE, sunrealtype t, N_Vector v); + SUNDIALS_EXPORT SUNErrCode SUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(SUNDomEigEstimator DEE, int num_iters); +SUNDIALS_EXPORT +SUNErrCode SUNDomEigEstimator_SetRelTol_Arnoldi(SUNDomEigEstimator DEE, + sunrealtype tol); + SUNDIALS_EXPORT SUNErrCode SUNDomEigEstimator_SetInitialGuess_Arnoldi(SUNDomEigEstimator DEE, N_Vector q); diff --git a/include/sundomeigest/sundomeigest_power.h b/include/sundomeigest/sundomeigest_power.h index 50755a19dc..421f05166a 100644 --- a/include/sundomeigest/sundomeigest_power.h +++ b/include/sundomeigest/sundomeigest_power.h @@ -32,6 +32,10 @@ extern "C" { #endif +#ifndef MAX_DQITERS +#define MAX_DQITERS 3 +#endif + /* ----------------------------------------------------- * Power Iteration Implementation of SUNDomEigEstimator * ----------------------------------------------------- */ @@ -41,16 +45,23 @@ struct SUNDomEigEstimatorContent_Power_ SUNATimesFn ATimes; /* User provided ATimes function */ void* ATdata; /* ATimes function data*/ - N_Vector V, q; /* workspace vectors */ + N_Vector V, q, q_prev, rhs_linY, Fy, work; /* workspace vectors */ - int num_warmups; /* Number of preprocessing iterations */ - long int max_iters; /* Maximum number of power iterations */ - long int num_iters; /* Number of iterations in last Estimate call */ + int num_warmups; /* Number of preprocessing iterations */ + long int max_iters; /* Maximum number of power iterations */ + long int num_iters; /* Number of iterations in last Estimate call */ + sunrealtype rhs_linT; /* Time value for linearization point */ long int num_ATimes; /* Number of ATimes calls */ sunrealtype rel_tol; /* Convergence criteria for the power iteration */ sunrealtype res; /* Residual from the last Estimate call */ + + SUNRhsFn rhsfn; /* User provided RHS function */ + void* rhs_data; /* RHS function data */ + long int nfevals; /* Number of RHS evaluations */ + + sunbooleantype complex; /* Flag for complex eigenvalue request */ }; typedef struct SUNDomEigEstimatorContent_Power_* SUNDomEigEstimatorContent_Power; @@ -68,6 +79,10 @@ SUNDIALS_EXPORT SUNErrCode SUNDomEigEstimator_SetATimes_Power(SUNDomEigEstimator DEE, void* A_data, SUNATimesFn ATimes); +SUNDIALS_EXPORT +SUNErrCode SUNDomEigEstimator_SetRhs_Power(SUNDomEigEstimator DEE, + void* rhs_data, SUNRhsFn RHSfn); + SUNDIALS_EXPORT SUNErrCode SUNDomEigEstimator_SetMaxIters_Power(SUNDomEigEstimator DEE, long int max_iters); @@ -84,6 +99,14 @@ SUNDIALS_EXPORT SUNErrCode SUNDomEigEstimator_SetInitialGuess_Power(SUNDomEigEstimator DEE, N_Vector q); +SUNDIALS_EXPORT +SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint_Power( + SUNDomEigEstimator DEE, sunrealtype t, N_Vector v); + +SUNDIALS_EXPORT +SUNErrCode SUNDomEigEstimator_SetIsReal_Power(SUNDomEigEstimator DEE, + sunbooleantype real); + SUNDIALS_EXPORT SUNErrCode SUNDomEigEstimator_Initialize_Power(SUNDomEigEstimator DEE); diff --git a/src/arkode/arkode_lsrkstep.c b/src/arkode/arkode_lsrkstep.c index 52c9d92d89..895cbd6761 100644 --- a/src/arkode/arkode_lsrkstep.c +++ b/src/arkode/arkode_lsrkstep.c @@ -562,8 +562,9 @@ int lsrkStep_TakeStepRKC(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) sunrealtype hmax, w0, w1, temp1, temp2, arg, bjm1, bjm2, mus, thjm1, thjm2, zjm1, zjm2, dzjm1, dzjm2, d2zjm1, d2zjm2, zj, dzj, d2zj, bj, ajm1, mu, nu, thj; - const sunrealtype onep54 = SUN_RCONST(1.54), c13 = SUN_RCONST(13.0), - p8 = SUN_RCONST(0.8), p4 = SUN_RCONST(0.4); + sunrealtype stability_norm; + + const sunrealtype p8 = SUN_RCONST(0.8), p4 = SUN_RCONST(0.4); ARKodeLSRKStepMem step_mem; /* initialize algebraic solver convergence flag to success, @@ -580,6 +581,9 @@ int lsrkStep_TakeStepRKC(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) N_Vector tmp1 = ark_mem->tempv1; N_Vector tmp2 = ark_mem->tempv2; + const sunrealtype coefz = + THREE / TWO / (ONE - TWO / SUN_RCONST(15.0) * step_mem->rkc_damping); + /* Initialize the current stage index */ step_mem->istage = 0; step_mem->req_stages = step_mem->stage_max_limit; @@ -591,10 +595,24 @@ int lsrkStep_TakeStepRKC(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) if (retval != ARK_SUCCESS) { return retval; } } - sunrealtype ss = - SUNRceil(SUNRsqrt(onep54 * SUNRabs(ark_mem->h) * step_mem->spectral_radius)); - ss = SUNMAX(ss, SUN_RCONST(2.0)); - + /* Compute the number of stages based on the current step size and dominant + eigenvalue using Eq. 2.7 in Verwer et al. (2004) + https://doi.org/10.1016/j.jcp.2004.05.002 + + Note beta(s) in Eq. 2.7 is positive (i.e., beta = -zR = h * lambdaR assuming + that h * lambdaR < 0) and we have incorporated the minus sign on zR below. We + rely on SUNRsqrt(x) returning 0 for x <= 0, to ensure we use the minimum number + of stages (ss = 2) for zR > 0. */ + sunrealtype zR = ark_mem->h * step_mem->lambdaR; + sunrealtype zI = ark_mem->h * step_mem->lambdaI; + sunrealtype ss = zR > ZERO ? SUN_RCONST(2.0) + : SUNRceil(SUNRsqrt(ONE - coefz * zR)); + ss = SUNMAX(ss, SUN_RCONST(2.0)); + + /* Check if number of stages exceeds maximum allowed. + If so, and if adaptive stepping is enabled, reduce step size + and return ARK_RETRY_STEP. If fixed step size, return + ARK_MAX_STAGE_LIMIT_FAIL error. */ if (ss >= step_mem->stage_max_limit) { SUNLogInfo(ARK_LOGGER, "compute-num-stages", @@ -605,8 +623,9 @@ int lsrkStep_TakeStepRKC(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) if (!ark_mem->fixedstep) { - hmax = ark_mem->hadapt_mem->safety * SUNSQR(step_mem->stage_max_limit) / - (onep54 * step_mem->spectral_radius); + hmax = ark_mem->hadapt_mem->safety * + (ONE - SUNSQR(step_mem->stage_max_limit)) / + (coefz * step_mem->lambdaR); ark_mem->eta = hmax / ark_mem->h; *nflagPtr = ARK_RETRY_STEP; ark_mem->hadapt_mem->nst_exp++; @@ -622,8 +641,82 @@ int lsrkStep_TakeStepRKC(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) } } - step_mem->req_stages = (int)ss; - step_mem->stage_max = SUNMAX(step_mem->req_stages, step_mem->stage_max); + int req_stages = (int)ss; + + if (zR < ZERO) + { + /* We first check whether the combination of ss, step size, and dominant + eigenvalue, is stable. If not, we then check whether it would be stable + when using ss = stage_max_limit -- if so, we increase ss until stability is + obtained. Otherwise, we reject the step, resulting either in method failure + when using fixed step sizes, or time step reduction when using adaptive + steps. */ + retval = lsrkStep_RKC_CheckStabilityNorm(step_mem, req_stages, ark_mem->h, + &stability_norm); + if (retval != ARK_SUCCESS) { return retval; } + + if (stability_norm > ONE - SUN_UNIT_ROUNDOFF) + { + sunrealtype initial_stability_norm = stability_norm; + sunbooleantype max_stage_is_stable = SUNFALSE; + + if (req_stages < step_mem->stage_max_limit) + { + retval = lsrkStep_RKC_CheckStabilityNorm(step_mem, + step_mem->stage_max_limit, + ark_mem->h, &stability_norm); + if (retval != ARK_SUCCESS) { return retval; } + + max_stage_is_stable = (stability_norm <= ONE - SUN_UNIT_ROUNDOFF); + stability_norm = initial_stability_norm; + } + + if (max_stage_is_stable) + { + while ((stability_norm > ONE - SUN_UNIT_ROUNDOFF) && + (req_stages < step_mem->stage_max_limit)) + { + req_stages += 1; + retval = lsrkStep_RKC_CheckStabilityNorm(step_mem, req_stages, + ark_mem->h, &stability_norm); + if (retval != ARK_SUCCESS) { return retval; } + } + } + + if (stability_norm > ONE - SUN_UNIT_ROUNDOFF) + { + if (!ark_mem->fixedstep) + { + /* For adaptive simulations, we adjust the step size by the ellipse approximation */ + const sunrealtype a = + (TWO / THREE) * (SUNSQR(ss) - ONE) * + (ONE - TWO / SUN_RCONST(15.0) * step_mem->rkc_damping) / TWO; + const sunrealtype b = a / (ss == 2 ? SUN_RCONST(0.6) + : SUN_RCONST(1.825) * ss); + + //TODO: Get opinions for negative eta if we remove the positivity restriction above for zR. + //Should we enforce a minimum eta value here? + ark_mem->eta = ark_mem->hadapt_mem->safety * (-TWO * a * b * b * zR) / + (SUNSQR(b * zR) + SUNSQR(a * zI)); + *nflagPtr = ARK_RETRY_STEP; + ark_mem->hadapt_mem->nst_exp++; + return ARK_RETRY_STEP; + } + else + { + arkProcessError(ark_mem, ARK_MAX_STAGE_LIMIT_FAIL, __LINE__, __func__, + __FILE__, + "Unable to achieve stable results: Either reduce the " + "step size or increase the stage_max_limit"); + return ARK_MAX_STAGE_LIMIT_FAIL; + } + } + } + } + + step_mem->req_stages = req_stages; + + step_mem->stage_max = SUNMAX(step_mem->req_stages, step_mem->stage_max); SUNLogInfo(ARK_LOGGER, "compute-num-stages", "spectral radius = " SUN_FORMAT_G @@ -671,7 +764,7 @@ int lsrkStep_TakeStepRKC(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) step_mem->step_nst = ark_mem->nst + 1; /* Initialize constants */ - w0 = (ONE + TWO / (c13 * SUNSQR((sunrealtype)(step_mem->req_stages)))); + w0 = (ONE + step_mem->rkc_damping / SUNSQR((sunrealtype)(step_mem->req_stages))); temp1 = SUNSQR(w0) - ONE; temp2 = SUNRsqrt(temp1); arg = step_mem->req_stages * SUNRlog(w0 + temp2); @@ -914,6 +1007,7 @@ int lsrkStep_TakeStepRKL(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) int retval; sunrealtype hmax, w1, bjm1, bjm2, mus, bj, ajm1, temj, cj, mu, nu; const sunrealtype p8 = SUN_RCONST(0.8), p4 = SUN_RCONST(0.4); + sunrealtype stability_norm; ARKodeLSRKStepMem step_mem; /* initialize algebraic solver convergence flag to success, @@ -941,14 +1035,23 @@ int lsrkStep_TakeStepRKL(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) if (retval != ARK_SUCCESS) { return retval; } } + /* Compute number of stages based on current step size and + dominant eigenvalue using Eq. 21 in Meyer et al. (2014) */ + sunrealtype zR = ark_mem->h * step_mem->lambdaR; + sunrealtype zI = ark_mem->h * step_mem->lambdaI; + sunrealtype zRabs = SUNRabs(zR); sunrealtype ss = - SUNRceil((SUNRsqrt(SUN_RCONST(9.0) + SUN_RCONST(8.0) * SUNRabs(ark_mem->h) * - step_mem->spectral_radius) - - ONE) / - TWO); + zR > ZERO + ? SUN_RCONST(2.0) + : SUNRceil((SUNRsqrt(SUN_RCONST(9.0) + SUN_RCONST(8.0) * zRabs) - ONE) / + TWO); ss = SUNMAX(ss, SUN_RCONST(2.0)); + /* Check if number of stages exceeds maximum allowed. + If so, and if adaptive stepping is enabled, reduce step size + and return ARK_RETRY_STEP. If fixed step size, return + ARK_MAX_STAGE_LIMIT_FAIL error. */ if (ss >= step_mem->stage_max_limit) { SUNLogInfo(ARK_LOGGER, "compute-num-stages", @@ -962,7 +1065,7 @@ int lsrkStep_TakeStepRKL(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) hmax = ark_mem->hadapt_mem->safety * (SUNSQR(step_mem->stage_max_limit) + step_mem->stage_max_limit - TWO) / - (TWO * step_mem->spectral_radius); + (-TWO * step_mem->lambdaR); ark_mem->eta = hmax / ark_mem->h; *nflagPtr = ARK_RETRY_STEP; ark_mem->hadapt_mem->nst_exp++; @@ -978,8 +1081,90 @@ int lsrkStep_TakeStepRKL(ARKodeMem ark_mem, sunrealtype* dsmPtr, int* nflagPtr) } } - step_mem->req_stages = (int)ss; - step_mem->stage_max = SUNMAX(step_mem->req_stages, step_mem->stage_max); + int req_stages = (int)ss; + + if (zR < ZERO) + { + /* To check stability, we evaluate the analytic stability function or an + inscribed ellipse approximation. If the stability norm is greater than + one, first check whether the method is stable at stage_max_limit. If so, + increase the number of stages until stability is obtained. Otherwise, + keep the existing fixed-step error and adaptive-step eta update logic. */ + retval = lsrkStep_RKL_CheckStabilityNorm(step_mem, req_stages, ark_mem->h, + &stability_norm); + if (retval != ARK_SUCCESS) { return retval; } + + if (stability_norm > ONE - SUN_UNIT_ROUNDOFF) + { + sunrealtype initial_stability_norm = stability_norm; + sunbooleantype max_stage_is_stable = SUNFALSE; + + if (req_stages < step_mem->stage_max_limit) + { + retval = lsrkStep_RKL_CheckStabilityNorm(step_mem, + step_mem->stage_max_limit, + ark_mem->h, &stability_norm); + if (retval != ARK_SUCCESS) { return retval; } + + max_stage_is_stable = (stability_norm <= ONE - SUN_UNIT_ROUNDOFF); + stability_norm = initial_stability_norm; + } + + if (max_stage_is_stable) + { + while ((stability_norm > ONE - SUN_UNIT_ROUNDOFF) && + (req_stages < step_mem->stage_max_limit)) + { + req_stages += 1; + retval = lsrkStep_RKL_CheckStabilityNorm(step_mem, req_stages, + ark_mem->h, &stability_norm); + if (retval != ARK_SUCCESS) { return retval; } + } + } + + if (stability_norm > ONE - SUN_UNIT_ROUNDOFF) + { + if (!ark_mem->fixedstep) + { + const sunrealtype aspect_ratio[7] = { + SUN_RCONST(0.3) * ss, /* s = 2 */ + SUN_RCONST(0.75) * ss, /* s = 3 */ + SUN_RCONST(0.665) * ss, /* s = 4 */ + SUN_RCONST(0.665) * ss, /* s = 5 */ + SUN_RCONST(0.635) * ss, /* s = 6 to 20 */ + SUN_RCONST(0.6) * ss, /* s >= 20 and odd */ + SUN_RCONST(0.53) * ss /* s >= 20 and even */ + }; + const sunrealtype a = + (((TWO * ss + ONE) * (TWO * ss + ONE)) - SUN_RCONST(9.0)) / + SUN_RCONST(16.0); + sunrealtype b; + + if (ss < 7) { b = a / aspect_ratio[(int)ss - 2]; } + else if (ss <= 20) { b = a / aspect_ratio[4]; } + else { b = a / aspect_ratio[6 - (int)ss % 2]; } + + ark_mem->eta = ark_mem->hadapt_mem->safety * (-TWO * a * b * b * zR) / + (SUNSQR(b * zR) + SUNSQR(a * zI)); + *nflagPtr = ARK_RETRY_STEP; + ark_mem->hadapt_mem->nst_exp++; + return ARK_RETRY_STEP; + } + else + { + arkProcessError(ark_mem, ARK_MAX_STAGE_LIMIT_FAIL, __LINE__, __func__, + __FILE__, + "Unable to achieve stable results: Either reduce the " + "step size or increase the stage_max_limit"); + return ARK_MAX_STAGE_LIMIT_FAIL; + } + } + } + } + + step_mem->req_stages = req_stages; + + step_mem->stage_max = SUNMAX(step_mem->req_stages, step_mem->stage_max); SUNLogInfo(ARK_LOGGER, "compute-num-stages", "spectral radius = " SUN_FORMAT_G @@ -2580,12 +2765,16 @@ void lsrkStep_PrintMem(ARKodeMem ark_mem, FILE* outfile) step_mem->spectral_radius_min); fprintf(outfile, "LSRKStep: dom_eig_safety = " SUN_FORMAT_G "\n", step_mem->dom_eig_safety); + fprintf(outfile, "LSRKStep: rkc_damping = " SUN_FORMAT_G "\n", + step_mem->rkc_damping); /* output sunbooleantype quantities */ fprintf(outfile, "LSRKStep: dom_eig_update = %d\n", step_mem->dom_eig_update); fprintf(outfile, "LSRKStep: dom_eig_is_current = %d\n", step_mem->dom_eig_is_current); + fprintf(outfile, "LSRKStep: use_ellipse = %d\n", + step_mem->use_ellipse); if (step_mem->DEE != NULL) { @@ -2755,7 +2944,7 @@ int lsrkStep_ComputeNewDomEig(ARKodeMem ark_mem, ARKodeLSRKStepMem step_mem) return ARK_DOMEIG_FAIL; } - if (step_mem->lambdaR * ark_mem->h > ZERO) + if (step_mem->lambdaR * ark_mem->h > SUNRsqrt(SUN_UNIT_ROUNDOFF)) { arkProcessError(NULL, ARK_DOMEIG_FAIL, __LINE__, __func__, __FILE__, "lambdaR*h must be nonpositive"); @@ -2793,6 +2982,250 @@ int lsrkStep_ComputeNewDomEig(ARKodeMem ark_mem, ARKodeLSRKStepMem step_mem) return retval; } +/*--------------------------------------------------------------- + lsrkStep_RKC_CheckStabilityNorm: + + This routine computes the stability norm for RKC methods. + If use_ellipse is SUNTRUE, we use a heuristic that approximates the stability region by an ellipse. + If use_ellipse is SUNFALSE, we compute the stability norm directly from the stability function using + the Chebyshev polynomial. + ---------------------------------------------------------------*/ +int lsrkStep_RKC_CheckStabilityNorm(ARKodeLSRKStepMem step_mem, int num_stages, + sunrealtype h, sunrealtype* stability_norm) +{ + sunrealtype ss = (sunrealtype)num_stages; + sunrealtype w0, w1, wr, wi, th, sh, ch, b_s, a_s, Ts, Ts_p, Ts_pp; + sunrealtype zR = SUNRabs(h) * step_mem->lambdaR; + sunrealtype zI = SUNRabs(h) * step_mem->lambdaI; + + if (step_mem->use_ellipse) + { + /* The stability region of the damped RKC method is approximated by an ellipse + centered at (-a,0), with horizontal semiaxis a and vertical semiaxis b, so + that its vertices are located at (0,0), (-2a,0), and (-a,+/-b). These + quantities depend on the damping parameter. Also, b is estimated + heuristically from the ellipse aspect ratio, taken as approximately 1.825s, + where s is the number of stages (for s=2, the ratio is approximated as 0.6). + This heuristic reflects the observed near-linear growth of the imaginary + extent with the number of stages. The numerical factors (1.825 and 0.6) + were obtained empirically from stability-region plots using the default + damping parameter and may change if the damping is modified. */ + const sunrealtype a = (TWO / THREE) * (SUNSQR(ss) - ONE) * + (ONE - TWO / SUN_RCONST(15.0) * step_mem->rkc_damping) / + TWO; + const sunrealtype b = a / + (ss == 2 ? SUN_RCONST(0.6) : SUN_RCONST(1.825) * ss); + + *stability_norm = SUNRsqrt(SUNSQR((zR + a) / a) + SUNSQR((zI) / b)); + } + else + { + w0 = ONE + step_mem->rkc_damping / (ss * ss); + th = SUNRacosh(w0); + sh = SUNRsinh(th); + ch = SUNRcosh(th); + + Ts = SUNRcosh(ss * th); + Ts_p = ss * SUNRsinh(ss * th) / sh; + Ts_pp = (ss * ss * SUNRcosh(ss * th) / (sh * sh)) - + ss * ch * SUNRsinh(ss * th) / (sh * sh * sh); + + b_s = Ts_pp / (Ts_p * Ts_p); + a_s = ONE - b_s * Ts; + w1 = Ts_p / Ts_pp; + + wr = w0 + w1 * zR; + wi = w1 * zI; + + sunrealtype TsR, TsI, Ps_ZR, Ps_ZI; + int retval = lsrkStep_cheb_T_complex(num_stages, wr, wi, &TsR, &TsI); + if (retval != ARK_SUCCESS) { return retval; } + + Ps_ZR = a_s + b_s * TsR; + Ps_ZI = b_s * TsI; + + *stability_norm = SUNRsqrt(SUNSQR(Ps_ZR) + SUNSQR(Ps_ZI)); + } + + return ARK_SUCCESS; +} + +/*--------------------------------------------------------------- + lsrkStep_RKL_CheckStabilityNorm: + + This routine computes the stability norm for RKL methods. + If use_ellipse is SUNTRUE, we use a heuristic that approximates the stability region by an ellipse. + If use_ellipse is SUNFALSE, we compute the stability norm directly from the stability function using + the Chebyshev polynomial. + ---------------------------------------------------------------*/ +int lsrkStep_RKL_CheckStabilityNorm(ARKodeLSRKStepMem step_mem, int num_stages, + sunrealtype h, sunrealtype* stability_norm) +{ + sunrealtype ss = (sunrealtype)num_stages; + sunrealtype w1, wr, wi, a_s, b_s; + sunrealtype zR = SUNRabs(h) * step_mem->lambdaR; + sunrealtype zI = SUNRabs(h) * step_mem->lambdaI; + + if (step_mem->use_ellipse) + { + /* The stability region of the RKL method is approximated by an ellipse + centered at (-a,0), with horizontal semiaxis a and vertical semiaxis b, + so that its vertices are located at (0,0), (-2a,0), and (-a,+/-b). + The half-height b is estimated heuristically from the ellipse aspect + ratio a/b based on the number of stages as follows: + s = 2 -> 0.3 s + s = 3 -> 0.75 s + s = 4 -> 0.665 s + s = 5 -> 0.665 s + s = 6 to 20 -> 0.635 s + s >= 20 and odd -> 0.6 s + s >= 20 and even -> 0.53 s */ + const sunrealtype aspect_ratio[7] = { + SUN_RCONST(0.3) * ss, /* s = 2 */ + SUN_RCONST(0.75) * ss, /* s = 3 */ + SUN_RCONST(0.665) * ss, /* s = 4 */ + SUN_RCONST(0.665) * ss, /* s = 5 */ + SUN_RCONST(0.635) * ss, /* s = 6 to 20 */ + SUN_RCONST(0.6) * ss, /* s >= 20 and odd */ + SUN_RCONST(0.53) * ss /* s >= 20 and even */ + }; + const sunrealtype a = + (((TWO * ss + ONE) * (TWO * ss + ONE)) - SUN_RCONST(9.0)) / + SUN_RCONST(16.0); + sunrealtype b; + + if (ss < 7) { b = a / (aspect_ratio[(int)ss - 2]); } + else + { + if (ss <= 20) { b = a / (aspect_ratio[4]); } + else { b = a / (aspect_ratio[6 - (int)ss % 2]); } + } + + *stability_norm = SUNRsqrt(SUNSQR((zR + a) / a) + SUNSQR(zI / b)); + } + else + { + b_s = (ss * ss + ss - TWO) / (TWO * ss * (ss + ONE)); + a_s = ONE - b_s; + w1 = FOUR / (ss * ss + ss - TWO); // Eq.(15) in Meyer et al. (2014) + wr = ONE + w1 * zR; + wi = w1 * zI; + + sunrealtype PsR, PsI, Ps_ZR, Ps_ZI; + int retval = lsrkStep_legendre_P_complex(num_stages, wr, wi, &PsR, &PsI); + if (retval != ARK_SUCCESS) { return retval; } + + Ps_ZR = a_s + b_s * PsR; + Ps_ZI = b_s * PsI; + + *stability_norm = SUNRsqrt(SUNSQR(Ps_ZR) + SUNSQR(Ps_ZI)); + } + + return ARK_SUCCESS; +} + +/*--------------------------------------------------------------- + lsrkStep_cheb_T_complex: + + This routine computes the Chebyshev polynomial of the first kind + T_s(z) for complex argument z = zR + i*zI using the + recurrence relation: + T_0(z) = 1 + T_1(z) = z + T_{k+1}(z) = 2*z*T_k(z) - T_{k-1}(z), k = 1,...,s-1 + ---------------------------------------------------------------*/ +int lsrkStep_cheb_T_complex(int s, sunrealtype zR, sunrealtype zI, + sunrealtype* TsR, sunrealtype* TsI) +{ + if (s < 0) + { + arkProcessError(NULL, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, + "s cannot be negative"); + return ARK_ILL_INPUT; + } + else if (s == 0) + { + *TsR = ONE; + *TsI = ZERO; + return ARK_SUCCESS; + } + else if (s == 1) + { + *TsR = zR; + *TsI = zI; + return ARK_SUCCESS; + } + else + { + sunrealtype Tkm1R = ONE, Tkm1I = ZERO; // T_0(z) + sunrealtype TkR = zR, TkI = zI; // T_1(z) + sunrealtype Tkp1R, Tkp1I; + for (int k = 1; k < s; k++) + { + Tkp1R = SUN_RCONST(2.0) * (zR * TkR - zI * TkI) - Tkm1R; + Tkp1I = SUN_RCONST(2.0) * (zR * TkI + zI * TkR) - Tkm1I; + Tkm1R = TkR; + Tkm1I = TkI; + TkR = Tkp1R; + TkI = Tkp1I; + } + *TsR = TkR; + *TsI = TkI; + } + return ARK_SUCCESS; +} + +/*--------------------------------------------------------------- + lsrkStep_legendre_P_complex: + + This routine computes the Legendre polynomial P_s(z) for complex + argument z = zR + i*zI using the recurrence relation: + P_0(z) = 1 + P_1(z) = z + P_{k+1}(z) = ((2*k+1)*z*P_k(z) - k*P_{k-1}(z))/(k+1), k = 1,...,s-1 + ---------------------------------------------------------------*/ + +int lsrkStep_legendre_P_complex(int s, sunrealtype zR, sunrealtype zI, + sunrealtype* PsR, sunrealtype* PsI) +{ + if (s < 0) + { + arkProcessError(NULL, ARK_ILL_INPUT, __LINE__, __func__, __FILE__, + "s cannot be negative"); + return ARK_ILL_INPUT; + } + else if (s == 0) + { + *PsR = ONE; + *PsI = ZERO; + return ARK_SUCCESS; + } + else if (s == 1) + { + *PsR = zR; + *PsI = zI; + return ARK_SUCCESS; + } + else + { + sunrealtype Pkm1R = ONE, Pkm1I = ZERO; // P_0(z) + sunrealtype PkR = zR, PkI = zI; // P_1(z) + sunrealtype Pkp1R, Pkp1I; + for (int k = 1; k < s; k++) + { + Pkp1R = ((TWO * k + ONE) * (zR * PkR - zI * PkI) - k * Pkm1R) / (k + ONE); + Pkp1I = ((TWO * k + ONE) * (zR * PkI + zI * PkR) - k * Pkm1I) / (k + ONE); + Pkm1R = PkR; + Pkm1I = PkI; + PkR = Pkp1R; + PkI = Pkp1I; + } + *PsR = PkR; + *PsI = PkI; + } + return ARK_SUCCESS; +} + /*--------------------------------------------------------------- lsrkStep_DQJtimes: diff --git a/src/arkode/arkode_lsrkstep_impl.h b/src/arkode/arkode_lsrkstep_impl.h index 54f1563641..47e002e721 100644 --- a/src/arkode/arkode_lsrkstep_impl.h +++ b/src/arkode/arkode_lsrkstep_impl.h @@ -32,6 +32,7 @@ extern "C" { #define STAGE_MAX_LIMIT_DEFAULT 200 #define DOM_EIG_SAFETY_DEFAULT SUN_RCONST(1.01) +#define RKC_DAMPING_DEFAULT SUN_RCONST(2.0) / SUN_RCONST(13.0) #define DOM_EIG_FREQ_DEFAULT 25 #define DOM_EIG_NUM_WARMUPS_DEFAULT 0 #define DOM_EIG_NUM_INIT_WARMUPS_DEFAULT -1 /* use DEE's default value */ @@ -119,6 +120,29 @@ extern "C" { #endif #endif +/* + * ----------------------------------------------------------------- + * Function : SUNRacosh + * ----------------------------------------------------------------- + * Usage : sunrealtype acosh_x; + * acosh_x = SUNRacosh(x); + * ----------------------------------------------------------------- + * SUNRacosh(x) returns acosh(x) (the hyperbolic arcosine of x). + * ----------------------------------------------------------------- + */ +#ifndef SUNRacosh +#if defined(SUNDIALS_DOUBLE_PRECISION) +#define SUNRacosh(x) (acosh((x))) +#elif defined(SUNDIALS_SINGLE_PRECISION) +#define SUNRacosh(x) (acoshf((x))) +#elif defined(SUNDIALS_EXTENDED_PRECISION) +#define SUNRacosh(x) (acoshl((x))) +#else +#error \ + "SUNDIALS precision not defined, report to github.com/LLNL/sundials/issues" +#endif +#endif + /*=============================================================== LSRK time step module data structure ===============================================================*/ @@ -161,6 +185,7 @@ typedef struct ARKodeLSRKStepMemRec sunrealtype spectral_radius_max; /* max spectral radius*/ sunrealtype spectral_radius_min; /* min spectral radius*/ sunrealtype dom_eig_safety; /* some safety factor for the user provided dom_eig*/ + sunrealtype rkc_damping; /* damping parameter for RKC methods*/ long int dom_eig_freq; /* indicates dom_eig update after dom_eig_freq successful steps*/ int num_init_warmups; /* number of warm-ups in the first DEE estimates */ int num_warmups; /* number of warm-ups in succeeding DEE estimates */ @@ -171,8 +196,9 @@ typedef struct ARKodeLSRKStepMemRec sunbooleantype dom_eig_update; /* flag indicating new dom_eig is needed */ sunbooleantype const_Jac; /* flag indicating Jacobian is constant */ sunbooleantype dom_eig_is_current; /* SUNTRUE if dom_eig has been evaluated at tn */ - sunbooleantype is_SSP; /* flag indicating SSP method*/ - sunbooleantype init_warmup; /* flag indicating initial warm-up*/ + sunbooleantype use_ellipse; /* flag indicating whether to use ellipse or exact stability region for stability checks */ + sunbooleantype is_SSP; /* flag indicating SSP method*/ + sunbooleantype init_warmup; /* flag indicating initial warm-up*/ /* Reusable fused vector operation arrays */ sunrealtype* cvals; @@ -220,6 +246,14 @@ int lsrkStep_AccessStepMem(ARKodeMem ark_mem, const char* fname, void lsrkStep_DomEigUpdateLogic(ARKodeMem ark_mem, ARKodeLSRKStepMem step_mem, sunrealtype dsm, N_Vector fnew); int lsrkStep_ComputeNewDomEig(ARKodeMem ark_mem, ARKodeLSRKStepMem step_mem); +int lsrkStep_RKC_CheckStabilityNorm(ARKodeLSRKStepMem step_mem, int num_stages, + sunrealtype h, sunrealtype* stability_norm); +int lsrkStep_RKL_CheckStabilityNorm(ARKodeLSRKStepMem step_mem, int num_stages, + sunrealtype h, sunrealtype* stability_norm); +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); /*=============================================================== diff --git a/src/arkode/arkode_lsrkstep_io.c b/src/arkode/arkode_lsrkstep_io.c index a92cd4930a..e425211343 100644 --- a/src/arkode/arkode_lsrkstep_io.c +++ b/src/arkode/arkode_lsrkstep_io.c @@ -316,6 +316,29 @@ int LSRKStepSetDomEigSafetyFactor(void* arkode_mem, sunrealtype dom_eig_safety) return ARK_SUCCESS; } +/*--------------------------------------------------------------- + LSRKStepSetUseAnalyticStabRegion sets whether to use the ellipse or the exact + stability region for stability checks. + ---------------------------------------------------------------*/ +int LSRKStepSetUseAnalyticStabRegion(void* arkode_mem, + sunbooleantype use_analytic_stab_region) +{ + ARKodeMem ark_mem; + ARKodeLSRKStepMem step_mem; + int retval; + + /* access ARKodeMem and ARKodeLSRKStepMem structures */ + retval = lsrkStep_AccessARKODEStepMem(arkode_mem, __func__, &ark_mem, + &step_mem); + if (retval != ARK_SUCCESS) { return retval; } + + step_mem->dom_eig_update = SUNTRUE; + step_mem->dom_eig_is_current = SUNFALSE; + step_mem->use_ellipse = !use_analytic_stab_region; + + return ARK_SUCCESS; +} + /*--------------------------------------------------------------- LSRKStepSetNumDomEigEstInitPreprocessIters sets the number of the preprocessing iterations before the very first estimate call. @@ -761,9 +784,11 @@ int lsrkStep_SetDefaults(ARKodeMem ark_mem) /* Spectral info */ step_mem->dom_eig_safety = DOM_EIG_SAFETY_DEFAULT; step_mem->dom_eig_freq = DOM_EIG_FREQ_DEFAULT; + step_mem->rkc_damping = RKC_DAMPING_DEFAULT; step_mem->const_Jac = SUNFALSE; step_mem->num_init_warmups = DOM_EIG_NUM_INIT_WARMUPS_DEFAULT; step_mem->num_warmups = DOM_EIG_NUM_WARMUPS_DEFAULT; + step_mem->use_ellipse = SUNTRUE; /* Load the default SUNAdaptController */ retval = arkReplaceAdaptController(ark_mem, NULL, SUNTRUE); @@ -895,6 +920,10 @@ int lsrkStep_WriteParameters(ARKodeMem ark_mem, FILE* fp) step_mem->spectral_radius); fprintf(fp, " Safety factor for the dom eig = " SUN_FORMAT_G "\n", step_mem->dom_eig_safety); + fprintf(fp, " Use elliptical stability region = %i\n", + step_mem->use_ellipse); + fprintf(fp, " Damping factor for RKC = " SUN_FORMAT_G "\n", + step_mem->rkc_damping); fprintf(fp, " Max num of successful steps before new dom eig update = %li\n", step_mem->dom_eig_freq); fprintf(fp, " Number of first preprocessing warmups = %i\n", diff --git a/src/arkode/fmod_int32/farkode_lsrkstep_mod.c b/src/arkode/fmod_int32/farkode_lsrkstep_mod.c index d310719af8..95b9a8f923 100644 --- a/src/arkode/fmod_int32/farkode_lsrkstep_mod.c +++ b/src/arkode/fmod_int32/farkode_lsrkstep_mod.c @@ -431,6 +431,20 @@ SWIGEXPORT int _wrap_FLSRKStepSetDomEigSafetyFactor(void *farg1, double const *f } +SWIGEXPORT int _wrap_FLSRKStepSetUseAnalyticStabRegion(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)LSRKStepSetUseAnalyticStabRegion(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FLSRKStepSetNumDomEigEstInitPreprocessIters(void *farg1, int const *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/arkode/fmod_int32/farkode_lsrkstep_mod.f90 b/src/arkode/fmod_int32/farkode_lsrkstep_mod.f90 index d943d4904c..ec7392a5ac 100644 --- a/src/arkode/fmod_int32/farkode_lsrkstep_mod.f90 +++ b/src/arkode/fmod_int32/farkode_lsrkstep_mod.f90 @@ -56,6 +56,7 @@ module farkode_lsrkstep_mod public :: FLSRKStepSetDomEigFrequency public :: FLSRKStepSetMaxNumStages public :: FLSRKStepSetDomEigSafetyFactor + public :: FLSRKStepSetUseAnalyticStabRegion public :: FLSRKStepSetNumDomEigEstInitPreprocessIters public :: FLSRKStepSetNumDomEigEstPreprocessIters public :: FLSRKStepSetNumSSPStages @@ -193,6 +194,15 @@ function swigc_FLSRKStepSetDomEigSafetyFactor(farg1, farg2) & integer(C_INT) :: fresult end function +function swigc_FLSRKStepSetUseAnalyticStabRegion(farg1, farg2) & +bind(C, name="_wrap_FLSRKStepSetUseAnalyticStabRegion") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + function swigc_FLSRKStepSetNumDomEigEstInitPreprocessIters(farg1, farg2) & bind(C, name="_wrap_FLSRKStepSetNumDomEigEstInitPreprocessIters") & result(fresult) @@ -513,6 +523,22 @@ function FLSRKStepSetDomEigSafetyFactor(arkode_mem, dom_eig_safety) & swig_result = fresult end function +function FLSRKStepSetUseAnalyticStabRegion(arkode_mem, analytic_stab_region) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: analytic_stab_region +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = analytic_stab_region +fresult = swigc_FLSRKStepSetUseAnalyticStabRegion(farg1, farg2) +swig_result = fresult +end function + function FLSRKStepSetNumDomEigEstInitPreprocessIters(arkode_mem, num_iters) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/arkode/fmod_int64/farkode_lsrkstep_mod.c b/src/arkode/fmod_int64/farkode_lsrkstep_mod.c index d310719af8..95b9a8f923 100644 --- a/src/arkode/fmod_int64/farkode_lsrkstep_mod.c +++ b/src/arkode/fmod_int64/farkode_lsrkstep_mod.c @@ -431,6 +431,20 @@ SWIGEXPORT int _wrap_FLSRKStepSetDomEigSafetyFactor(void *farg1, double const *f } +SWIGEXPORT int _wrap_FLSRKStepSetUseAnalyticStabRegion(void *farg1, int const *farg2) { + int fresult ; + void *arg1 = (void *) 0 ; + int arg2 ; + int result; + + arg1 = (void *)(farg1); + arg2 = (int)(*farg2); + result = (int)LSRKStepSetUseAnalyticStabRegion(arg1,arg2); + fresult = (int)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FLSRKStepSetNumDomEigEstInitPreprocessIters(void *farg1, int const *farg2) { int fresult ; void *arg1 = (void *) 0 ; diff --git a/src/arkode/fmod_int64/farkode_lsrkstep_mod.f90 b/src/arkode/fmod_int64/farkode_lsrkstep_mod.f90 index d943d4904c..ec7392a5ac 100644 --- a/src/arkode/fmod_int64/farkode_lsrkstep_mod.f90 +++ b/src/arkode/fmod_int64/farkode_lsrkstep_mod.f90 @@ -56,6 +56,7 @@ module farkode_lsrkstep_mod public :: FLSRKStepSetDomEigFrequency public :: FLSRKStepSetMaxNumStages public :: FLSRKStepSetDomEigSafetyFactor + public :: FLSRKStepSetUseAnalyticStabRegion public :: FLSRKStepSetNumDomEigEstInitPreprocessIters public :: FLSRKStepSetNumDomEigEstPreprocessIters public :: FLSRKStepSetNumSSPStages @@ -193,6 +194,15 @@ function swigc_FLSRKStepSetDomEigSafetyFactor(farg1, farg2) & integer(C_INT) :: fresult end function +function swigc_FLSRKStepSetUseAnalyticStabRegion(farg1, farg2) & +bind(C, name="_wrap_FLSRKStepSetUseAnalyticStabRegion") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +integer(C_INT), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + function swigc_FLSRKStepSetNumDomEigEstInitPreprocessIters(farg1, farg2) & bind(C, name="_wrap_FLSRKStepSetNumDomEigEstInitPreprocessIters") & result(fresult) @@ -513,6 +523,22 @@ function FLSRKStepSetDomEigSafetyFactor(arkode_mem, dom_eig_safety) & swig_result = fresult end function +function FLSRKStepSetUseAnalyticStabRegion(arkode_mem, analytic_stab_region) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(C_PTR) :: arkode_mem +integer(C_INT), intent(in) :: analytic_stab_region +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +integer(C_INT) :: farg2 + +farg1 = arkode_mem +farg2 = analytic_stab_region +fresult = swigc_FLSRKStepSetUseAnalyticStabRegion(farg1, farg2) +swig_result = fresult +end function + function FLSRKStepSetNumDomEigEstInitPreprocessIters(arkode_mem, num_iters) & result(swig_result) use, intrinsic :: ISO_C_BINDING diff --git a/src/sundials/fmod_int32/fsundials_core_mod.c b/src/sundials/fmod_int32/fsundials_core_mod.c index 008f283619..c690e87f87 100644 --- a/src/sundials/fmod_int32/fsundials_core_mod.c +++ b/src/sundials/fmod_int32/fsundials_core_mod.c @@ -3856,6 +3856,38 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetATimes(SUNDomEigEstimator farg1, voi } +SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhs(SUNDomEigEstimator farg1, void *farg2, SUNRhsFn farg3) { + int fresult ; + SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; + void *arg2 = (void *) 0 ; + SUNRhsFn arg3 = (SUNRhsFn) 0 ; + SUNErrCode result; + + arg1 = (SUNDomEigEstimator)(farg1); + arg2 = (void *)(farg2); + arg3 = (SUNRhsFn)(farg3); + result = (SUNErrCode)SUNDomEigEstimator_SetRhs(arg1,arg2,arg3); + fresult = (SUNErrCode)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint(SUNDomEigEstimator farg1, double const *farg2, N_Vector farg3) { + int fresult ; + SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; + sunrealtype arg2 ; + N_Vector arg3 = (N_Vector) 0 ; + SUNErrCode result; + + arg1 = (SUNDomEigEstimator)(farg1); + arg2 = (sunrealtype)(*farg2); + arg3 = (N_Vector)(farg3); + result = (SUNErrCode)SUNDomEigEstimator_SetRhsLinearizationPoint(arg1,arg2,arg3); + 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 e85111d008..78e071862e 100644 --- a/src/sundials/fmod_int32/fsundials_core_mod.f90 +++ b/src/sundials/fmod_int32/fsundials_core_mod.f90 @@ -690,6 +690,8 @@ module fsundials_core_mod ! struct struct SUNDomEigEstimator_Ops_ type, bind(C), public :: SUNDomEigEstimator_Ops type(C_FUNPTR), public :: setatimes + type(C_FUNPTR), public :: setrhs + type(C_FUNPTR), public :: setrhslinearizationpoint type(C_FUNPTR), public :: setoptions type(C_FUNPTR), public :: setmaxiters type(C_FUNPTR), public :: setnumpreprocessiters @@ -713,6 +715,8 @@ module fsundials_core_mod public :: FSUNDomEigEstimator_NewEmpty public :: FSUNDomEigEstimator_FreeEmpty public :: FSUNDomEigEstimator_SetATimes + public :: FSUNDomEigEstimator_SetRhs + public :: FSUNDomEigEstimator_SetRhsLinearizationPoint public :: FSUNDomEigEstimator_SetMaxIters public :: FSUNDomEigEstimator_SetNumPreprocessIters public :: FSUNDomEigEstimator_SetRelTol @@ -2975,6 +2979,26 @@ function swigc_FSUNDomEigEstimator_SetATimes(farg1, farg2, farg3) & integer(C_INT) :: fresult end function +function swigc_FSUNDomEigEstimator_SetRhs(farg1, farg2, farg3) & +bind(C, name="_wrap_FSUNDomEigEstimator_SetRhs") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +type(C_FUNPTR), value :: farg3 +integer(C_INT) :: fresult +end function + +function swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint(farg1, farg2, farg3) & +bind(C, name="_wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +type(C_PTR), value :: farg3 +integer(C_INT) :: fresult +end function + function swigc_FSUNDomEigEstimator_SetMaxIters(farg1, farg2) & bind(C, name="_wrap_FSUNDomEigEstimator_SetMaxIters") & result(fresult) @@ -7229,6 +7253,44 @@ function FSUNDomEigEstimator_SetATimes(dee, a_data, atimes) & swig_result = fresult end function +function FSUNDomEigEstimator_SetRhs(dee, rhs_data, rhsfn) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(SUNDomEigEstimator), target, intent(inout) :: dee +type(C_PTR) :: rhs_data +type(C_FUNPTR), intent(in), value :: rhsfn +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 +type(C_FUNPTR) :: farg3 + +farg1 = c_loc(dee) +farg2 = rhs_data +farg3 = rhsfn +fresult = swigc_FSUNDomEigEstimator_SetRhs(farg1, farg2, farg3) +swig_result = fresult +end function + +function FSUNDomEigEstimator_SetRhsLinearizationPoint(dee, t, v) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(SUNDomEigEstimator), target, intent(inout) :: dee +real(C_DOUBLE), intent(in) :: t +type(N_Vector), target, intent(inout) :: v +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 +type(C_PTR) :: farg3 + +farg1 = c_loc(dee) +farg2 = t +farg3 = c_loc(v) +fresult = swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint(farg1, farg2, farg3) +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 0e57bd84a1..da329cd45e 100644 --- a/src/sundials/fmod_int64/fsundials_core_mod.c +++ b/src/sundials/fmod_int64/fsundials_core_mod.c @@ -3856,6 +3856,38 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetATimes(SUNDomEigEstimator farg1, voi } +SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhs(SUNDomEigEstimator farg1, void *farg2, SUNRhsFn farg3) { + int fresult ; + SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; + void *arg2 = (void *) 0 ; + SUNRhsFn arg3 = (SUNRhsFn) 0 ; + SUNErrCode result; + + arg1 = (SUNDomEigEstimator)(farg1); + arg2 = (void *)(farg2); + arg3 = (SUNRhsFn)(farg3); + result = (SUNErrCode)SUNDomEigEstimator_SetRhs(arg1,arg2,arg3); + fresult = (SUNErrCode)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint(SUNDomEigEstimator farg1, double const *farg2, N_Vector farg3) { + int fresult ; + SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; + sunrealtype arg2 ; + N_Vector arg3 = (N_Vector) 0 ; + SUNErrCode result; + + arg1 = (SUNDomEigEstimator)(farg1); + arg2 = (sunrealtype)(*farg2); + arg3 = (N_Vector)(farg3); + result = (SUNErrCode)SUNDomEigEstimator_SetRhsLinearizationPoint(arg1,arg2,arg3); + 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 4c52de03de..9434bee137 100644 --- a/src/sundials/fmod_int64/fsundials_core_mod.f90 +++ b/src/sundials/fmod_int64/fsundials_core_mod.f90 @@ -690,6 +690,8 @@ module fsundials_core_mod ! struct struct SUNDomEigEstimator_Ops_ type, bind(C), public :: SUNDomEigEstimator_Ops type(C_FUNPTR), public :: setatimes + type(C_FUNPTR), public :: setrhs + type(C_FUNPTR), public :: setrhslinearizationpoint type(C_FUNPTR), public :: setoptions type(C_FUNPTR), public :: setmaxiters type(C_FUNPTR), public :: setnumpreprocessiters @@ -713,6 +715,8 @@ module fsundials_core_mod public :: FSUNDomEigEstimator_NewEmpty public :: FSUNDomEigEstimator_FreeEmpty public :: FSUNDomEigEstimator_SetATimes + public :: FSUNDomEigEstimator_SetRhs + public :: FSUNDomEigEstimator_SetRhsLinearizationPoint public :: FSUNDomEigEstimator_SetMaxIters public :: FSUNDomEigEstimator_SetNumPreprocessIters public :: FSUNDomEigEstimator_SetRelTol @@ -2975,6 +2979,26 @@ function swigc_FSUNDomEigEstimator_SetATimes(farg1, farg2, farg3) & integer(C_INT) :: fresult end function +function swigc_FSUNDomEigEstimator_SetRhs(farg1, farg2, farg3) & +bind(C, name="_wrap_FSUNDomEigEstimator_SetRhs") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +type(C_FUNPTR), value :: farg3 +integer(C_INT) :: fresult +end function + +function swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint(farg1, farg2, farg3) & +bind(C, name="_wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +type(C_PTR), value :: farg3 +integer(C_INT) :: fresult +end function + function swigc_FSUNDomEigEstimator_SetMaxIters(farg1, farg2) & bind(C, name="_wrap_FSUNDomEigEstimator_SetMaxIters") & result(fresult) @@ -7229,6 +7253,44 @@ function FSUNDomEigEstimator_SetATimes(dee, a_data, atimes) & swig_result = fresult end function +function FSUNDomEigEstimator_SetRhs(dee, rhs_data, rhsfn) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(SUNDomEigEstimator), target, intent(inout) :: dee +type(C_PTR) :: rhs_data +type(C_FUNPTR), intent(in), value :: rhsfn +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 +type(C_FUNPTR) :: farg3 + +farg1 = c_loc(dee) +farg2 = rhs_data +farg3 = rhsfn +fresult = swigc_FSUNDomEigEstimator_SetRhs(farg1, farg2, farg3) +swig_result = fresult +end function + +function FSUNDomEigEstimator_SetRhsLinearizationPoint(dee, t, v) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(SUNDomEigEstimator), target, intent(inout) :: dee +real(C_DOUBLE), intent(in) :: t +type(N_Vector), target, intent(inout) :: v +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 +type(C_PTR) :: farg3 + +farg1 = c_loc(dee) +farg2 = t +farg3 = c_loc(v) +fresult = swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint(farg1, farg2, farg3) +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 54efcc9002..aba7b5b66a 100644 --- a/src/sundials/sundials_domeigestimator.c +++ b/src/sundials/sundials_domeigestimator.c @@ -193,6 +193,31 @@ SUNErrCode SUNDomEigEstimator_SetATimes(SUNDomEigEstimator DEE, void* A_data, return (ier); } +SUNErrCode SUNDomEigEstimator_SetRhs(SUNDomEigEstimator DEE, void* rhs_data, + SUNRhsFn RHSfn) +{ + SUNErrCode ier; + SUNDIALS_MARK_FUNCTION_BEGIN(getSUNProfiler(DEE)); + if (DEE->ops->setrhs) { ier = DEE->ops->setrhs(DEE, rhs_data, RHSfn); } + else { ier = SUN_SUCCESS; } + SUNDIALS_MARK_FUNCTION_END(getSUNProfiler(DEE)); + return (ier); +} + +SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint(SUNDomEigEstimator DEE, + sunrealtype t, N_Vector v) +{ + SUNErrCode ier; + SUNDIALS_MARK_FUNCTION_BEGIN(getSUNProfiler(DEE)); + if (DEE->ops->setrhslinearizationpoint) + { + ier = DEE->ops->setrhslinearizationpoint(DEE, t, v); + } + 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 3c0799fc5a..295c851da4 100644 --- a/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.c +++ b/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.c @@ -396,6 +396,78 @@ SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Arnoldi__q_get(SwigClassWrap } +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 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhs_linY", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (N_Vector)(farg2); + if (arg1) (arg1)->rhs_linY = arg2; +} + + +SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_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_::rhs_linY", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (N_Vector) ((arg1)->rhs_linY); + fresult = result; + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_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_::Fy", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (N_Vector)(farg2); + if (arg1) (arg1)->Fy = arg2; +} + + +SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_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_::Fy", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (N_Vector) ((arg1)->Fy); + fresult = result; + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__work_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_::work", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (N_Vector)(farg2); + if (arg1) (arg1)->work = arg2; +} + + +SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Arnoldi__work_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_::work", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (N_Vector) ((arg1)->work); + fresult = result; + return fresult; +} + + SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__kry_dim_set(SwigClassWrapper const *farg1, int const *farg2) { struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; int arg2 ; @@ -468,6 +540,78 @@ SWIGEXPORT long _wrap_SUNDomEigEstimatorContent_Arnoldi__num_iters_get(SwigClass } +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_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_::warmup_to_tol", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (int)(*farg2); + if (arg1) (arg1)->warmup_to_tol = arg2; +} + + +SWIGEXPORT int _wrap_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_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_::warmup_to_tol", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (int) ((arg1)->warmup_to_tol); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set(SwigClassWrapper const *farg1, double const *farg2) { + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + sunrealtype arg2 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::tol_preprocess", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (sunrealtype)(*farg2); + if (arg1) (arg1)->tol_preprocess = arg2; +} + + +SWIGEXPORT double _wrap_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get(SwigClassWrapper const *farg1) { + double fresult ; + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + sunrealtype result; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::tol_preprocess", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (sunrealtype) ((arg1)->tol_preprocess); + fresult = (sunrealtype)(result); + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set(SwigClassWrapper const *farg1, double const *farg2) { + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + sunrealtype arg2 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhs_linT", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (sunrealtype)(*farg2); + if (arg1) (arg1)->rhs_linT = arg2; +} + + +SWIGEXPORT double _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get(SwigClassWrapper const *farg1) { + double fresult ; + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + sunrealtype result; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhs_linT", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (sunrealtype) ((arg1)->rhs_linT); + fresult = (sunrealtype)(result); + return fresult; +} + + SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_set(SwigClassWrapper const *farg1, long const *farg2) { struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; long arg2 ; @@ -492,6 +636,78 @@ SWIGEXPORT long _wrap_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_get(SwigClas } +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set(SwigClassWrapper const *farg1, SUNRhsFn farg2) { + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + SUNRhsFn arg2 = (SUNRhsFn) 0 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhsfn", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (SUNRhsFn)(farg2); + if (arg1) (arg1)->rhsfn = arg2; +} + + +SWIGEXPORT SUNRhsFn _wrap_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get(SwigClassWrapper const *farg1) { + SUNRhsFn fresult ; + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + SUNRhsFn result; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhsfn", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (SUNRhsFn) ((arg1)->rhsfn); + fresult = result; + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set(SwigClassWrapper const *farg1, void *farg2) { + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + void *arg2 = (void *) 0 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhs_data", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (void *)(farg2); + if (arg1) (arg1)->rhs_data = arg2; +} + + +SWIGEXPORT void * _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get(SwigClassWrapper const *farg1) { + void * fresult ; + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + void *result = 0 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhs_data", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (void *) ((arg1)->rhs_data); + fresult = result; + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__nfevals_set(SwigClassWrapper const *farg1, long const *farg2) { + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + long arg2 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::nfevals", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (long)(*farg2); + if (arg1) (arg1)->nfevals = arg2; +} + + +SWIGEXPORT long _wrap_SUNDomEigEstimatorContent_Arnoldi__nfevals_get(SwigClassWrapper const *farg1) { + long fresult ; + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + long result; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::nfevals", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (long) ((arg1)->nfevals); + fresult = (long)(result); + return fresult; +} + + SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_set(SwigClassWrapper const *farg1, double *farg2) { struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; sunrealtype *arg2 = (sunrealtype *) 0 ; @@ -723,6 +939,38 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetATimes_Arnoldi(SUNDomEigEstimator fa } +SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhs_Arnoldi(SUNDomEigEstimator farg1, void *farg2, SUNRhsFn farg3) { + int fresult ; + SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; + void *arg2 = (void *) 0 ; + SUNRhsFn arg3 = (SUNRhsFn) 0 ; + SUNErrCode result; + + arg1 = (SUNDomEigEstimator)(farg1); + arg2 = (void *)(farg2); + arg3 = (SUNRhsFn)(farg3); + result = (SUNErrCode)SUNDomEigEstimator_SetRhs_Arnoldi(arg1,arg2,arg3); + fresult = (SUNErrCode)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(SUNDomEigEstimator farg1, double const *farg2, N_Vector farg3) { + int fresult ; + SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; + sunrealtype arg2 ; + N_Vector arg3 = (N_Vector) 0 ; + SUNErrCode result; + + arg1 = (SUNDomEigEstimator)(farg1); + arg2 = (sunrealtype)(*farg2); + arg3 = (N_Vector)(farg3); + result = (SUNErrCode)SUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(arg1,arg2,arg3); + fresult = (SUNErrCode)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(SUNDomEigEstimator farg1, int const *farg2) { int fresult ; SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; @@ -737,6 +985,20 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(SUNDomEig } +SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRelTol_Arnoldi(SUNDomEigEstimator farg1, double const *farg2) { + int fresult ; + SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; + sunrealtype arg2 ; + SUNErrCode result; + + arg1 = (SUNDomEigEstimator)(farg1); + arg2 = (sunrealtype)(*farg2); + result = (SUNErrCode)SUNDomEigEstimator_SetRelTol_Arnoldi(arg1,arg2); + fresult = (SUNErrCode)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetInitialGuess_Arnoldi(SUNDomEigEstimator farg1, N_Vector 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 787459fa26..18adbe9fd2 100644 --- a/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.f90 +++ b/src/sundomeigest/arnoldi/fmod_int32/fsundomeigest_arnoldi_mod.f90 @@ -28,6 +28,7 @@ module fsundomeigest_arnoldi_mod private ! DECLARATION CONSTRUCTS + integer(C_INT), parameter, public :: MAX_DQITERS = 3_C_INT integer, parameter :: swig_cmem_own_bit = 0 integer, parameter :: swig_cmem_rvalue_bit = 1 @@ -48,14 +49,32 @@ module fsundomeigest_arnoldi_mod 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 + procedure :: get_Fy => swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_get + procedure :: set_work => swigf_SUNDomEigEstimatorContent_Arnoldi__work_set + procedure :: get_work => swigf_SUNDomEigEstimatorContent_Arnoldi__work_get procedure :: set_kry_dim => swigf_SUNDomEigEstimatorContent_Arnoldi__kry_dim_set procedure :: get_kry_dim => swigf_SUNDomEigEstimatorContent_Arnoldi__kry_dim_get procedure :: set_num_warmups => swigf_SUNDomEigEstimatorContent_Arnoldi__num_warmups_set procedure :: get_num_warmups => swigf_SUNDomEigEstimatorContent_Arnoldi__num_warmups_get procedure :: set_num_iters => swigf_SUNDomEigEstimatorContent_Arnoldi__num_iters_set 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_tol_preprocess => swigf_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set + procedure :: get_tol_preprocess => swigf_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get + procedure :: set_rhs_linT => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set + procedure :: get_rhs_linT => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get procedure :: set_num_ATimes => swigf_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_set procedure :: get_num_ATimes => swigf_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_get + procedure :: set_rhsfn => swigf_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set + procedure :: get_rhsfn => swigf_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get + procedure :: set_rhs_data => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set + procedure :: get_rhs_data => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get + procedure :: set_nfevals => swigf_SUNDomEigEstimatorContent_Arnoldi__nfevals_set + procedure :: get_nfevals => swigf_SUNDomEigEstimatorContent_Arnoldi__nfevals_get procedure :: set_LAPACK_A => swigf_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_set procedure :: get_LAPACK_A => swigf_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_get procedure :: set_LAPACK_wr => swigf_SUNDomEigEstimatorContent_Arnoldi__LAPACK_wr_set @@ -79,7 +98,10 @@ module fsundomeigest_arnoldi_mod end interface public :: FSUNDomEigEstimator_Arnoldi public :: FSUNDomEigEstimator_SetATimes_Arnoldi + public :: FSUNDomEigEstimator_SetRhs_Arnoldi + public :: FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi public :: FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi + public :: FSUNDomEigEstimator_SetRelTol_Arnoldi public :: FSUNDomEigEstimator_SetInitialGuess_Arnoldi public :: FSUNDomEigEstimator_Initialize_Arnoldi public :: FSUNDomEigEstimator_Estimate_Arnoldi @@ -158,6 +180,57 @@ function swigc_SUNDomEigEstimatorContent_Arnoldi__q_get(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 +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR), value :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR), value :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__work_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__work_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR), value :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__work_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__work_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: fresult +end function + subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__kry_dim_set(farg1, farg2) & bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__kry_dim_set") use, intrinsic :: ISO_C_BINDING @@ -209,6 +282,57 @@ function swigc_SUNDomEigEstimatorContent_Arnoldi__num_iters_get(farg1) & integer(C_LONG) :: fresult end function +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +integer(C_INT), intent(in) :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +integer(C_INT) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE) :: fresult +end function + subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_set(farg1, farg2) & bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_set") use, intrinsic :: ISO_C_BINDING @@ -226,6 +350,57 @@ function swigc_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_get(farg1) & integer(C_LONG) :: fresult end function +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_FUNPTR), value :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_FUNPTR) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR), value :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__nfevals_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__nfevals_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +integer(C_LONG), intent(in) :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__nfevals_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__nfevals_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +integer(C_LONG) :: fresult +end function + subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_set(farg1, farg2) & bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_set") use, intrinsic :: ISO_C_BINDING @@ -388,6 +563,26 @@ function swigc_FSUNDomEigEstimator_SetATimes_Arnoldi(farg1, farg2, farg3) & integer(C_INT) :: fresult end function +function swigc_FSUNDomEigEstimator_SetRhs_Arnoldi(farg1, farg2, farg3) & +bind(C, name="_wrap_FSUNDomEigEstimator_SetRhs_Arnoldi") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +type(C_FUNPTR), value :: farg3 +integer(C_INT) :: fresult +end function + +function swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(farg1, farg2, farg3) & +bind(C, name="_wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +type(C_PTR), value :: farg3 +integer(C_INT) :: fresult +end function + function swigc_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(farg1, farg2) & bind(C, name="_wrap_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi") & result(fresult) @@ -397,6 +592,15 @@ function swigc_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(farg1, farg2) & integer(C_INT) :: fresult end function +function swigc_FSUNDomEigEstimator_SetRelTol_Arnoldi(farg1, farg2) & +bind(C, name="_wrap_FSUNDomEigEstimator_SetRelTol_Arnoldi") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + function swigc_FSUNDomEigEstimator_SetInitialGuess_Arnoldi(farg1, farg2) & bind(C, name="_wrap_FSUNDomEigEstimator_SetInitialGuess_Arnoldi") & result(fresult) @@ -564,6 +768,81 @@ function swigf_SUNDomEigEstimatorContent_Arnoldi__q_get(self) & 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 +type(N_Vector), target, intent(inout) :: rhs_liny +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: farg2 + +farg1 = self%swigdata +farg2 = c_loc(rhs_liny) +call swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_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__rhs_linY_get(farg1) +call c_f_pointer(fresult, swig_result) +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_set(self, fy) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(N_Vector), target, intent(inout) :: fy +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: farg2 + +farg1 = self%swigdata +farg2 = c_loc(fy) +call swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_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__Fy_get(farg1) +call c_f_pointer(fresult, swig_result) +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__work_set(self, work) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(N_Vector), target, intent(inout) :: work +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: farg2 + +farg1 = self%swigdata +farg2 = c_loc(work) +call swigc_SUNDomEigEstimatorContent_Arnoldi__work_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__work_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__work_get(farg1) +call c_f_pointer(fresult, swig_result) +end function + subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__kry_dim_set(self, kry_dim) use, intrinsic :: ISO_C_BINDING class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self @@ -639,6 +918,81 @@ function swigf_SUNDomEigEstimatorContent_Arnoldi__num_iters_get(self) & swig_result = fresult end function +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_set(self, warmup_to_tol) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +integer(C_INT), intent(in) :: warmup_to_tol +type(SwigClassWrapper) :: farg1 +integer(C_INT) :: farg2 + +farg1 = self%swigdata +farg2 = warmup_to_tol +call swigc_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_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__warmup_to_tol_get(farg1) +swig_result = fresult +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set(self, tol_preprocess) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +real(C_DOUBLE), intent(in) :: tol_preprocess +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = self%swigdata +farg2 = tol_preprocess +call swigc_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +real(C_DOUBLE) :: swig_result +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +real(C_DOUBLE) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get(farg1) +swig_result = fresult +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set(self, rhs_lint) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +real(C_DOUBLE), intent(in) :: rhs_lint +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = self%swigdata +farg2 = rhs_lint +call swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +real(C_DOUBLE) :: swig_result +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +real(C_DOUBLE) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get(farg1) +swig_result = fresult +end function + subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_set(self, num_atimes) use, intrinsic :: ISO_C_BINDING class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self @@ -664,6 +1018,81 @@ function swigf_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_get(self) & swig_result = fresult end function +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set(self, rhsfn) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(C_FUNPTR), intent(in), value :: rhsfn +type(SwigClassWrapper) :: farg1 +type(C_FUNPTR) :: farg2 + +farg1 = self%swigdata +farg2 = rhsfn +call swigc_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +type(C_FUNPTR) :: swig_result +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(C_FUNPTR) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get(farg1) +swig_result = fresult +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set(self, rhs_data) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(C_PTR) :: rhs_data +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: farg2 + +farg1 = self%swigdata +farg2 = rhs_data +call swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +type(C_PTR) :: swig_result +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(C_PTR) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get(farg1) +swig_result = fresult +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__nfevals_set(self, nfevals) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +integer(C_LONG), intent(in) :: nfevals +type(SwigClassWrapper) :: farg1 +integer(C_LONG) :: farg2 + +farg1 = self%swigdata +farg2 = nfevals +call swigc_SUNDomEigEstimatorContent_Arnoldi__nfevals_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__nfevals_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_LONG) :: swig_result +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +integer(C_LONG) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__nfevals_get(farg1) +swig_result = fresult +end function + subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_set(self, lapack_a) use, intrinsic :: ISO_C_BINDING class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self @@ -914,6 +1343,44 @@ function FSUNDomEigEstimator_SetATimes_Arnoldi(dee, a_data, atimes) & swig_result = fresult end function +function FSUNDomEigEstimator_SetRhs_Arnoldi(dee, rhs_data, rhsfn) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(SUNDomEigEstimator), target, intent(inout) :: dee +type(C_PTR) :: rhs_data +type(C_FUNPTR), intent(in), value :: rhsfn +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 +type(C_FUNPTR) :: farg3 + +farg1 = c_loc(dee) +farg2 = rhs_data +farg3 = rhsfn +fresult = swigc_FSUNDomEigEstimator_SetRhs_Arnoldi(farg1, farg2, farg3) +swig_result = fresult +end function + +function FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(dee, t, v) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(SUNDomEigEstimator), target, intent(inout) :: dee +real(C_DOUBLE), intent(in) :: t +type(N_Vector), target, intent(inout) :: v +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 +type(C_PTR) :: farg3 + +farg1 = c_loc(dee) +farg2 = t +farg3 = c_loc(v) +fresult = swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(farg1, farg2, farg3) +swig_result = fresult +end function + function FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(dee, num_iters) & result(swig_result) use, intrinsic :: ISO_C_BINDING @@ -930,6 +1397,22 @@ function FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(dee, num_iters) & swig_result = fresult end function +function FSUNDomEigEstimator_SetRelTol_Arnoldi(dee, tol) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(SUNDomEigEstimator), target, intent(inout) :: dee +real(C_DOUBLE), intent(in) :: tol +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = c_loc(dee) +farg2 = tol +fresult = swigc_FSUNDomEigEstimator_SetRelTol_Arnoldi(farg1, farg2) +swig_result = fresult +end function + function FSUNDomEigEstimator_SetInitialGuess_Arnoldi(dee, q) & 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 6a49443ec8..98bddb6137 100644 --- a/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.c +++ b/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.c @@ -396,6 +396,78 @@ SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Arnoldi__q_get(SwigClassWrap } +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 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhs_linY", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (N_Vector)(farg2); + if (arg1) (arg1)->rhs_linY = arg2; +} + + +SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_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_::rhs_linY", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (N_Vector) ((arg1)->rhs_linY); + fresult = result; + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_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_::Fy", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (N_Vector)(farg2); + if (arg1) (arg1)->Fy = arg2; +} + + +SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_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_::Fy", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (N_Vector) ((arg1)->Fy); + fresult = result; + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__work_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_::work", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (N_Vector)(farg2); + if (arg1) (arg1)->work = arg2; +} + + +SWIGEXPORT N_Vector _wrap_SUNDomEigEstimatorContent_Arnoldi__work_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_::work", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (N_Vector) ((arg1)->work); + fresult = result; + return fresult; +} + + SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__kry_dim_set(SwigClassWrapper const *farg1, int const *farg2) { struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; int arg2 ; @@ -468,6 +540,78 @@ SWIGEXPORT long _wrap_SUNDomEigEstimatorContent_Arnoldi__num_iters_get(SwigClass } +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_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_::warmup_to_tol", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (int)(*farg2); + if (arg1) (arg1)->warmup_to_tol = arg2; +} + + +SWIGEXPORT int _wrap_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_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_::warmup_to_tol", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (int) ((arg1)->warmup_to_tol); + fresult = (int)(result); + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set(SwigClassWrapper const *farg1, double const *farg2) { + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + sunrealtype arg2 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::tol_preprocess", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (sunrealtype)(*farg2); + if (arg1) (arg1)->tol_preprocess = arg2; +} + + +SWIGEXPORT double _wrap_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get(SwigClassWrapper const *farg1) { + double fresult ; + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + sunrealtype result; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::tol_preprocess", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (sunrealtype) ((arg1)->tol_preprocess); + fresult = (sunrealtype)(result); + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set(SwigClassWrapper const *farg1, double const *farg2) { + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + sunrealtype arg2 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhs_linT", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (sunrealtype)(*farg2); + if (arg1) (arg1)->rhs_linT = arg2; +} + + +SWIGEXPORT double _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get(SwigClassWrapper const *farg1) { + double fresult ; + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + sunrealtype result; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhs_linT", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (sunrealtype) ((arg1)->rhs_linT); + fresult = (sunrealtype)(result); + return fresult; +} + + SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_set(SwigClassWrapper const *farg1, long const *farg2) { struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; long arg2 ; @@ -492,6 +636,78 @@ SWIGEXPORT long _wrap_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_get(SwigClas } +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set(SwigClassWrapper const *farg1, SUNRhsFn farg2) { + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + SUNRhsFn arg2 = (SUNRhsFn) 0 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhsfn", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (SUNRhsFn)(farg2); + if (arg1) (arg1)->rhsfn = arg2; +} + + +SWIGEXPORT SUNRhsFn _wrap_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get(SwigClassWrapper const *farg1) { + SUNRhsFn fresult ; + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + SUNRhsFn result; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhsfn", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (SUNRhsFn) ((arg1)->rhsfn); + fresult = result; + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set(SwigClassWrapper const *farg1, void *farg2) { + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + void *arg2 = (void *) 0 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhs_data", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (void *)(farg2); + if (arg1) (arg1)->rhs_data = arg2; +} + + +SWIGEXPORT void * _wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get(SwigClassWrapper const *farg1) { + void * fresult ; + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + void *result = 0 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::rhs_data", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (void *) ((arg1)->rhs_data); + fresult = result; + return fresult; +} + + +SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__nfevals_set(SwigClassWrapper const *farg1, long const *farg2) { + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + long arg2 ; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::nfevals", return ); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + arg2 = (long)(*farg2); + if (arg1) (arg1)->nfevals = arg2; +} + + +SWIGEXPORT long _wrap_SUNDomEigEstimatorContent_Arnoldi__nfevals_get(SwigClassWrapper const *farg1) { + long fresult ; + struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; + long result; + + SWIG_check_mutable_nonnull(*farg1, "struct SUNDomEigEstimatorContent_Arnoldi_ *", "SUNDomEigEstimatorContent_Arnoldi_", "SUNDomEigEstimatorContent_Arnoldi_::nfevals", return 0); + arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *)(farg1->cptr); + result = (long) ((arg1)->nfevals); + fresult = (long)(result); + return fresult; +} + + SWIGEXPORT void _wrap_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_set(SwigClassWrapper const *farg1, double *farg2) { struct SUNDomEigEstimatorContent_Arnoldi_ *arg1 = (struct SUNDomEigEstimatorContent_Arnoldi_ *) 0 ; sunrealtype *arg2 = (sunrealtype *) 0 ; @@ -723,6 +939,38 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetATimes_Arnoldi(SUNDomEigEstimator fa } +SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhs_Arnoldi(SUNDomEigEstimator farg1, void *farg2, SUNRhsFn farg3) { + int fresult ; + SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; + void *arg2 = (void *) 0 ; + SUNRhsFn arg3 = (SUNRhsFn) 0 ; + SUNErrCode result; + + arg1 = (SUNDomEigEstimator)(farg1); + arg2 = (void *)(farg2); + arg3 = (SUNRhsFn)(farg3); + result = (SUNErrCode)SUNDomEigEstimator_SetRhs_Arnoldi(arg1,arg2,arg3); + fresult = (SUNErrCode)(result); + return fresult; +} + + +SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(SUNDomEigEstimator farg1, double const *farg2, N_Vector farg3) { + int fresult ; + SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; + sunrealtype arg2 ; + N_Vector arg3 = (N_Vector) 0 ; + SUNErrCode result; + + arg1 = (SUNDomEigEstimator)(farg1); + arg2 = (sunrealtype)(*farg2); + arg3 = (N_Vector)(farg3); + result = (SUNErrCode)SUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(arg1,arg2,arg3); + fresult = (SUNErrCode)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(SUNDomEigEstimator farg1, int const *farg2) { int fresult ; SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; @@ -737,6 +985,20 @@ SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(SUNDomEig } +SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetRelTol_Arnoldi(SUNDomEigEstimator farg1, double const *farg2) { + int fresult ; + SUNDomEigEstimator arg1 = (SUNDomEigEstimator) 0 ; + sunrealtype arg2 ; + SUNErrCode result; + + arg1 = (SUNDomEigEstimator)(farg1); + arg2 = (sunrealtype)(*farg2); + result = (SUNErrCode)SUNDomEigEstimator_SetRelTol_Arnoldi(arg1,arg2); + fresult = (SUNErrCode)(result); + return fresult; +} + + SWIGEXPORT int _wrap_FSUNDomEigEstimator_SetInitialGuess_Arnoldi(SUNDomEigEstimator farg1, N_Vector 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 21ed6d482c..9ac8f69d94 100644 --- a/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.f90 +++ b/src/sundomeigest/arnoldi/fmod_int64/fsundomeigest_arnoldi_mod.f90 @@ -28,6 +28,7 @@ module fsundomeigest_arnoldi_mod private ! DECLARATION CONSTRUCTS + integer(C_INT), parameter, public :: MAX_DQITERS = 3_C_INT integer, parameter :: swig_cmem_own_bit = 0 integer, parameter :: swig_cmem_rvalue_bit = 1 @@ -48,14 +49,32 @@ module fsundomeigest_arnoldi_mod 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 + procedure :: get_Fy => swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_get + procedure :: set_work => swigf_SUNDomEigEstimatorContent_Arnoldi__work_set + procedure :: get_work => swigf_SUNDomEigEstimatorContent_Arnoldi__work_get procedure :: set_kry_dim => swigf_SUNDomEigEstimatorContent_Arnoldi__kry_dim_set procedure :: get_kry_dim => swigf_SUNDomEigEstimatorContent_Arnoldi__kry_dim_get procedure :: set_num_warmups => swigf_SUNDomEigEstimatorContent_Arnoldi__num_warmups_set procedure :: get_num_warmups => swigf_SUNDomEigEstimatorContent_Arnoldi__num_warmups_get procedure :: set_num_iters => swigf_SUNDomEigEstimatorContent_Arnoldi__num_iters_set 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_tol_preprocess => swigf_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set + procedure :: get_tol_preprocess => swigf_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get + procedure :: set_rhs_linT => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set + procedure :: get_rhs_linT => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get procedure :: set_num_ATimes => swigf_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_set procedure :: get_num_ATimes => swigf_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_get + procedure :: set_rhsfn => swigf_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set + procedure :: get_rhsfn => swigf_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get + procedure :: set_rhs_data => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set + procedure :: get_rhs_data => swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get + procedure :: set_nfevals => swigf_SUNDomEigEstimatorContent_Arnoldi__nfevals_set + procedure :: get_nfevals => swigf_SUNDomEigEstimatorContent_Arnoldi__nfevals_get procedure :: set_LAPACK_A => swigf_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_set procedure :: get_LAPACK_A => swigf_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_get procedure :: set_LAPACK_wr => swigf_SUNDomEigEstimatorContent_Arnoldi__LAPACK_wr_set @@ -79,7 +98,10 @@ module fsundomeigest_arnoldi_mod end interface public :: FSUNDomEigEstimator_Arnoldi public :: FSUNDomEigEstimator_SetATimes_Arnoldi + public :: FSUNDomEigEstimator_SetRhs_Arnoldi + public :: FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi public :: FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi + public :: FSUNDomEigEstimator_SetRelTol_Arnoldi public :: FSUNDomEigEstimator_SetInitialGuess_Arnoldi public :: FSUNDomEigEstimator_Initialize_Arnoldi public :: FSUNDomEigEstimator_Estimate_Arnoldi @@ -158,6 +180,57 @@ function swigc_SUNDomEigEstimatorContent_Arnoldi__q_get(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 +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR), value :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR), value :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__Fy_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__work_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__work_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR), value :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__work_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__work_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: fresult +end function + subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__kry_dim_set(farg1, farg2) & bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__kry_dim_set") use, intrinsic :: ISO_C_BINDING @@ -209,6 +282,57 @@ function swigc_SUNDomEigEstimatorContent_Arnoldi__num_iters_get(farg1) & integer(C_LONG) :: fresult end function +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +integer(C_INT), intent(in) :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +integer(C_INT) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE) :: fresult +end function + subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_set(farg1, farg2) & bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_set") use, intrinsic :: ISO_C_BINDING @@ -226,6 +350,57 @@ function swigc_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_get(farg1) & integer(C_LONG) :: fresult end function +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_FUNPTR), value :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_FUNPTR) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR), value :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: fresult +end function + +subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__nfevals_set(farg1, farg2) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__nfevals_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +integer(C_LONG), intent(in) :: farg2 +end subroutine + +function swigc_SUNDomEigEstimatorContent_Arnoldi__nfevals_get(farg1) & +bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__nfevals_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: farg1 +integer(C_LONG) :: fresult +end function + subroutine swigc_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_set(farg1, farg2) & bind(C, name="_wrap_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_set") use, intrinsic :: ISO_C_BINDING @@ -388,6 +563,26 @@ function swigc_FSUNDomEigEstimator_SetATimes_Arnoldi(farg1, farg2, farg3) & integer(C_INT) :: fresult end function +function swigc_FSUNDomEigEstimator_SetRhs_Arnoldi(farg1, farg2, farg3) & +bind(C, name="_wrap_FSUNDomEigEstimator_SetRhs_Arnoldi") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +type(C_PTR), value :: farg2 +type(C_FUNPTR), value :: farg3 +integer(C_INT) :: fresult +end function + +function swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(farg1, farg2, farg3) & +bind(C, name="_wrap_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +type(C_PTR), value :: farg3 +integer(C_INT) :: fresult +end function + function swigc_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(farg1, farg2) & bind(C, name="_wrap_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi") & result(fresult) @@ -397,6 +592,15 @@ function swigc_FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(farg1, farg2) & integer(C_INT) :: fresult end function +function swigc_FSUNDomEigEstimator_SetRelTol_Arnoldi(farg1, farg2) & +bind(C, name="_wrap_FSUNDomEigEstimator_SetRelTol_Arnoldi") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +type(C_PTR), value :: farg1 +real(C_DOUBLE), intent(in) :: farg2 +integer(C_INT) :: fresult +end function + function swigc_FSUNDomEigEstimator_SetInitialGuess_Arnoldi(farg1, farg2) & bind(C, name="_wrap_FSUNDomEigEstimator_SetInitialGuess_Arnoldi") & result(fresult) @@ -564,6 +768,81 @@ function swigf_SUNDomEigEstimatorContent_Arnoldi__q_get(self) & 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 +type(N_Vector), target, intent(inout) :: rhs_liny +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: farg2 + +farg1 = self%swigdata +farg2 = c_loc(rhs_liny) +call swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linY_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__rhs_linY_get(farg1) +call c_f_pointer(fresult, swig_result) +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_set(self, fy) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(N_Vector), target, intent(inout) :: fy +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: farg2 + +farg1 = self%swigdata +farg2 = c_loc(fy) +call swigc_SUNDomEigEstimatorContent_Arnoldi__Fy_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__Fy_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__Fy_get(farg1) +call c_f_pointer(fresult, swig_result) +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__work_set(self, work) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(N_Vector), target, intent(inout) :: work +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: farg2 + +farg1 = self%swigdata +farg2 = c_loc(work) +call swigc_SUNDomEigEstimatorContent_Arnoldi__work_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__work_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__work_get(farg1) +call c_f_pointer(fresult, swig_result) +end function + subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__kry_dim_set(self, kry_dim) use, intrinsic :: ISO_C_BINDING class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self @@ -639,6 +918,81 @@ function swigf_SUNDomEigEstimatorContent_Arnoldi__num_iters_get(self) & swig_result = fresult end function +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_set(self, warmup_to_tol) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +integer(C_INT), intent(in) :: warmup_to_tol +type(SwigClassWrapper) :: farg1 +integer(C_INT) :: farg2 + +farg1 = self%swigdata +farg2 = warmup_to_tol +call swigc_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__warmup_to_tol_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__warmup_to_tol_get(farg1) +swig_result = fresult +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set(self, tol_preprocess) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +real(C_DOUBLE), intent(in) :: tol_preprocess +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = self%swigdata +farg2 = tol_preprocess +call swigc_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +real(C_DOUBLE) :: swig_result +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +real(C_DOUBLE) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__tol_preprocess_get(farg1) +swig_result = fresult +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set(self, rhs_lint) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +real(C_DOUBLE), intent(in) :: rhs_lint +type(SwigClassWrapper) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = self%swigdata +farg2 = rhs_lint +call swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +real(C_DOUBLE) :: swig_result +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +real(C_DOUBLE) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_linT_get(farg1) +swig_result = fresult +end function + subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_set(self, num_atimes) use, intrinsic :: ISO_C_BINDING class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self @@ -664,6 +1018,81 @@ function swigf_SUNDomEigEstimatorContent_Arnoldi__num_ATimes_get(self) & swig_result = fresult end function +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set(self, rhsfn) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(C_FUNPTR), intent(in), value :: rhsfn +type(SwigClassWrapper) :: farg1 +type(C_FUNPTR) :: farg2 + +farg1 = self%swigdata +farg2 = rhsfn +call swigc_SUNDomEigEstimatorContent_Arnoldi__rhsfn_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +type(C_FUNPTR) :: swig_result +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(C_FUNPTR) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__rhsfn_get(farg1) +swig_result = fresult +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set(self, rhs_data) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(C_PTR) :: rhs_data +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: farg2 + +farg1 = self%swigdata +farg2 = rhs_data +call swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_data_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +type(C_PTR) :: swig_result +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +type(C_PTR) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__rhs_data_get(farg1) +swig_result = fresult +end function + +subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__nfevals_set(self, nfevals) +use, intrinsic :: ISO_C_BINDING +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +integer(C_LONG), intent(in) :: nfevals +type(SwigClassWrapper) :: farg1 +integer(C_LONG) :: farg2 + +farg1 = self%swigdata +farg2 = nfevals +call swigc_SUNDomEigEstimatorContent_Arnoldi__nfevals_set(farg1, farg2) +end subroutine + +function swigf_SUNDomEigEstimatorContent_Arnoldi__nfevals_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_LONG) :: swig_result +class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self +integer(C_LONG) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_SUNDomEigEstimatorContent_Arnoldi__nfevals_get(farg1) +swig_result = fresult +end function + subroutine swigf_SUNDomEigEstimatorContent_Arnoldi__LAPACK_A_set(self, lapack_a) use, intrinsic :: ISO_C_BINDING class(SUNDomEigEstimatorContent_Arnoldi_), intent(in) :: self @@ -914,6 +1343,44 @@ function FSUNDomEigEstimator_SetATimes_Arnoldi(dee, a_data, atimes) & swig_result = fresult end function +function FSUNDomEigEstimator_SetRhs_Arnoldi(dee, rhs_data, rhsfn) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(SUNDomEigEstimator), target, intent(inout) :: dee +type(C_PTR) :: rhs_data +type(C_FUNPTR), intent(in), value :: rhsfn +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +type(C_PTR) :: farg2 +type(C_FUNPTR) :: farg3 + +farg1 = c_loc(dee) +farg2 = rhs_data +farg3 = rhsfn +fresult = swigc_FSUNDomEigEstimator_SetRhs_Arnoldi(farg1, farg2, farg3) +swig_result = fresult +end function + +function FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(dee, t, v) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(SUNDomEigEstimator), target, intent(inout) :: dee +real(C_DOUBLE), intent(in) :: t +type(N_Vector), target, intent(inout) :: v +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 +type(C_PTR) :: farg3 + +farg1 = c_loc(dee) +farg2 = t +farg3 = c_loc(v) +fresult = swigc_FSUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi(farg1, farg2, farg3) +swig_result = fresult +end function + function FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(dee, num_iters) & result(swig_result) use, intrinsic :: ISO_C_BINDING @@ -930,6 +1397,22 @@ function FSUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(dee, num_iters) & swig_result = fresult end function +function FSUNDomEigEstimator_SetRelTol_Arnoldi(dee, tol) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(SUNDomEigEstimator), target, intent(inout) :: dee +real(C_DOUBLE), intent(in) :: tol +integer(C_INT) :: fresult +type(C_PTR) :: farg1 +real(C_DOUBLE) :: farg2 + +farg1 = c_loc(dee) +farg2 = tol +fresult = swigc_FSUNDomEigEstimator_SetRelTol_Arnoldi(farg1, farg2) +swig_result = fresult +end function + function FSUNDomEigEstimator_SetInitialGuess_Arnoldi(dee, q) & 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 47d89e8b10..1d7be20587 100644 --- a/src/sundomeigest/arnoldi/sundomeigest_arnoldi.c +++ b/src/sundomeigest/arnoldi/sundomeigest_arnoldi.c @@ -38,10 +38,12 @@ #error Incompatible sunrealtype for LAPACK; disable LAPACK and rebuild #endif -#define ONE SUN_RCONST(1.0) +#define ZERO SUN_RCONST(0.0) +#define ONE SUN_RCONST(1.0) /* Default estimator parameters */ #define DEE_NUM_OF_WARMUPS_ARNOLDI_DEFAULT 100 +#define DEE_TOL_OF_WARMUPS_ARNOLDI_DEFAULT SUN_RCONST(0.005) /* Default Arnoldi Iteration parameters */ #define DEE_KRYLOV_DIM_DEFAULT 3 @@ -62,6 +64,8 @@ int sundomeigest_Compare(const void* a, const void* b); +SUNErrCode dee_DQJtimes_Arnoldi(void* voidstarDEE, N_Vector v, N_Vector Jv); + /* * ----------------------------------------------------------------- * exported functions @@ -102,8 +106,12 @@ SUNDomEigEstimator SUNDomEigEstimator_Arnoldi(N_Vector q, int kry_dim, /* Attach operations */ DEE->ops->setatimes = SUNDomEigEstimator_SetATimes_Arnoldi; + DEE->ops->setrhs = SUNDomEigEstimator_SetRhs_Arnoldi; + DEE->ops->setrhslinearizationpoint = + SUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi; DEE->ops->setnumpreprocessiters = SUNDomEigEstimator_SetNumPreprocessIters_Arnoldi; + DEE->ops->setreltol = SUNDomEigEstimator_SetRelTol_Arnoldi; DEE->ops->setinitialguess = SUNDomEigEstimator_SetInitialGuess_Arnoldi; DEE->ops->initialize = SUNDomEigEstimator_Initialize_Arnoldi; DEE->ops->estimate = SUNDomEigEstimator_Estimate_Arnoldi; @@ -121,21 +129,30 @@ SUNDomEigEstimator SUNDomEigEstimator_Arnoldi(N_Vector q, int kry_dim, DEE->content = content; /* Fill content */ - content->ATimes = NULL; - content->ATdata = NULL; - content->V = NULL; - content->q = NULL; - content->kry_dim = kry_dim; - content->num_warmups = DEE_NUM_OF_WARMUPS_ARNOLDI_DEFAULT; - content->num_iters = 0; - content->num_ATimes = 0; - content->LAPACK_A = NULL; - content->LAPACK_wr = NULL; - content->LAPACK_wi = NULL; - content->LAPACK_work = NULL; - content->LAPACK_lwork = 0; - content->LAPACK_arr = NULL; - content->Hes = NULL; + content->ATimes = NULL; + content->ATdata = NULL; + content->V = NULL; + content->q = NULL; + content->rhs_linY = NULL; + content->rhs_linT = ZERO; + content->Fy = NULL; + content->work = NULL; + content->kry_dim = kry_dim; + content->num_warmups = DEE_NUM_OF_WARMUPS_ARNOLDI_DEFAULT; + content->num_iters = 0; + content->num_ATimes = 0; + content->warmup_to_tol = SUNFALSE; + content->tol_preprocess = DEE_TOL_OF_WARMUPS_ARNOLDI_DEFAULT; + content->rhsfn = NULL; + content->rhs_data = NULL; + content->nfevals = 0; + content->LAPACK_A = NULL; + content->LAPACK_wr = NULL; + content->LAPACK_wi = NULL; + content->LAPACK_work = NULL; + content->LAPACK_lwork = 0; + content->LAPACK_arr = NULL; + content->Hes = NULL; /* Allocate content */ content->q = N_VClone(q); @@ -171,6 +188,48 @@ SUNErrCode SUNDomEigEstimator_SetATimes_Arnoldi(SUNDomEigEstimator DEE, return SUN_SUCCESS; } +SUNErrCode SUNDomEigEstimator_SetRhs_Arnoldi(SUNDomEigEstimator DEE, + void* rhs_data, SUNRhsFn RHSfn) +{ + SUNFunctionBegin(DEE->sunctx); + + SUNAssert(DEE, SUN_ERR_ARG_CORRUPT); + SUNAssert(Arnoldi_CONTENT(DEE), SUN_ERR_ARG_CORRUPT); + + /* set function pointers to integrator-supplied RHS routine + and data, and return with success */ + Arnoldi_CONTENT(DEE)->rhsfn = RHSfn; + Arnoldi_CONTENT(DEE)->rhs_data = rhs_data; + + DEE->ops->setatimes(DEE, (void*)DEE, dee_DQJtimes_Arnoldi); + SUNCheckLastErr(); + + return SUN_SUCCESS; +} + +SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint_Arnoldi( + SUNDomEigEstimator DEE, sunrealtype t, N_Vector y) +{ + SUNFunctionBegin(DEE->sunctx); + + SUNAssert(DEE, SUN_ERR_ARG_CORRUPT); + SUNAssert(Arnoldi_CONTENT(DEE), SUN_ERR_ARG_CORRUPT); + SUNAssert(y, SUN_ERR_ARG_CORRUPT); + + if (Arnoldi_CONTENT(DEE)->rhs_linY == NULL) + { + Arnoldi_CONTENT(DEE)->rhs_linY = N_VClone(y); + SUNCheckLastErr(); + } + + Arnoldi_CONTENT(DEE)->rhs_linT = t; + + N_VScale(ONE, y, Arnoldi_CONTENT(DEE)->rhs_linY); + SUNCheckLastErr(); + + return SUN_SUCCESS; +} + SUNErrCode SUNDomEigEstimator_Initialize_Arnoldi(SUNDomEigEstimator DEE) { SUNFunctionBegin(DEE->sunctx); @@ -282,6 +341,33 @@ SUNErrCode SUNDomEigEstimator_SetNumPreprocessIters_Arnoldi(SUNDomEigEstimator D /* set the number of warmups */ Arnoldi_CONTENT(DEE)->num_warmups = num_iters; + + return SUN_SUCCESS; +} + +SUNErrCode SUNDomEigEstimator_SetRelTol_Arnoldi(SUNDomEigEstimator DEE, + sunrealtype tol) +{ + SUNFunctionBegin(DEE->sunctx); + + SUNAssert(DEE, SUN_ERR_ARG_CORRUPT); + SUNAssert(Arnoldi_CONTENT(DEE), SUN_ERR_ARG_CORRUPT); + + /* set the tolerance for preprocessing iterations */ + if (tol < ZERO) + { + Arnoldi_CONTENT(DEE)->warmup_to_tol = SUNFALSE; + return SUN_SUCCESS; + } + else if (tol == ZERO || tol > ONE - SUN_UNIT_ROUNDOFF) + { + tol = DEE_TOL_OF_WARMUPS_ARNOLDI_DEFAULT; + } + Arnoldi_CONTENT(DEE)->tol_preprocess = tol; + + /* set the type of warmup iterations */ + Arnoldi_CONTENT(DEE)->warmup_to_tol = SUNTRUE; + return SUN_SUCCESS; } @@ -327,6 +413,10 @@ SUNErrCode SUNDomEigEstimator_Estimate_Arnoldi(SUNDomEigEstimator DEE, Arnoldi_CONTENT(DEE)->num_ATimes = 0; Arnoldi_CONTENT(DEE)->num_iters = 0; + sunrealtype res; + sunrealtype newnormlambda = ZERO; + sunrealtype oldnormlambda = ZERO; + /* Set the initial q = A^{num_warmups}q/||A^{num_warmups}q|| */ for (int i = 0; i < Arnoldi_CONTENT(DEE)->num_warmups; i++) { @@ -337,12 +427,29 @@ SUNErrCode SUNDomEigEstimator_Estimate_Arnoldi(SUNDomEigEstimator DEE, Arnoldi_CONTENT(DEE)->num_iters++; if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; } + if (Arnoldi_CONTENT(DEE)->warmup_to_tol) + { + newnormlambda = N_VDotProd(Arnoldi_CONTENT(DEE)->V[0], + Arnoldi_CONTENT(DEE)->q); //Rayleigh quotient + SUNCheckLastErr(); + } + normq = N_VDotProd(Arnoldi_CONTENT(DEE)->q, Arnoldi_CONTENT(DEE)->q); SUNCheckLastErr(); normq = SUNRsqrt(normq); N_VScale(ONE / normq, Arnoldi_CONTENT(DEE)->q, Arnoldi_CONTENT(DEE)->V[0]); SUNCheckLastErr(); + + if (Arnoldi_CONTENT(DEE)->warmup_to_tol) + { + res = SUNRabs(newnormlambda - oldnormlambda); + oldnormlambda = newnormlambda; + if (res <= Arnoldi_CONTENT(DEE)->tol_preprocess * SUNRabs(newnormlambda)) + { + break; + } + } } for (int i = 0; i < n; i++) @@ -488,6 +595,21 @@ SUNErrCode SUNDomEigEstimator_Destroy_Arnoldi(SUNDomEigEstimator* DEEptr) N_VDestroy(Arnoldi_CONTENT(DEE)->q); Arnoldi_CONTENT(DEE)->q = NULL; } + if (Arnoldi_CONTENT(DEE)->rhs_linY) + { + N_VDestroy(Arnoldi_CONTENT(DEE)->rhs_linY); + Arnoldi_CONTENT(DEE)->rhs_linY = NULL; + } + if (Arnoldi_CONTENT(DEE)->Fy) + { + N_VDestroy(Arnoldi_CONTENT(DEE)->Fy); + Arnoldi_CONTENT(DEE)->Fy = NULL; + } + if (Arnoldi_CONTENT(DEE)->work) + { + N_VDestroy(Arnoldi_CONTENT(DEE)->work); + Arnoldi_CONTENT(DEE)->work = NULL; + } if (Arnoldi_CONTENT(DEE)->V) { N_VDestroyVectorArray(Arnoldi_CONTENT(DEE)->V, @@ -561,3 +683,83 @@ int sundomeigest_Compare(const void* a, const void* b) sunrealtype mag_b = SUNRsqrt(cplx_b[0] * cplx_b[0] + cplx_b[1] * cplx_b[1]); return (mag_b > mag_a) - (mag_b < mag_a); // Descending order } + +/*--------------------------------------------------------------- + dee_DQJtimes_Arnoldi: + + 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. + ---------------------------------------------------------------*/ +SUNErrCode dee_DQJtimes_Arnoldi(void* voidstarDEE, N_Vector v, N_Vector Jv) +{ + SUNDomEigEstimator DEE = (SUNDomEigEstimator)voidstarDEE; + SUNFunctionBegin(DEE->sunctx); + + SUNAssert(DEE, SUN_ERR_ARG_CORRUPT); + SUNAssert(Arnoldi_CONTENT(DEE), SUN_ERR_ARG_CORRUPT); + SUNAssert(v, SUN_ERR_ARG_CORRUPT); + SUNAssert(Jv, SUN_ERR_ARG_CORRUPT); + SUNAssert(Arnoldi_CONTENT(DEE)->rhsfn, SUN_ERR_ARG_CORRUPT); + SUNAssert(Arnoldi_CONTENT(DEE)->rhs_linY, SUN_ERR_ARG_CORRUPT); + + sunrealtype vdotv = N_VDotProd(v, v); + if (vdotv <= SUN_SMALL_REAL) + { + N_VScale(ZERO, v, Jv); + SUNCheckLastErr(); + return SUN_SUCCESS; + } + + if (Arnoldi_CONTENT(DEE)->work == NULL) + { + Arnoldi_CONTENT(DEE)->work = N_VClone(v); + SUNCheckLastErr(); + } + if (Arnoldi_CONTENT(DEE)->Fy == NULL) + { + Arnoldi_CONTENT(DEE)->Fy = N_VClone(v); + SUNCheckLastErr(); + } + + 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; + + 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; } + + /* Initialize perturbation */ + sunrealtype ydotv = N_VDotProd(y, v); + sunrealtype sq1norm = N_VL1Norm(v); + sunrealtype sign = (ydotv >= ZERO) ? ONE : -ONE; + sunrealtype sqrteps = SUNRsqrt(SUN_UNIT_ROUNDOFF); + sig = sign * sqrteps * SUNMAX(SUNRabs(ydotv), sq1norm) / vdotv; + + for (iter = 0; iter < MAX_DQITERS; iter++) + { + /* Set work = y + sig*v */ + 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++; + if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; } + + /* If f failed recoverably, shrink sig and retry */ + sig *= SUN_RCONST(0.25); + } + + /* Replace Jv by (Jv - fn)/sig */ + siginv = ONE / sig; + N_VLinearSum(siginv, Jv, -siginv, Fy, Jv); + + return SUN_SUCCESS; +} \ No newline at end of file diff --git a/src/sundomeigest/power/sundomeigest_power.c b/src/sundomeigest/power/sundomeigest_power.c index ad4b30ec47..568b63a6e5 100644 --- a/src/sundomeigest/power/sundomeigest_power.c +++ b/src/sundomeigest/power/sundomeigest_power.c @@ -51,6 +51,18 @@ * ----------------------------------------------------------------- */ +/* +* -------------------------------------------------------------------------- +* private functions +* -------------------------------------------------------------------------- +*/ + +SUNErrCode sundomeigestimator_complex_dom_eigs_from_PI( + SUNDomEigEstimator DEE, sunrealtype lambdaR, N_Vector v_prev, N_Vector v, + sunrealtype* lambdaR_out, sunrealtype* lambdaI_out); + +SUNErrCode dee_DQJtimes_Power(void* voidstarDEE, N_Vector v, N_Vector Jv); + /* ---------------------------------------------------------------------------- * Function to create a new PI estimator */ @@ -79,8 +91,11 @@ SUNDomEigEstimator SUNDomEigEstimator_Power(N_Vector q, long int max_iters, /* check for max_iters values; if illegal use defaults */ if (max_iters <= 0) { max_iters = DEE_MAX_ITER_DEFAULT; } - /* Check if rel_tol > 0 */ - if (rel_tol < SUN_SMALL_REAL) { rel_tol = DEE_TOL_DEFAULT; } + /* Check if rel_tol > 0 and < 1 */ + if (rel_tol < SUN_SMALL_REAL || rel_tol > ONE - SUN_UNIT_ROUNDOFF) + { + rel_tol = DEE_TOL_DEFAULT; + } /* Create dominant eigenvalue estimator */ DEE = NULL; @@ -88,7 +103,10 @@ SUNDomEigEstimator SUNDomEigEstimator_Power(N_Vector q, long int max_iters, SUNCheckLastErrNull(); /* Attach operations */ - DEE->ops->setatimes = SUNDomEigEstimator_SetATimes_Power; + DEE->ops->setatimes = SUNDomEigEstimator_SetATimes_Power; + DEE->ops->setrhs = SUNDomEigEstimator_SetRhs_Power; + DEE->ops->setrhslinearizationpoint = + SUNDomEigEstimator_SetRhsLinearizationPoint_Power; DEE->ops->setmaxiters = SUNDomEigEstimator_SetMaxIters_Power; DEE->ops->setnumpreprocessiters = SUNDomEigEstimator_SetNumPreprocessIters_Power; DEE->ops->setreltol = SUNDomEigEstimator_SetRelTol_Power; @@ -114,10 +132,19 @@ SUNDomEigEstimator SUNDomEigEstimator_Power(N_Vector q, long int max_iters, 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->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; @@ -156,6 +183,72 @@ SUNErrCode SUNDomEigEstimator_SetATimes_Power(SUNDomEigEstimator DEE, return SUN_SUCCESS; } +SUNErrCode SUNDomEigEstimator_SetRhs_Power(SUNDomEigEstimator DEE, + void* rhs_data, SUNRhsFn RHSfn) +{ + SUNFunctionBegin(DEE->sunctx); + + SUNAssert(DEE, SUN_ERR_ARG_CORRUPT); + SUNAssert(PI_CONTENT(DEE), SUN_ERR_ARG_CORRUPT); + SUNAssert(RHSfn, SUN_ERR_ARG_CORRUPT); + + /* set function pointers to integrator-supplied RHS routine + and data, and return with success */ + PI_CONTENT(DEE)->rhsfn = RHSfn; + PI_CONTENT(DEE)->rhs_data = rhs_data; + + DEE->ops->setatimes(DEE, (void*)DEE, dee_DQJtimes_Power); + SUNCheckLastErr(); + + return SUN_SUCCESS; +} + +SUNErrCode SUNDomEigEstimator_SetRhsLinearizationPoint_Power( + SUNDomEigEstimator DEE, sunrealtype t, N_Vector y) +{ + SUNFunctionBegin(DEE->sunctx); + + SUNAssert(DEE, SUN_ERR_ARG_CORRUPT); + SUNAssert(PI_CONTENT(DEE), SUN_ERR_ARG_CORRUPT); + SUNAssert(y, SUN_ERR_ARG_CORRUPT); + + if (PI_CONTENT(DEE)->rhs_linY == NULL) + { + PI_CONTENT(DEE)->rhs_linY = N_VClone(y); + SUNCheckLastErr(); + } + + PI_CONTENT(DEE)->rhs_linT = t; + + N_VScale(ONE, y, PI_CONTENT(DEE)->rhs_linY); + SUNCheckLastErr(); + + return SUN_SUCCESS; +} + +SUNErrCode SUNDomEigEstimator_SetIsReal_Power(SUNDomEigEstimator DEE, + sunbooleantype real) +{ + SUNFunctionBegin(DEE->sunctx); + + SUNAssert(DEE, SUN_ERR_ARG_CORRUPT); + SUNAssert(PI_CONTENT(DEE), SUN_ERR_ARG_CORRUPT); + + /* set the complex flag to the opposite of the real flag */ + PI_CONTENT(DEE)->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)->complex) && PI_CONTENT(DEE)->q_prev) + { + N_VDestroy(PI_CONTENT(DEE)->q_prev); + PI_CONTENT(DEE)->q_prev = NULL; + } + + return SUN_SUCCESS; +} + SUNErrCode SUNDomEigEstimator_Initialize_Power(SUNDomEigEstimator DEE) { SUNFunctionBegin(DEE->sunctx); @@ -163,7 +256,8 @@ SUNErrCode SUNDomEigEstimator_Initialize_Power(SUNDomEigEstimator DEE) SUNAssert(DEE, SUN_ERR_ARG_CORRUPT); SUNAssert(PI_CONTENT(DEE), SUN_ERR_ARG_CORRUPT); - if (PI_CONTENT(DEE)->rel_tol < SUN_SMALL_REAL) + if (PI_CONTENT(DEE)->rel_tol < SUN_SMALL_REAL || + PI_CONTENT(DEE)->rel_tol > ONE - SUN_UNIT_ROUNDOFF) { PI_CONTENT(DEE)->rel_tol = DEE_TOL_DEFAULT; } @@ -180,6 +274,11 @@ SUNErrCode SUNDomEigEstimator_Initialize_Power(SUNDomEigEstimator DEE) SUNAssert(PI_CONTENT(DEE)->V, SUN_ERR_ARG_CORRUPT); SUNAssert(PI_CONTENT(DEE)->q, SUN_ERR_ARG_CORRUPT); + if (PI_CONTENT(DEE)->complex) + { + SUNAssert(PI_CONTENT(DEE)->q_prev == NULL, SUN_ERR_ARG_CORRUPT); + } + /* Initialize the vector V */ sunrealtype normq = N_VDotProd(PI_CONTENT(DEE)->q, PI_CONTENT(DEE)->q); SUNCheckLastErr(); @@ -216,8 +315,11 @@ SUNErrCode SUNDomEigEstimator_SetRelTol_Power(SUNDomEigEstimator DEE, SUNAssert(DEE, SUN_ERR_ARG_CORRUPT); SUNAssert(PI_CONTENT(DEE), SUN_ERR_ARG_CORRUPT); - /* Check if rel_tol > 0 */ - if (rel_tol < SUN_SMALL_REAL) { rel_tol = DEE_TOL_DEFAULT; } + /* Check if rel_tol > 0 and < 1 */ + if (rel_tol < SUN_SMALL_REAL || rel_tol > ONE - SUN_UNIT_ROUNDOFF) + { + rel_tol = DEE_TOL_DEFAULT; + } /* set the tolerance */ PI_CONTENT(DEE)->rel_tol = rel_tol; @@ -275,6 +377,12 @@ SUNErrCode SUNDomEigEstimator_Estimate_Power(SUNDomEigEstimator DEE, 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)->complex && (PI_CONTENT(DEE)->q_prev == NULL)) + { + /* allocate q_prev vector */ + PI_CONTENT(DEE)->q_prev = N_VClone(PI_CONTENT(DEE)->q); + SUNCheckLastErr(); + } sunrealtype newlambdaR = ZERO; sunrealtype oldlambdaR = ZERO; @@ -284,6 +392,7 @@ SUNErrCode SUNDomEigEstimator_Estimate_Power(SUNDomEigEstimator DEE, PI_CONTENT(DEE)->num_ATimes = 0; PI_CONTENT(DEE)->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++) { retval = PI_CONTENT(DEE)->ATimes(PI_CONTENT(DEE)->ATdata, @@ -302,6 +411,12 @@ SUNErrCode SUNDomEigEstimator_Estimate_Power(SUNDomEigEstimator DEE, for (int k = 0; k < PI_CONTENT(DEE)->max_iters; k++) { + if (PI_CONTENT(DEE)->complex) + { + N_VScale(ONE, PI_CONTENT(DEE)->V, PI_CONTENT(DEE)->q_prev); + SUNCheckLastErr(); + } + retval = PI_CONTENT(DEE)->ATimes(PI_CONTENT(DEE)->ATdata, PI_CONTENT(DEE)->V, PI_CONTENT(DEE)->q); PI_CONTENT(DEE)->num_ATimes++; @@ -312,10 +427,6 @@ SUNErrCode SUNDomEigEstimator_Estimate_Power(SUNDomEigEstimator DEE, PI_CONTENT(DEE)->q); //Rayleigh quotient SUNCheckLastErr(); - PI_CONTENT(DEE)->res = SUNRabs(newlambdaR - oldlambdaR) / SUNRabs(newlambdaR); - - if (PI_CONTENT(DEE)->res < PI_CONTENT(DEE)->rel_tol) { break; } - normq = N_VDotProd(PI_CONTENT(DEE)->q, PI_CONTENT(DEE)->q); SUNCheckLastErr(); @@ -323,11 +434,28 @@ SUNErrCode SUNDomEigEstimator_Estimate_Power(SUNDomEigEstimator DEE, N_VScale(ONE / normq, PI_CONTENT(DEE)->q, PI_CONTENT(DEE)->V); SUNCheckLastErr(); + PI_CONTENT(DEE)->res = SUNRabs(newlambdaR - oldlambdaR); + if (PI_CONTENT(DEE)->res <= PI_CONTENT(DEE)->rel_tol * SUNRabs(newlambdaR)) + { + break; + } + oldlambdaR = newlambdaR; } - *lambdaI = ZERO; - *lambdaR = newlambdaR; + if (PI_CONTENT(DEE)->complex) + { + retval = sundomeigestimator_complex_dom_eigs_from_PI(DEE, newlambdaR, + PI_CONTENT(DEE)->q_prev, + PI_CONTENT(DEE)->V, + lambdaR, lambdaI); + if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; } + } + else + { + *lambdaR = newlambdaR; + *lambdaI = ZERO; + } return SUN_SUCCESS; } @@ -401,6 +529,92 @@ SUNErrCode SUNDomEigEstimator_Write_Power(SUNDomEigEstimator DEE, FILE* outfile) return SUN_SUCCESS; } +SUNErrCode sundomeigestimator_complex_dom_eigs_from_PI( + SUNDomEigEstimator DEE, sunrealtype lambdaR, N_Vector v_prev, N_Vector v, + sunrealtype* lambdaR_out, sunrealtype* lambdaI_out) +{ + SUNFunctionBegin(DEE->sunctx); + + int retval; + sunrealtype cos_qs, gram_det, det_G_inv, h11, h12, h21, 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 */ + sunrealtype gram_det_tol = SUNMAX(SUN_RCONST(10.0) * SUN_UNIT_ROUNDOFF, + SUN_RCONST(10.0) * PI_CONTENT(DEE)->rel_tol); + cos_qs = N_VDotProd(v_prev, v); + SUNCheckLastErr(); + + /* Safety against roundoff in dot product */ + if (cos_qs > ONE) { cos_qs = ONE; } + if (cos_qs < -ONE) { cos_qs = -ONE; } + + /* Use Gram determinant as the near-dependence measure: + G = [ [1, cos_qs], [cos_qs, 1] ], det(G) = 1 - cos_qs^2 + This assumes v_prev and v are normalized. */ + gram_det = ONE - cos_qs * cos_qs; + + if (gram_det <= gram_det_tol) + { + /* Dominant eigenvalue is real */ + *lambdaR_out = lambdaR; + *lambdaI_out = ZERO; + return SUN_SUCCESS; + } + else + { + det_G_inv = ONE / gram_det; + + /* Solve for G = [v_prev v]' * [v_prev v] and compute + projected matrix P = G^{-1} * [v_prev v]' * A * [v_prev v] */ + + retval = PI_CONTENT(DEE)->ATimes(PI_CONTENT(DEE)->ATdata, v_prev, + PI_CONTENT(DEE)->q); + PI_CONTENT(DEE)->num_ATimes++; + PI_CONTENT(DEE)->num_iters++; + if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; } + + h11 = N_VDotProd(v_prev, PI_CONTENT(DEE)->q); + SUNCheckLastErr(); + h21 = N_VDotProd(v, PI_CONTENT(DEE)->q); + SUNCheckLastErr(); + + retval = PI_CONTENT(DEE)->ATimes(PI_CONTENT(DEE)->ATdata, v, + PI_CONTENT(DEE)->q); + PI_CONTENT(DEE)->num_ATimes++; + PI_CONTENT(DEE)->num_iters++; + if (retval != 0) { return SUN_ERR_USER_FCN_FAIL; } + + h12 = N_VDotProd(v_prev, PI_CONTENT(DEE)->q); + SUNCheckLastErr(); + h22 = N_VDotProd(v, PI_CONTENT(DEE)->q); + SUNCheckLastErr(); + + p11 = det_G_inv * (h11 - cos_qs * h21); + p12 = det_G_inv * (h12 - cos_qs * h22); + p21 = det_G_inv * (h21 - cos_qs * h11); + p22 = det_G_inv * (h22 - cos_qs * h12); + + /* Compute eigenvalues of P */ + sunrealtype traceP = p11 + p22; + sunrealtype detP = p11 * p22 - p12 * p21; + sunrealtype discrim = traceP * traceP - SUN_RCONST(4.0) * detP; + if (discrim >= ZERO) + { + /* Dominant eigenvalue is real */ + *lambdaR_out = (traceP + SUNRsqrt(discrim)) / SUN_RCONST(2.0); + *lambdaI_out = ZERO; + } + else + { + /* Dominant eigenvalue is complex */ + *lambdaR_out = traceP / SUN_RCONST(2.0); + *lambdaI_out = SUNRsqrt(-discrim) / SUN_RCONST(2.0); + } + } + + return SUN_SUCCESS; +} + SUNErrCode SUNDomEigEstimator_Destroy_Power(SUNDomEigEstimator* DEEptr) { SUNFunctionBegin((*DEEptr)->sunctx); @@ -417,12 +631,31 @@ SUNErrCode SUNDomEigEstimator_Destroy_Power(SUNDomEigEstimator* DEEptr) N_VDestroy(PI_CONTENT(DEE)->q); PI_CONTENT(DEE)->q = NULL; } + if (PI_CONTENT(DEE)->q_prev) + { + N_VDestroy(PI_CONTENT(DEE)->q_prev); + PI_CONTENT(DEE)->q_prev = NULL; + } if (PI_CONTENT(DEE)->V) { N_VDestroy(PI_CONTENT(DEE)->V); PI_CONTENT(DEE)->V = NULL; } - + if (PI_CONTENT(DEE)->rhs_linY) + { + N_VDestroy(PI_CONTENT(DEE)->rhs_linY); + PI_CONTENT(DEE)->rhs_linY = NULL; + } + if (PI_CONTENT(DEE)->Fy) + { + N_VDestroy(PI_CONTENT(DEE)->Fy); + PI_CONTENT(DEE)->Fy = NULL; + } + if (PI_CONTENT(DEE)->work) + { + N_VDestroy(PI_CONTENT(DEE)->work); + PI_CONTENT(DEE)->work = NULL; + } free(DEE->content); DEE->content = NULL; } @@ -435,3 +668,84 @@ SUNErrCode SUNDomEigEstimator_Destroy_Power(SUNDomEigEstimator* DEEptr) *DEEptr = NULL; return SUN_SUCCESS; } + +/*--------------------------------------------------------------- + dee_DQJtimes_Power: + + 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. + ---------------------------------------------------------------*/ +SUNErrCode dee_DQJtimes_Power(void* voidstarDEE, N_Vector v, N_Vector Jv) +{ + SUNDomEigEstimator DEE = (SUNDomEigEstimator)voidstarDEE; + SUNFunctionBegin(DEE->sunctx); + + SUNAssert(DEE, SUN_ERR_ARG_CORRUPT); + SUNAssert(PI_CONTENT(DEE), SUN_ERR_ARG_CORRUPT); + SUNAssert(v, SUN_ERR_ARG_CORRUPT); + SUNAssert(Jv, SUN_ERR_ARG_CORRUPT); + SUNAssert(PI_CONTENT(DEE)->rhsfn, SUN_ERR_ARG_CORRUPT); + SUNAssert(PI_CONTENT(DEE)->rhs_linY, SUN_ERR_ARG_CORRUPT); + + sunrealtype vdotv = N_VDotProd(v, v); + if (vdotv <= SUN_SMALL_REAL) + { + N_VScale(ZERO, v, Jv); + SUNCheckLastErr(); + return SUN_SUCCESS; + } + + if (PI_CONTENT(DEE)->work == NULL) + { + PI_CONTENT(DEE)->work = N_VClone(v); + SUNCheckLastErr(); + } + if (PI_CONTENT(DEE)->Fy == NULL) + { + PI_CONTENT(DEE)->Fy = N_VClone(v); + SUNCheckLastErr(); + } + + 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; + + 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; } + + /* Initialize perturbation */ + sunrealtype ydotv = N_VDotProd(y, v); + sunrealtype sq1norm = N_VL1Norm(v); + sunrealtype sign = (ydotv >= ZERO) ? ONE : -ONE; + sunrealtype sqrteps = SUNRsqrt(SUN_UNIT_ROUNDOFF); + sig = sign * sqrteps * SUNMAX(SUNRabs(ydotv), sq1norm) / vdotv; + + for (iter = 0; iter < MAX_DQITERS; iter++) + { + /* Set work = y + sig*v */ + 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++; + if (retval == 0) { break; } + if (retval < 0) { return (-1); } + + /* If f failed recoverably, shrink sig and retry */ + sig *= SUN_RCONST(0.25); + } + + /* Replace Jv by (Jv - fn)/sig */ + siginv = ONE / sig; + N_VLinearSum(siginv, Jv, -siginv, Fy, Jv); + + return SUN_SUCCESS; +} diff --git a/test/answers b/test/answers index 49ef12abaf..1de5f12bd0 160000 --- a/test/answers +++ b/test/answers @@ -1 +1 @@ -Subproject commit 49ef12abaf84ae79a8d6f6c1abea53a06e0e9771 +Subproject commit 1de5f12bd0444d044cbe30d308af0ec348e45f41 diff --git a/test/unit_tests/sundomeigest/Arnoldi/test_sundomeigest_arnoldi.c b/test/unit_tests/sundomeigest/Arnoldi/test_sundomeigest_arnoldi.c index 4631632648..3cbb48dbdb 100644 --- a/test/unit_tests/sundomeigest/Arnoldi/test_sundomeigest_arnoldi.c +++ b/test/unit_tests/sundomeigest/Arnoldi/test_sundomeigest_arnoldi.c @@ -220,7 +220,7 @@ int main(int argc, char* argv[]) rel_error /= norm_of_dom_eig; - if (rel_error < rel_tol) + if (rel_error < SUN_RCONST(10.0) * rel_tol) { printf("\n\nPASS: relative error = " SUN_FORMAT_G " \n\n", rel_error); } diff --git a/test/unit_tests/sundomeigest/Power/CMakeLists.txt b/test/unit_tests/sundomeigest/Power/CMakeLists.txt index fe89578067..fb93caf1fa 100644 --- a/test/unit_tests/sundomeigest/Power/CMakeLists.txt +++ b/test/unit_tests/sundomeigest/Power/CMakeLists.txt @@ -25,7 +25,11 @@ set(sundomeigest_power_examples "test_sundomeigest_power\;100 100 0 0\;" "test_sundomeigest_power\;1000 100 0 0\;" "test_sundomeigest_power\;10000 100 0 0\;" - "test_sundomeigest_power\;100000 100 0 0\;") + "test_sundomeigest_power\;100000 100 0 0\;" + "test_sundomeigest_power_complex\;100 100 0 0\;" + "test_sundomeigest_power_complex\;1000 100 0 0\;" + "test_sundomeigest_power_complex\;10000 100 100 0\;" + "test_sundomeigest_power_complex\;100000 100 100 0\;") # Dependencies for dominant eigenvalue examples set(sundomeigest_power_dependencies test_sundomeigest) diff --git a/test/unit_tests/sundomeigest/Power/test_sundomeigest_power_complex.c b/test/unit_tests/sundomeigest/Power/test_sundomeigest_power_complex.c new file mode 100644 index 0000000000..c0eab6ed86 --- /dev/null +++ b/test/unit_tests/sundomeigest/Power/test_sundomeigest_power_complex.c @@ -0,0 +1,298 @@ +/* ----------------------------------------------------------------- + * Programmer(s): Mustafa Aggul @ SMU + * ----------------------------------------------------------------- + * SUNDIALS Copyright Start + * Copyright (c) 2025, Lawrence Livermore National Security, + * University of Maryland Baltimore County, and the SUNDIALS contributors. + * Copyright (c) 2013-2025, Lawrence Livermore National Security + * and Southern Methodist University. + * Copyright (c) 2002-2013, Lawrence Livermore National Security. + * All rights reserved. + * + * See the top-level LICENSE and NOTICE files for details. + * + * SPDX-License-Identifier: BSD-3-Clause + * SUNDIALS Copyright End + * ----------------------------------------------------------------- + * These test functions check some components of Power Iteration + * module implementation. + * ----------------------------------------------------------------- + */ + +#include +#include +#include "../test_sundomeigest.h" + +/* constants */ +#define ZERO SUN_RCONST(0.0) + +#define factor SUN_RCONST(-10.0) +#define realpart SUN_RCONST(-30000.0) +#define imagpart SUN_RCONST(+40000.0) + +/* user data structure */ +typedef struct +{ + sunindextype N; /* problem size */ + N_Vector diag; /* matrix diagonal */ + + /* nondiagonal entries of the matrix that lead to the complex conjugate eigenvalues */ + sunrealtype real_part; + sunrealtype imag_part; +} UserData; + +/* private functions */ +/* matrix-vector product */ +int ATimes(void* ProbData, N_Vector v, N_Vector z); +/* checks function return values */ +int check_flag(void* flagvalue, const char* funcname, int opt); + +/* ---------------------------------------------------------------------- + * DomEig Module Testing Routine + * --------------------------------------------------------------------*/ +int main(int argc, char* argv[]) +{ + int fails = 0; /* counter for test failures */ + int passfail = 0; /* overall pass/fail flag */ + SUNDomEigEstimator DEE = NULL; /* domeig estimator object */ + UserData ProbData; /* problem data structure */ + int num_warmups; /* number of preprocessing iters */ + long int max_iters; /* max power iteration */ + long int num_iters; /* cur. number of iterations */ + long int num_ATimes; /* number of ATimes calls */ + int print_timing; /* timing output flag */ + sunrealtype res; /* current residual */ + sunrealtype lambdaR, lambdaI; /* computed domeig parts */ + sunrealtype tlambdaR, tlambdaI; /* true domeig parts */ + SUNContext sunctx; + sunrealtype rel_tol = SUN_RCONST(1.0e-2); /* relative tol for pass/fail */ + sunrealtype rel_error; + N_Vector q; /* random initial eigenvector */ + + if (SUNContext_Create(SUN_COMM_NULL, &sunctx)) + { + printf("ERROR: SUNContext_Create failed\n"); + return (-1); + } + + /* check inputs: local problem size, max iters, num preprocessing, timing flag */ + if (argc < 5) + { + printf("ERROR: FOUR (4) Inputs required:\n"); + printf(" Problem size should be >= 2\n"); + printf(" Maximum number of power iterations should be > 0\n"); + printf(" Number of preprocessing iters should be >= 0\n"); + printf(" Include timers for calculation (0=off, 1=on)\n"); + return 1; + } + ProbData.N = (sunindextype)atol(argv[1]); + if (ProbData.N <= 0) + { + printf("ERROR: Problem size must be a positive integer\n"); + return 1; + } + max_iters = atoi(argv[2]); + if (max_iters <= 0) + { + printf( + "ERROR: Maximum number of power iterations must be a positive integer\n"); + return 1; + } + num_warmups = atoi(argv[3]); + if (num_warmups < 0) + { + printf("ERROR: Number of preprocessing must be a nonnegative integer\n"); + return 1; + } + print_timing = atoi(argv[4]); + SetTiming(print_timing); + + printf("\nDomEig module test:\n"); + printf(" Problem size = %ld\n", (long int)ProbData.N); + printf(" Number of power iterations = %ld\n", (long int)max_iters); + printf(" Number of preprocessing iters = %i\n", num_warmups); + printf(" Timing output flag = %i\n\n", print_timing); + + /* Create vectors */ + ProbData.diag = N_VNew_Serial(ProbData.N, sunctx); + if (check_flag(ProbData.diag, "N_VNew_Serial", 0)) { return 1; } + + q = N_VClone(ProbData.diag); + if (check_flag(q, "N_VClone", 0)) { return 1; } + + sunrealtype* qd = N_VGetArrayPointer(q); + for (int i = 0; i < ProbData.N; i++) + { + qd[i] = (sunrealtype)rand() / (sunrealtype)RAND_MAX; + } + + /* Fill matrix diagonal and problem data */ + // real diag is [3 4 5 ... N 0 0]*factor + // 2x2 block matrix attached to the last two diagonals is + // [ realpart imagpart; + // [-imagpart realpart] + // This setup allows two types of dominant eigenvalues (real and complex) + // based on the "factor" and the problem dimension N. + sunrealtype* v = N_VGetArrayPointer(ProbData.diag); + for (int i = 0; i < ProbData.N - 2; i++) { v[i] = factor * (i + 3); } + + // Set the problem data corresponding to 2x2 block matrix + ProbData.real_part = realpart; + ProbData.imag_part = imagpart; + + /* Create Power Iteration Dominant Eigvalue Estimator (DEE)*/ + DEE = SUNDomEigEstimator_Power(q, max_iters, rel_tol, sunctx); + if (check_flag(DEE, "SUNDomEigEstimator_Power", 0)) { return 1; } + + fails += Test_SUNDomEigEstimator_SetATimes(DEE, &ProbData, ATimes, 0); + fails += Test_SUNDomEigEstimator_SetMaxIters(DEE, max_iters, 0); + fails += Test_SUNDomEigEstimator_SetNumPreprocessIters(DEE, num_warmups, 0); + fails += Test_SUNDomEigEstimator_SetRelTol(DEE, rel_tol, 0); + fails += Test_SUNDomEigEstimator_SetInitialGuess(DEE, q, 0); + fails += Test_SUNDomEigEstimator_Initialize(DEE, 0); + fails += Test_SUNDomEigEstimator_Estimate(DEE, &lambdaR, &lambdaI, 0); + fails += Test_SUNDomEigEstimator_GetRes(DEE, &res, 0); + if (res < SUN_SMALL_REAL) + { + printf(" >>> FAILED test -- SUNDomEigEstimator_GetRes return value\n"); + fails++; + } + fails += Test_SUNDomEigEstimator_GetNumIters(DEE, &num_iters, 0); + if (num_iters <= 0) + { + printf( + " >>> FAILED test -- SUNDomEigEstimator_GetNumIters return value\n"); + fails++; + } + fails += Test_SUNDomEigEstimator_GetNumATimesCalls(DEE, &num_ATimes, 0); + fails += Test_SUNDomEigEstimator_Write(DEE, 0); + + if (fails) + { + printf("FAIL: SUNDomEigEstimator_Power module failed %i initialization " + "tests\n\n", + fails); + return 1; + } + else + { + printf("SUCCESS: SUNDomEigEstimator_Power module passed all initialization " + "tests\n\n"); + } + + /* First check if the computed eigenvalue has a nonzero magnitute */ + sunrealtype norm_of_dom_eig = SUNRsqrt(lambdaR * lambdaR + lambdaI * lambdaI); + if (norm_of_dom_eig < SUN_SMALL_REAL) + { + printf("FAIL: Dominant Eigenvalue Test Failed\n\n"); + return 1; + } + + /* Identify the tlambdaR and tlambdaI based on given parameters*/ + if (SUNRsqrt(realpart * realpart + imagpart * imagpart) > -factor * ProbData.N) + { + /* Dominant eigenvalue corresponds to the 2x2 block matrix */ + tlambdaR = realpart; + tlambdaI = imagpart; + } + else + { + /* Dominant eigenvalue corresponds to the maximum real value at the diagonal */ + tlambdaR = factor * ProbData.N; + tlambdaI = ZERO; + } + + printf("\ncomputed dominant eigenvalue = " SUN_FORMAT_G " + " SUN_FORMAT_G + " i\n", + lambdaR, lambdaI); + printf(" true dominant eigenvalue = " SUN_FORMAT_G " + " SUN_FORMAT_G + " i\n", + tlambdaR, tlambdaI); + + /* Compare the estimated dom_eig with the tlambdaR and tlambdaI*/ + rel_error = SUNRsqrt((lambdaR - tlambdaR) * (lambdaR - tlambdaR) + + (lambdaI - tlambdaI) * (lambdaI - tlambdaI)); + + rel_error /= norm_of_dom_eig; + + if (rel_error < SUN_RCONST(10.0) * rel_tol) + { + printf("\n\nPASS: relative error = " SUN_FORMAT_G " \n\n", rel_error); + } + else + { + printf("\n\nFAIL: relative error = " SUN_FORMAT_G " \n\n", rel_error); + passfail += 1; + } + + /* Free solver and vectors */ + N_VDestroy(ProbData.diag); + SUNContext_Free(&sunctx); + N_VDestroy(q); + SUNDomEigEstimator_Destroy(&DEE); + + return (passfail); +} + +/* ---------------------------------------------------------------------- + * Private helper functions + * --------------------------------------------------------------------*/ + +/* matrix-vector product */ +int ATimes(void* Data, N_Vector v_vec, N_Vector z_vec) +{ + /* local variables */ + sunrealtype *v, *z, *diag, real_part, imag_part; + sunindextype i, N; + UserData* ProbData; + + /* access user data structure and vector data */ + ProbData = (UserData*)Data; + v = N_VGetArrayPointer(v_vec); + if (check_flag(v, "N_VGetArrayPointer", 0)) { return 1; } + z = N_VGetArrayPointer(z_vec); + if (check_flag(z, "N_VGetArrayPointer", 0)) { return 1; } + N = ProbData->N; + real_part = ProbData->real_part; + imag_part = ProbData->imag_part; + diag = N_VGetArrayPointer(ProbData->diag); + if (check_flag(diag, "N_VGetArrayPointer", 0)) { return 1; } + + /* perform product on the diagonal part of the matrix */ + for (i = 0; i < N - 2; i++) { z[i] = diag[i] * v[i]; } + + /* perform product at the non-diagonal last two rows */ + z[N - 2] = v[N - 2] * real_part + v[N - 1] * imag_part; + z[N - 1] = v[N - 1] * real_part - v[N - 2] * imag_part; + /* return with success */ + return 0; +} + +/* Check function return value based on "opt" input: + 0: function allocates memory so check for NULL pointer + 1: function returns a flag so check for flag != 0 */ +int check_flag(void* flagvalue, const char* funcname, int opt) +{ + int* errflag; + + /* Check if function returned NULL pointer - no memory allocated */ + if (opt == 0 && flagvalue == NULL) + { + fprintf(stderr, "\nERROR: %s() failed - returned NULL pointer\n\n", funcname); + return 1; + } + + /* Check if flag != 0 */ + if (opt == 1) + { + errflag = (int*)flagvalue; + if (*errflag != 0) + { + fprintf(stderr, "\nERROR: %s() failed with flag = %d\n\n", funcname, + *errflag); + return 1; + } + } + + return 0; +}