12.1. Introduction to Dominant Eigenvalue Estimators

For problems that require the dominant eigenvalue of a matrix (i.e., the Jacobian), the SUNDIALS packages operate using generic dominant eigenvalue estimator modules defined through the SUNDomEigEstimator class. This allows SUNDIALS packages to utilize any valid SUNDomEigEstimator implementation that provides a set of required functions. These functions can be divided into three categories. The first are the core estimator functions. The second group consists of “set” routines to supply the dominant eigenvalue estimator object with functions provided by the SUNDIALS package, or for modification of estimator parameters. The last group consists of “get” routines for retrieving artifacts (statistics, residual, etc.) from the estimator. All of these functions are defined in the header file sundials/sundials_domeigestimator.h.

The implementations provided with SUNDIALS work in coordination with the SUNDIALS N_Vector modules to provide a set of compatible data structures for the estimator. Moreover, advanced users can provide a customized SUNDomEigEstimator implementation to any SUNDIALS package, particularly in cases where they provide their own 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, \(\tau\),

(12.1)\[\frac{\left|\lambda_k - \lambda_{k-1}\right|}{\left|\lambda_k \right|} < \tau.\]

For users interested in providing their own SUNDomEigEstimator(), the following section presents the SUNDomEigEstimator class and its implementation beginning with the definition of SUNDomEigEstimator functions in §12.2.1 – §12.2.3. This is followed by the definition of functions supplied to an estimator implementation in §12.2.4. The SUNDomEigEstimator type is defined §12.2.5. The section that then follows describes the SUNDomEigEstimator functions required by this SUNDIALS package, and provides additional package specific details. Then the remaining sections of this chapter present the SUNDomEigEstimator modules provided with SUNDIALS.

12.2. The SUNDomEigEstimator API

Added in version 7.5.0.

The SUNDomEigEst API defines several dominant eigenvalue estimation operations that enable SUNDIALS packages to utilize this API. These functions can be divided into three categories. The first are the core dominant eigenvalue estimation functions. The second consist of “set” routines to supply the dominant eigenvalue estimator with functions provided by the SUNDIALS packages and to modify estimator parameters. The final group consists of “get” routines for retrieving dominant eigenvalue estimation statistics. All of these functions are defined in the header file sundials/sundials_domeigestimator.h.

12.2.1. SUNDomEigEstimator core functions

The SUNDomEigEstimator base class provides two utility routines for implementers, SUNDomEigEstimator_NewEmpty() and SUNDomEigEstimator_FreeEmpty().

Implementations of SUNDomEigEstimators must include a required SUNDomEigEstimator_Estimate() function to estimate the dominant eigenvalue.

SUNDomEigEstimator SUNDomEigEstimator_NewEmpty(SUNContext sunctx)

This function allocates a new SUNDomEigEstimator object and initializes its content pointer and the function pointers in the operations structure to NULL.

Arguments:

• sunctx – the SUNContext object (see §1.3).

Return value:

If successful, this function returns a SUNDomEigEstimator object. If an error occurs when allocating the object, then this routine will return NULL.

SUNErrCode SUNDomEigEstimator_Initialize(SUNDomEigEstimator DEE)

This optional function performs dominant eigenvalue estimator initialization (assuming that all estimator-specific options have been set).

Arguments:

• DEE – a SUNDomEigEstimator object.

Return value:

A SUNErrCode.

SUNErrCode SUNDomEigEstimator_Estimate(SUNDomEigEstimator DEE, sunrealtype *lambdaR, sunrealtype *lambdaI)

This required function estimates the dominant eigenvalue, \(\lambda_{\max} = \lambda_{R} + \lambda_{I}i\) such that \(|\lambda| = \max\{|\lambda_i| : A \vec{v_i} = \lambda_i \vec{v_i}, \ \vec{v_i} \neq \vec{0} \}\).

Arguments:

• DEE – a SUNDomEigEstimator object.

• lambdaR – The real part of the dominant eigenvalue.

• lambdaI – The imaginary part of the dominant eigenvalue.

Return value:

SUN_SUCCESS for a successful call, or a relevant error code from SUNErrCode upon failure.

Note

When the estimator is used in a time-dependent context, an implementation may reuse the same initial guess as the initial call to SUNDomEigEstimator_Estimate() or use an improved guess based on the result of the most recent SUNDomEigEstimator_Estimate() call. See the documentation of the specific SUNDomEigEstimator implementation for more details.

SUNErrCode SUNDomEigEstimator_FreeEmpty(SUNDomEigEstimator DEE)

This routine frees the SUNDomEigEstimator object, under the assumption that any implementation-specific data that was allocated within the underlying content structure has already been freed. It will additionally test whether the ops pointer is NULL, and, if it is not, it will free it as well.

Arguments:

• DEE – a SUNDomEigEstimator object.

Return value:

A SUNErrCode.

SUNErrCode SUNDomEigEstimator_Destroy(SUNDomEigEstimator *DEEptr)

Frees memory allocated by the dominant eigenvalue estimator.

Arguments:

• DEEptr – a SUNDomEigEstimator object pointer.

Return value:

A SUNErrCode.

12.2.2. SUNDomEigEstimator “set” functions

The following functions supply dominant eigenvalue estimator modules with functions defined by the SUNDIALS packages and modify estimator parameters. When using the matrix-vector product routine provided by a SUNDIALS integration, the SetATimes is required. Otherwise, all set functions are optional. SUNDomEigEst implementations that do not provide the functionality for any optional routine should leave the corresponding function pointer NULL instead of supplying a dummy routine.

SUNErrCode SUNDomEigEstimator_SetOptions(SUNDomEigEstimator DEE, const char *Did, const char *file_name, int argc, char *argv[])

Sets SUNDomEigEstimator options from an array of strings or a file.

Parameters:

• DEE – the SUNDomEigEstimator object.

• Did – the prefix for options to read. The default is “sundomeigestimator”.

• file_name – the name of a file containing options to read. If this is NULL or an empty string, "", then no file is read.

• argc – number of command-line arguments passed to executable.

• argv – an array of strings containing the options to set and their values.

Returns:

SUNErrCode indicating success or failure.

Note

The argc and argv arguments are typically those supplied to the user’s main routine however, this is not required. The inputs are left unchanged by SUNDomEigEstimator_SetOptions().

If the Did argument is NULL then the default prefix, sundomeigestimator, must be used for all SUNDomEigEstimator options. Whether Did is supplied or not, a "." must be used to separate an option key from the prefix. For example, when using the default Did, the option sundomeigestimator.max_iters followed by the value can be used to set the maximum number of iterations.

SUNDomEigEstimator options set via SUNDomEigEstimator_SetOptions() will overwrite any previously-set values. Options are set in the order they are given in argv and, if an option with the same prefix appears multiple times in argv, the value of the last occurrence will used.

The supported option names are noted within the documentation for the corresponding SUNDomEigEstimator functions. For options that take a sunbooleantype as input, use 1 to indicate true and 0 for false.

Warning

This function is not available in the Fortran interface.

File-based options are not yet supported, so the file_name argument should be set to either NULL or the empty string "".

Added in version 7.5.0.

SUNErrCode SUNDomEigEstimator_SetATimes(SUNDomEigEstimator DEE, void *A_data, SUNATimesFn ATimes)

This function provides a SUNATimesFn function for performing matrix-vector products, as well as a void* pointer to a data structure used by this routine, to the dominant eigenvalue estimator. This function is required when using the matrix-vector product function provided by a SUNDIALS integrator, otherwise the function is optional.

Arguments:

• DEE – a SUNDomEigEstimator object.

• A_data – pointer to structure for ATimes.

• ATimes – function pointer to perform \(Av\) product.

Return value:

A SUNErrCode.

SUNErrCode SUNDomEigEstimator_SetNumPreprocessIters(SUNDomEigEstimator DEE, int num_iters)

This optional routine sets the number of preprocessing matrix-vector multiplications, performed at the beginning of each SUNDomEigEstimator_Estimate() evaluation.

Applying preprocessing iterations may be useful if the initial guess used in SUNDomEigEstimator_Estimate() is not a good approximation of the dominant eigenvector and can help reduce some computational overhead.

Arguments:

• DEE – a SUNDomEigEstimator object.

• num_iters – the number of preprocessing iterations. Supplying a value \(< 0\), will reset the value to the implementation default.

Return value:

A SUNErrCode.

Note

When the estimator is used in a time-dependent context, different numbers of preprocessing iterations may be desired for the initial estimate than on subsequent estimations. Thus, when the estimator is used with LSRKStep (see LSRKStepSetDomEigEstimator()), the initial value of num_iters should be set with LSRKStepSetNumDomEigEstInitPreprocessIters() while the number of preprocessing iterations for subsequent estimates should be set with LSRKStepSetNumDomEigEstPreprocessIters().

This routine will be called by SUNDomEigEstimator_SetOptions() when using the key “Did.num_preprocess_iters”.

SUNErrCode SUNDomEigEstimator_SetRelTol(SUNDomEigEstimator DEE, sunrealtype rel_tol)

This optional routine sets the estimator’s relative tolerance.

Arguments:

• DEE – a SUNDomEigEstimator object.

• rel_tol – the requested eigenvalue accuracy.

Return value:

A SUNErrCode.

Note

This routine will be called by SUNDomEigEstimator_SetOptions() when using the key “Did.rel_tol”.

SUNErrCode SUNDomEigEstimator_SetMaxIters(SUNDomEigEstimator DEE, long int max_iters)

This optional routine sets the maximum number of iterations.

Arguments:

• DEE – a SUNDomEigEstimator object.

• max_iters – the maximum number of iterations.

Return value:

A SUNErrCode.

Note

This routine will be called by SUNDomEigEstimator_SetOptions() when using the key “Did.max_iters”.

SUNErrCode SUNDomEigEstimator_SetInitialGuess(SUNDomEigEstimator DEE, N_Vector q)

This optional routine sets the initial vector guess to start with.

The vector q does not need to be normalized before this set routine.

Arguments:

• DEE – a SUNDomEigEstimator object.

• q – the initial guess vector.

Return value:

A SUNErrCode.

12.2.3. SUNDomEigEstimator “get” functions

The following functions allow SUNDIALS packages to retrieve results from a dominant eigenvalue estimator. All routines are optional.

SUNErrCode SUNDomEigEstimator_GetRes(SUNDomEigEstimator DEE, sunrealtype *cur_res)

This optional routine should return the final residual from the most-recent call to SUNDomEigEstimator_Estimate().

Arguments:

• DEE – a SUNDomEigEstimator object.

• cur_res – the residual.

Return value:

A SUNErrCode.

Usage:

sunrealtype cur_res;
retval = SUNDomEigEstimator_GetRes(DEE, &cur_res);

SUNErrCode SUNDomEigEstimator_GetNumIters(SUNDomEigEstimator DEE, long int *num_iters)

This optional routine should return the number of estimator iterations performed in the most-recent call to SUNDomEigEstimator_Estimate().

Arguments:

• DEE – a SUNDomEigEstimator object.

• num_iters – the number of iterations.

Return value:

A SUNErrCode.

Usage:

long int num_iters;
retval = SUNDomEigEstimator_GetNumIters(DEE, &num_iters);

SUNErrCode SUNDomEigEstimator_GetNumATimesCalls(SUNDomEigEstimator DEE, long int *num_ATimes)

This optional routine should return the number of calls to the SUNATimesFn function.

Arguments:

• DEE – a SUNDomEigEstimator object.

• num_ATimes – the number of calls to the Atimes function.

Return value:

A SUNErrCode.

Usage:

long int num_ATimes;
retval = SUNDomEigEstimator_GetNumATimesCalls(DEE, &num_ATimes);

SUNErrCode SUNDomEigEstimator_Write(SUNDomEigEstimator DEE, FILE *outfile)

This optional routine prints the dominant eigenvalue estimator settings to the file pointer.

Arguments:

• DEE – a SUNDomEigEstimator object.

• outfile – the output stream.

Return value:

A SUNErrCode.

12.2.4. Functions provided by SUNDIALS packages

To interface with SUNDomEigEst modules, the SUNDIALS packages supply a 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.

12.2.5. The generic SUNDomEigEstimator module

SUNDIALS packages interact with dominant eigenvalue estimator implementations through the SUNDomEigEstimator class. A SUNDomEigEstimator is a pointer to the SUNDomEigEstimator_ structure:

typedef struct SUNDomEigEstimator_ *SUNDomEigEstimator

struct SUNDomEigEstimator_

The structure defining the SUNDIALS dominant eigenvalue estimator class.

void *content

Pointer to the dominant eigenvalue estimator-specific member data

SUNDomEigEstimator_Ops ops

A virtual table of dominant eigenvalue estimator operations provided by a specific implementation

SUNContext sunctx

The SUNDIALS simulation context

The virtual table structure is defined as

typedef struct SUNDomEigEstimator_Ops_ *SUNDomEigEstimator_Ops

struct SUNDomEigEstimator_Ops_

The structure defining SUNDomEigEstimator operations.

SUNErrCode (*setatimes)(SUNDomEigEstimator, void*, SUNATimesFn)

The function implementing SUNDomEigEstimator_SetATimes()

SUNErrCode (*setmaxiters)(SUNDomEigEstimator, int)

The function implementing SUNDomEigEstimator_SetMaxIters()

SUNErrCode (*setnumpreprocessiters)(SUNDomEigEstimator, int)

The function implementing SUNDomEigEstimator_SetNumPreprocessIters()

SUNErrCode (*setreltol)(SUNDomEigEstimator, sunrealtype)

The function implementing SUNDomEigEstimator_SetRelTol()

SUNErrCode (*setinitialguess)(SUNDomEigEstimator, N_Vector)

The function implementing SUNDomEigEstimator_SetInitialGuess()

SUNErrCode (*initialize)(SUNDomEigEstimator)

The function implementing SUNDomEigEstimator_Initialize()

SUNErrCode (*estimate)(SUNDomEigEstimator, sunrealtype*, sunrealtype*)

The function implementing SUNDomEigEstimator_Estimate()

sunrealtype (*getres)(SUNDomEigEstimator)

The function implementing SUNDomEigEstimator_GetRes()

int (*getnumiters)(SUNDomEigEstimator)

The function implementing SUNDomEigEstimator_GetNumIters()

long int (*getnumatimescalls)(SUNDomEigEstimator)

The function implementing SUNDomEigEstimator_GetNumATimesCalls()

SUNErrCode (*write)(SUNDomEigEstimator, FILE*)

The function implementing SUNDomEigEstimator_Write()

SUNErrCode (*destroy)(SUNDomEigEstimator*)

The function implementing SUNDomEigEstimator_Destroy()

The generic SUNDomEigEst class defines and implements the dominant eigenvalue estimator operations defined in §12.2.1 – §12.2.3. These routines are in fact only wrappers to the dominant eigenvalue estimator operations defined by a particular SUNDomEigEst implementation, which are accessed through the ops field of the SUNDomEigEstimator structure. To illustrate this point we show below the implementation of a typical dominant eigenvalue estimator operation from the SUNDomEigEstimator base class, namely SUNDomEigEstimator_Initialize(), that initializes a SUNDomEigEstimator object for use after it has been created and configured, and returns a flag denoting a successful or failed operation:

SUNErrCode SUNDomEigEstimator_Initialize(SUNDomEigEstimator DEE)
{
  return (DEE->ops->initialize(DEE));
}

Additionally, a SUNDomEigEstimator implementation may do the following:

• Define and implement additional user-callable “set” routines acting on the SUNDomEigEstimator, e.g., for setting various configuration options to tune the dominant eigenvalue estimator for a particular problem.

• Provide additional user-callable “get” routines acting on the SUNDomEigEstimator object, e.g., for returning various estimator statistics.

12.3. SUNDIALS modules SUNDomEigEstimator interface

In Table 12.1, we list the SUNDomEigEst module functions used within SUNDIALS packages. We emphasize that the user does not need to know detailed usage of dominant eigenvalue estimator functions by a SUNDIALS package in order to use it. The information is presented as an implementation detail for the interested reader.

Table 12.1 List of SUNDomEigEst functions called by a SUNDIALS module dominant eigenvalue estimator interface. Functions marked with “X” are required; functions marked with “O” are only called if they are non-NULL and functions marked with “N/A” are not applicable in the SUNDomEigEstimator implementation that is being used.

Routine	Power Iteration	Arnoldi Iteration
SUNDomEigEstimator_SetATimes()	X	X
SUNDomEigEstimator_SetMaxIters()1	O	N/A
SUNDomEigEstimator_SetNumPreprocessIters()	O	O
SUNDomEigEstimator_SetRelTol()1	O	N/A
SUNDomEigEstimator_SetInitialGuess()	O	O
SUNDomEigEstimator_Initialize()	X	X
SUNDomEigEstimator_Estimate()	X	X
SUNDomEigEstimator_GetRes()2	O	O
SUNDomEigEstimator_GetNumIters()	O	O
SUNDomEigEstimator_GetNumATimesCalls()	O	O
SUNDomEigEstimator_Write()	O	O
SUNDomEigEstimator_Destroy()3

Notes:

1. SUNDomEigEstimator_SetMaxIters() and SUNDomEigEstimator_SetRelTol() might or might not be required depending on SUNDomEigEstimator implementation that is being used. These operations should be left as NULL if it is not applicable for an estimator.

2. Although SUNDomEigEstimator_GetRes() is optional, if it is not implemented by the SUNDomEigEstimator then the interface will consider all estimates a being exact.

3. Although the interface does not call SUNDomEigEstimator_Destroy() directly, this routine should be available for users to call when cleaning up from a simulation.

12.4. The SUNDomEigEstimator_Power Module

Added in version 7.5.0.

The SUNDomEigEstimator_Power implementation of the SUNDomEigEstimator class performs the Power Iteration (PI) method [156]; this is an iterative dominant eigenvalue estimator that is designed to be compatible with any N_Vector implementation that supports a minimal subset of operations (N_VClone(), N_VDotProd(), N_VScale(), and N_VDestroy()).

Power iteration is useful for large, sparse matrices whose dominant eigenvalue is real-valued and has algebraic multiplicity one. The algorithm starts with a non-zero vector \(\mathbf{v}_{0}\). It then iteratively updates this via

\[\mathbf{v}_{k+1} = \frac{A \mathbf{v}_k}{\|A \mathbf{v}_k\|},\]

where \(\| \cdot \|\) denotes the Euclidean norm. Over successive iterations, \(\mathbf{v}_k\) converges to the eigenvector corresponding to the dominant eigenvalue of \(A\). At each step, the corresponding eigenvalue can be approximated using the Rayleigh quotient

\[\lambda_k = \frac{\mathbf{v}_k^T A \mathbf{v}_k}{\|\mathbf{v}_k\|^2}.\]

The iteration continues until the two successive eigenvalue approximations are relatively close enough to one another. That is, for some relative tolerance.

Power iteration works for the matrices that have a real dominant eigenvalue. If it is strictly greater than all others (in magnitude), convergence is guaranteed. The speed of convergence depends on the ratios of the magnitude of the first two dominant eigenvalues.

The matrix \(A\) is not required explicitly; only a routine that provides the matrix-vector product \(Av\) is required. Also, PI requires a fixed amount of memory regardless of the number of iterations.

12.4.1. SUNDomEigEstimator_Power Usage

To use SUNDomEigEstimator_Arnoldi include the header file sundomeigest/sundomeigest_power.h, and link to the library libsundials_sundomeigestpower.

The module SUNDomEigEstimator_Power provides the following user-callable routines:

SUNDomEigEstimator SUNDomEigEstimator_Power(N_Vector q, long int max_iters, sunrealtype rel_tol, SUNContext sunctx)

This constructor function creates and allocates memory for the Power iteration implementation of a SUNDomEigEstimator.

Consistency checks are performed to ensure the input vector is non-zero and supplies the necessary operations.

Parameters:

• q – the initial guess for the dominant eigenvector; this should not be a non-dominant eigenvector of the Jacobian.

• max_iters – maximum number of iterations (default 100). Supplying a value \(\leq 0\) will result in using the default value. Although this default number is not high for large matrices, it is reasonable since (1) most solvers do not need too tight tolerances and consider a safety factor, and (2) an early (less costly) termination will be a good indicator whether the power iteration is compatible.

• rel_tol – relative tolerance for convergence checks (default 0.005). A value \(\leq 0\) will result in the default value. The default has been found to small enough for many internal applications.

• sunctx – the SUNContext object.

Returns:

If successful, a SUNDomEigEstimator otherwise NULL.

Note

When used in a time-dependent context, the initial guess supplied to the constructor, q, is used only for the first SUNDomEigEstimator_Estimate() call and is overwritten with the result of the next to last Power iteration from the most recent SUNDomEigEstimator_Estimate() call. This new value is used as the initial guess for subsequent estimates.

The initial guess can be reset with SUNDomEigEstimator_SetInitialGuess().

12.4.2. SUNDomEigEstimator_Power Description

The SUNDomEigEstimator_Power module defines the content field of a SUNDomEigEstimator to be the following structure:

struct SUNDomEigEstimatorContent_Power_ {
  SUNATimesFn ATimes;
  void* ATdata;
  N_Vector* V;
  N_Vector q;
  int num_warmups;
  long int max_iters;
  long int num_iters;
  long int num_ATimes;
  sunrealtype rel_tol;
  sunrealtype res;
};

These entries of the content field contain the following information:

• ATimes - function pointer to perform the product \(Av\),

• ATData - pointer to structure for ATimes,

• V, q - N_Vector used for workspace by the PI algorithm.

• num_warmups - number of preprocessing iterations (default is 100),

• max_iters - maximum number of iterations (default is 100),

• num_iters - number of iterations (preprocessing and estimation) in the last SUNDomEigEstimator_Estimate() call,

• num_ATimes - number of calls to the ATimes function,

• rel_tol - relative tolerance for the convergence criteria (default is 0.005),

• res - the residual from the last SUNDomEigEstimator_Estimate() call.

This estimator is constructed to perform the following operations:

• During construction all N_Vector estimator data is allocated, with vectors cloned from a template N_Vector that is input, and default generic estimator parameters are set.

• User-facing “set” routines may be called to modify default estimator parameters.

• SUNDIALS packages will call SUNDomEigEstimator_SetATimes() to supply the ATimes function pointer and the related data ATData.

• In SUNDomEigEstimator_Initialize(), the estimator parameters are checked for validity and the initial eigenvector is normalized.

• In SUNDomEigEstimator_Estimate(), the initial nonzero vector \(q_0\) is preprocessed with some fixed number of Power iterations,

  \[q_1 = \frac{Aq_0}{||Aq_0||} \quad \cdots \quad q_k = \frac{Aq_{k-1}}{||Aq_{k-1}||},\]

  (see LSRKStepSetNumDomEigEstInitPreprocessIters() and LSRKStepSetNumDomEigEstPreprocessIters() for setting the number of preprocessing iterations) before computing the estimate.

The SUNDomEigEstimator_Power module defines implementations of all dominant eigenvalue estimator operations listed in §12.2:

• SUNDomEigEstimator_SetATimes_Power

• SUNDomEigEstimator_SetMaxIters_Power

• SUNDomEigEstimator_SetNumPreprocessIters_Power

• SUNDomEigEstimator_SetRelTol_Power

• SUNDomEigEstimator_Initialize_Power

• SUNDomEigEstimator_Estimate_Power

• SUNDomEigEstimator_SetInitialGuess_Power

• SUNDomEigEstimator_GetRes_Power

• SUNDomEigEstimator_GetNumIters_Power

• SUNDomEigEstimator_GetNumATimesCalls_Power

• SUNDomEigEstimator_Write_Power

• SUNDomEigEstimator_Destroy_Power

12.5. The SUNDomEigEstimator_Arnoldi Module

Added in version 7.5.0.

The SUNDomEigEstimator_Arnoldi implementation of the SUNDomEigEstimator class performs the Arnoldi Iteration method [13]; this is an iterative dominant eigenvalue estimator that is designed to be compatible with any N_Vector implementation that supports a minimal subset of operations (N_VClone(), N_VDotProd(), N_VScale(), and N_VDestroy()).

Arnoldi iteration is particularly effective for large, sparse matrices where only the dominant eigenvalue is needed. It constructs an orthonormal basis of the Krylov subspace

\[\mathcal{K}_m(A, \mathbf{v}) = \text{span}\{\mathbf{v}, A \mathbf{v}, A^2 \mathbf{v}, \dots, A^{m-1} \mathbf{v}\}\]

using the Gram-Schmidt process. The matrix \(A\) is projected onto this subspace to form a small upper Hessenberg matrix \(H_m\). The eigenvalues of \(H_m\) approximate some of the eigenvalues of \(A\); the dominant eigenvalue of \(A\) is well-approximated by the dominant eigenvalue of \(H_m\).

Arnoldi iteration works for matrices with both real and complex eigenvalues. It supports estimations with a user-specified fixed Krylov subspace dimension (at least 3). While the choice of dimension results in a prefixed amount of memory, it strictly determines the quality of the estimate. To improve the estimation accuracy, we have found that preprocessing with a number of Power iterations is particularly useful. This operation is free from any additional memory requirement and is further explained below.

The matrix \(A\) is not required explicitly; only a routine that provides an approximation of the matrix-vector product, \(Av\), is required.

12.5.1. SUNDomEigEstimator_Arnoldi Usage

To use SUNDomEigEstimator_Arnoldi include the header file sundomeigest/sundomeigest_arnoldi.h, and link to the library libsundials_sundomeigestarnoldi.

The module SUNDomEigEstimator_Arnoldi provides the following user-callable routines:

SUNDomEigEstimator SUNDomEigEstimator_Arnoldi(N_Vector q, int kry_dim, SUNContext sunctx);

This constructor function creates and allocates memory for the Arnoldi iteration implementation of a SUNDomEigEstimator.

Consistency checks are performed to ensure the input vector is non-zero and supplies the necessary operations.

Parameters:

• q – the initial guess for the dominant eigenvector; this should not be a non-dominant eigenvector of the Jacobian.

• kry_dim – the dimension of the Krylov subspace (default 3). A value \(\leq 2\) will result in using default value. This default is chosen to minimize the memory footprint.

• sunctx – the SUNContext object.

Returns:

If successful, a SUNDomEigEstimator otherwise NULL.

Note

When used in a time-dependent context, the initial guess supplied to the constructor, q, is used only in the first SUNDomEigEstimator_Estimate() call and is overwritten with the result of the most recent preprocessing iterations (see SUNDomEigEstimator_SetNumPreprocessIters()). As an initial guess too close to the dominant eigenvector may cause a breakdown in the Gram–Schmidt process within the Arnoldi iteration, users should account for this when setting the number of initial and subsequent preprocessing iterations (e.g., with LSRKStep see LSRKStepSetNumDomEigEstInitPreprocessIters() and LSRKStepSetNumDomEigEstPreprocessIters()).

The initial guess can be reset with SUNDomEigEstimator_SetInitialGuess().

12.5.2. SUNDomEigEstimator_Arnoldi Description

The SUNDomEigEstimator_Arnoldi module defines the content field of a SUNDomEigEstimator to be the following structure:

struct SUNDomEigEstimatorContent_Arnoldi_ {
  SUNATimesFn ATimes;
  void* ATdata;
  N_Vector* V;
  N_Vector q;
  int kry_dim;
  int num_warmups;
  long int num_iters;
  long int num_ATimes;
  sunrealtype* LAPACK_A;
  sunrealtype* LAPACK_wr;
  sunrealtype* LAPACK_wi;
  sunrealtype* LAPACK_work;
  snuindextype LAPACK_lwork;
  sunrealtype** LAPACK_arr;
  sunrealtype** Hes;
};

These entries of the content field contain the following information:

• ATimes - function pointer to perform the product \(Av\),

• ATData - pointer to structure for ATimes,

• V, q - vectors used for workspace by the Arnoldi algorithm.

• kry_dim - dimension of Krylov subspaces (default is 3),

• num_warmups - number of preprocessing iterations (default is 100),

• num_iters - number of iterations (preprocessing and estimation) in the last SUNDomEigEstimator_Estimate() call,

• num_ATimes - number of calls to the ATimes function,

• LAPACK_A, LAPACK_wr, LAPACK_wi, LAPACK_work - sunrealtype used for workspace by LAPACK,

• LAPACK_lwork - the size of the LAPACK_work requested by LAPACK,

• LAPACK_arr - storage for the estimated dominant eigenvalues,

• Hes - Hessenberg matrix,

This estimator is constructed to perform the following operations:

• During construction all N_Vector estimator data is allocated, with vectors cloned from a template N_Vector that is input, and default generic estimator parameters are set.

• User-facing “set” routines may be called to modify default estimator parameters.

• SUNDIALS packages will call SUNDomEigEstimator_SetATimes() to supply the ATimes function pointer and the related data ATData.

• In SUNDomEigEstimator_Initialize(), the estimator parameters are checked for validity and the remaining Arnoldi estimator memory such as LAPACK workspace is allocated.

• In SUNDomEigEstimator_Estimate(), the initial nonzero vector \(q_0\) is preprocessed with some fixed number of Power iterations,

  \[q_1 = \frac{Aq_0}{||Aq_0||} \quad \cdots \quad q_k = \frac{Aq_{k-1}}{||Aq_{k-1}||},\]

  (see LSRKStepSetNumDomEigEstInitPreprocessIters() and LSRKStepSetNumDomEigEstPreprocessIters() for setting the number of preprocessing iterations). Then, the Arnoldi iteration is performed to compute the estimate.

The SUNDomEigEstimator_Arnoldi module defines implementations of all dominant eigenvalue estimator operations listed in §12.2:

• SUNDomEigEstimator_SetATimes_Arnoldi

• SUNDomEigEstimator_SetNumPreprocessIters_Arnoldi

• SUNDomEigEstimator_Initialize_Arnoldi

• SUNDomEigEstimator_Estimate_Arnoldi

• SUNDomEigEstimator_GetNumIters_Arnoldi

• SUNDomEigEstimator_GetNumATimesCalls_Arnoldi

• SUNDomEigEstimator_Write_Arnoldi

• SUNDomEigEstimator_Destroy_Arnoldi