NAG CL Interface f12juc (feast_poly_symm_solve)
Note:this function usesoptional parametersto define choices in the problem specification. If you wish to use default settings for all of the optional parameters, then the option setting function f12jbc need not be called.
If, however, you wish to reset some or all of the settings please refer to Section 11 in f12jbc for a detailed description of the specification of the optional parameters.
f12juc is an iterative solver used to find some of the eigenvalues and the corresponding eigenvectors of a polynomial eigenvalue problem defined by complex symmetric matrices. This is part of a suite of functions that also includes f12jac,f12jbc,f12jfcandf12jgc.
The function may be called by the names: f12juc or nag_sparseig_feast_poly_symm_solve.
3Description
The suite of functions is designed to calculate some of the eigenvalues, $\lambda $, and the corresponding eigenvectors, $x$, of a polynomial eigenvalue problem ${\sum}_{i=0}^{p}{\lambda}^{i}{A}_{i}x=0$ where the coefficient matrices ${A}_{i}$ are sparse, complex and symmetric. The suite can also be used to find selected eigenvalues/eigenvectors of smaller scale, dense problems.
f12juc is a reverse communication function, based on the FEAST eigensolver, described in Polizzi (2009), which finds eigenvalues using contour integration. Prior to calling f12juc, one of the contour definition functions f12jfcorf12jgc must be used to define nodes and weights for a contour around a region in the complex plane within which eigenvalues will be sought. f12jfcorf12jgc uses this interval to define nodes and weights for the contour to be used by f12juc.
The setup function f12jac and one of the contour definition functions f12jfcorf12jgc must be called before f12juc. Between the calls to f12jac and f12jfcorf12jgc, options may be set by calls to the option setting function f12jbc.
f12juc uses reverse communication, i.e., it returns repeatedly to the calling program with the argument irevcm (see Section 5) set to specified values which require the calling program to carry out one of the following tasks:
–compute a factorization of the matrix ${\sum}_{i=0}^{p}{{\mathbf{ze}}}^{i}{A}_{i}$, where ${\mathbf{ze}}$ is a point on the search contour;
–solve a linear system involving ${\sum}_{i=0}^{p}{{\mathbf{ze}}}^{i}{A}_{i}$ using the factorizations above;
–compute the matrix product $x={A}_{k}z$, for a given value of $k$;
–notify the completion of the computation.
The number of contour points, the number of iterations, and other options can all be set using the option setting function f12jbc (see Section 11.1 in f12jbc for details on setting options and of the default settings). The search contour itself is defined by a call to f12jfcorf12jgc.
4References
Gavin B, Międlar A and Polizzi E (2018) FEAST Eigensolver for Nonlinear Eigenvalue Problems J. Comput. Phys.27 107
Polizzi E (2009) Density-Matrix-Based Algorithms for Solving Eigenvalue Problems Phys. Rev. B. 79 115112
5Arguments
Note: this function uses reverse communication. Its use involves an initial entry, intermediate exits and re-entries, and a final exit, as indicated by the argument irevcm. Between intermediate exits and re-entries, all arguments other thanx and y must remain unchanged.
1: $\mathbf{handle}$ – void *Input
On entry: the handle to the internal data structure used by the NAG FEAST suite. It needs to be initialized by f12jac. It must not be changed between calls to the NAG FEAST suite.
2: $\mathbf{irevcm}$ – Integer *Input/Output
On initial entry: ${\mathbf{irevcm}}=0$, otherwise an error condition will be raised.
On intermediate re-entry: must be unchanged from its previous exit value. Changing irevcm to any other value between calls will result in an error.
On intermediate exit:
has the following meanings.
${\mathbf{irevcm}}=1$
The calling program must compute a factorization of the matrix ${\sum}_{i=0}^{p}{{\mathbf{ze}}}^{i}{A}_{i}$ suitable for solving a linear system, for example using f11dnc which computes an incomplete $LU$ factorization of a complex sparse matrix. All arguments to the routine must remain unchanged.
Note: the factorization can be computed in single precision.
${\mathbf{irevcm}}=2$
The calling program must compute the solution to the linear system $\left({\sum}_{i=0}^{p}{{\mathbf{ze}}}^{i}{A}_{i}\right)w=y$, overwriting $y$ with the result $w$. The matrix ${\sum}_{i=0}^{p}{{\mathbf{ze}}}^{i}{A}_{i}$ has previously been factorized (when ${\mathbf{irevcm}}=1$ was returned) and this factorization can be reused here.
Note: the solve can be performed in single precision.
${\mathbf{irevcm}}=3$
The calling program must compute ${A}_{k}z$, storing the result in $x$.
On final exit: ${\mathbf{irevcm}}=0$: f12juc has completed its tasks. The value of fail determines whether the iteration has been successfully completed, or whether errors have been detected.
Constraint:
on initial entry, ${\mathbf{irevcm}}=0$; on re-entry irevcm must remain unchanged.
Note: the matrices $x$, $y$ and $z$ referred to in this section are all of size ${\mathbf{n}}\times {\mathbf{m0}}$ and are stored in the arrays x, y and z, respectively.
Note: any values you return to f12juc as part of the reverse communication procedure should not include floating-point NaN (Not a Number) or infinity values, since these are not handled by f12juc. If your code inadvertently does return any NaNs or infinities, f12juc is likely to produce unexpected results.
3: $\mathbf{deg}$ – IntegerInput
On entry: the degree, $p$, of the polynomial eigenvalue problem.
Constraint:
${\mathbf{deg}}>0$.
4: $\mathbf{ze}$ – Complex *Input/Output
On initial entry: need not be set.
On intermediate exit:
contains the current point on the contour.
If ${\mathbf{irevcm}}=1$, then this must be used by the calling program to form a factorization of the matrix ${\sum}_{i=0}^{p}{{\mathbf{ze}}}^{i}{A}_{i}$.
5: $\mathbf{n}$ – IntegerInput
On entry: the order of the matrix $A$ (and the order of the matrix $B$ for the generalized problem) that defines the eigenvalue problem.
Note: the dimension, dim, of the array x
must be at least
${\mathbf{pdx}}\times {\mathbf{m0}}$.
The $(i,j)$th element of the matrix $X$ is stored in ${\mathbf{x}}\left[(j-1)\times {\mathbf{pdx}}+i-1\right]$.
On initial entry: need not be set.
On intermediate exit:
if ${\mathbf{irevcm}}=3$, the calling program must compute ${A}_{k}z$, storing the result in $x$ prior to re-entry.
Note: the matrices $x$, $y$ and $z$ referred to in this section are all of size ${\mathbf{n}}\times {\mathbf{m0}}$ and are stored in the arrays x, y and z, respectively.
the matrices $x$ and $z$ are stored in the first m0 columns of the arrays x and z, respectively.
7: $\mathbf{pdx}$ – IntegerInput
On entry: the stride separating matrix row elements in the array x.
Note: the dimension, dim, of the array y
must be at least
${\mathbf{pdy}}\times {\mathbf{m0}}$.
The $(i,j)$th element of the matrix $Y$ is stored in ${\mathbf{y}}\left[(j-1)\times {\mathbf{pdy}}+i-1\right]$.
On initial entry: need not be set.
On intermediate exit:
if ${\mathbf{irevcm}}=2$, the calling program must compute the solution to the linear system $\left({\sum}_{i=0}^{p}{{\mathbf{ze}}}^{i}{A}_{i}\right)w=y$, overwriting $y$ with the result $w$, prior to re-entry. The linear system has m0 right-hand sides.
9: $\mathbf{pdy}$ – IntegerInput
On entry: the stride separating matrix row elements in the array y.
Constraint:
${\mathbf{pdy}}\ge {\mathbf{n}}$.
10: $\mathbf{k}$ – Integer *Output
On intermediate exit:
if ${\mathbf{irevcm}}=3$, denotes which of the coefficient matrices ${A}_{k}$ needs to be used to form ${A}_{k}z$.
11: $\mathbf{m0}$ – Integer *Input/Output
On initial entry: the size of the search subspace used to find the eigenvalues. This should exceed the number of eigenvalues within the search contour. See Section 9 for further details.
On intermediate re-entry: m0 must remain unchanged.
On exit: if the initial search subspace was found by f12juc to be too large, then a new smaller suitable choice is returned.
Constraint:
$0<{\mathbf{m0}}\le {\mathbf{n}}$.
12: $\mathbf{nconv}$ – Integer *Output
On exit: the number of eigenvalues found within the search contour.
Note: the dimension, dim, of the array d
must be at least
${\mathbf{m0}}$.
On initial entry: if the option ${\mathbf{Subspace}}=\mathrm{Yes}$ was set using the option setting function f12jbc, then d should contain an initial guess at the eigenvalues lying within the eigenvector search subspace (this subspace should be specified by z), otherwise d need not be set.
On final exit: the first nconv entries in d contain the eigenvalues.
Note: if the option ${\mathbf{Subspace}}=\mathrm{Yes}$ was set using the option setting function f12jbc, then on final exit d contains an estimate of the eigenvalues after a single contour integral.
Note: the dimension, dim, of the array z
must be at least
${\mathbf{pdz}}\times {\mathbf{m0}}$.
The $(i,j)$th element of the matrix $Z$ is stored in ${\mathbf{z}}\left[(j-1)\times {\mathbf{pdz}}+i-1\right]$.
On initial entry: if the option ${\mathbf{Subspace}}=\mathrm{Yes}$ was set using the option setting function f12jbc, then z should contain an initial guess at the eigenvector search subspace, otherwise z need not be set.
On intermediate exit:
must not be changed.
On final exit: the first nconv columns of z contain the eigenvectors corresponding to the eigenvalues found within the contour.
Note: if the option ${\mathbf{Execution\; Mode}}=\mathrm{Subspace}$ was set using the option setting function f12jbc, then on final exit columns $1:{\mathbf{m0}}$ of z contain the current search subspace after one contour integral.
15: $\mathbf{pdz}$ – IntegerInput
On entry: the stride separating matrix row elements in the array z.
Constraint:
${\mathbf{pdz}}\ge {\mathbf{n}}$.
16: $\mathbf{eps}$ – double *Input/Output
On initial entry: need not be set.
On exit: the relative error on the trace. At iteration $k$, eps is given by the expression $|{\mathrm{trace}}_{k}-{\mathrm{trace}}_{k-1}|/(\left|{E}_{\mathrm{mid}}\right|+r)$, where ${\mathrm{trace}}_{k}$ is the sum of the eigenvalues found at the $k$th iteration, ${E}_{\mathrm{mid}}$ is the centre of mass of the contour and $r$ is the radius of the circle, centred at ${E}_{\mathrm{mid}}$ which just encloses the contour.
17: $\mathbf{iter}$ – Integer *Input/Output
On initial entry: need not be set.
On exit: the number of subspace iterations performed.
Note: the dimension, dim, of the array resid
must be at least
${\mathbf{m0}}$.
On initial entry: need not be set.
On final exit: for $i=1,\dots ,{\mathbf{nconv}}$, ${\mathbf{resid}}\left[i-1\right]$ contains the relative residual, in the $1$-norm, of the $i$th eigenpair found, that is ${\mathbf{resid}}\left[i-1\right]=\Vert \left({\sum}_{j=0}^{p}{\lambda}_{i}^{j}{A}_{j}\right){z}_{i}\Vert /(\Vert {\lambda}^{p}{A}_{p}{z}_{i}\Vert \times (\left|{E}_{\mathrm{mid}}\right|+r))$, where ${E}_{\mathrm{mid}}$ is the centre of mass of the contour and $r$ is the radius of the circle, centred at ${E}_{\mathrm{mid}}$ which just encloses the contour.
19: $\mathbf{fail}$ – NagError *Input/Output
The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).
6Error Indicators and Warnings
NE_ALLOC_FAIL
Dynamic memory allocation failed.
See Section 3.1.2 in the Introduction to the NAG Library CL Interface for further information.
NE_BAD_PARAM
On entry, argument $\u27e8\mathit{\text{value}}\u27e9$ had an illegal value.
NE_EIGENVALUES
No eigenvalues were found within the search contour.
NE_HANDLE
Either the contour setting function f12jec has not been called prior to the first call of this function or the supplied handle has become corrupted.
NE_INT
On entry, ${\mathbf{deg}}=\u27e8\mathit{\text{value}}\u27e9$. Constraint: ${\mathbf{deg}}>0$.
On entry, ${\mathbf{n}}=\u27e8\mathit{\text{value}}\u27e9$. Constraint: ${\mathbf{n}}\ge 1$.
On initial entry, ${\mathbf{irevcm}}=\u27e8\mathit{\text{value}}\u27e9$. Constraint: ${\mathbf{irevcm}}=0$.
On intermediate entry, ${\mathbf{irevcm}}=\u27e8\mathit{\text{value}}\u27e9$. Constraint: ${\mathbf{irevcm}}=1$, $2$, $3$, $4$, $5$ or $6$.
NE_INT_2
On entry, ${\mathbf{m0}}=\u27e8\mathit{\text{value}}\u27e9$ and ${\mathbf{n}}=\u27e8\mathit{\text{value}}\u27e9$. Constraint: $0<{\mathbf{m0}}\le {\mathbf{n}}$.
On entry, ${\mathbf{pdx}}=\u27e8\mathit{\text{value}}\u27e9$ and ${\mathbf{n}}=\u27e8\mathit{\text{value}}\u27e9$. Constraint: ${\mathbf{pdx}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{pdy}}=\u27e8\mathit{\text{value}}\u27e9$ and ${\mathbf{n}}=\u27e8\mathit{\text{value}}\u27e9$. Constraint: ${\mathbf{pdy}}\ge {\mathbf{n}}$.
On entry, ${\mathbf{pdz}}=\u27e8\mathit{\text{value}}\u27e9$ and ${\mathbf{n}}=\u27e8\mathit{\text{value}}\u27e9$. Constraint: ${\mathbf{pdz}}\ge {\mathbf{n}}$.
NE_INTERNAL_EIGVAL_FAIL
An internal error occurred in the reduced eigenvalue solver. Please contact NAG.
NE_INTERNAL_ERROR
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
See Section 7.5 in the Introduction to the NAG Library CL Interface for further information.
NE_NO_LICENCE
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library CL Interface for further information.
NE_NO_SOLUTION
The optional parameter ${\mathbf{Execution\; Mode}}=\mathrm{Estimate}$ was set using f12jbc.
nconv contains a stochastic estimate of the number of eigenvalues
within the contour.
The optional parameter ${\mathbf{Execution\; Mode}}=\mathrm{Subspace}$ was set using f12jbc. Columns $1:{\mathbf{m0}}$ of z contain the search subspace after one contour integral and d contains an estimate of the eigenvalues.
NE_TOO_MANY_ITER
The function did not converge after the maximum number of iterations. The results returned may still be useful, however they might be improved by increasing the maximum number of iterations using the option setting routine f12jbc, increasing the size of the eigenvector search subspace, m0, or experimenting with the choice of contour. Note that the returned eigenvalues and eigenvectors, together with the returned value of m0, can be used as the initial estimates for a new iteration of the solver.
NE_TOO_SMALL
The size of the eigenvector search subspace, m0, is too small.
NE_ZERO_VALUE
The option ${\mathbf{Subspace}}=\mathrm{Yes}$ was set using the option setting function f12jbc but no nonzero elements were found in the supplied subspace.
7Accuracy
A gauge on the accuracy of the computation can be obtained by looking at eps, the relative error on the trace, and the residuals, stored in resid.
Note: the factorizations and linear system solves required when ${\mathbf{irevcm}}=1$ or $2$ can be performed in single precision, without any loss of accuracy in the final eigenvalues and eigenvectors.
8Parallelism and Performance
Background information to multithreading can be found in the Multithreading documentation.
f12juc is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f12juc makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this function. Please also consult the Users' Note for your implementation for any additional implementation-specific information.
9Further Comments
Ideally, when using f12juc you should have an idea of the distribution of the eigenvalue spectrum to allow good choices of the search contour and m0 to be made. For best performance, m0 should exceed the number of eigenvalues within the search contour by a factor of approximately $1.5$. Note that a polynomial eigenvalue problem of order $n$ and degree $p$ can have up to $n\times p$ eigenvalues. However m0 must not exceed n.
The complex allocatable memory required by f12juc is approximately $2\times {\mathbf{m0}}\times {\mathbf{m0}}\times {\mathbf{deg}}\times {\mathbf{deg}}$.
9.1Additional Licensor
Parts of the code for f12juc are distributed under the BSD software License. Please refer to Library Licensors for further details.
10Example
This example solves the symmetric quadratic eigenproblem $({\lambda}^{2}{A}_{2}+\lambda {A}_{1}+{A}_{0})x=0$, where