f12jv: FL CL CPP AD

NAG FL Interface
f12jvf (feast_poly_gen_solve)

Note: this routine uses optional parameters to define choices in the problem specification. If you wish to use default settings for all of the optional parameters, then the option setting routine f12jbf need not be called. If, however, you wish to reset some or all of the settings please refer to Section 11 in f12jbf for a detailed description of the specification of the optional parameters.

Keyword Search:

NAG Library Manual, Mark 28.3

Interfaces: FL CL CPP AD

NAG FL Interface Introduction

F12 (Sparseig) Chapter Contents

F12 (Sparseig) Chapter Introduction

f12jv: FL CL CPP AD

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

▸▿ 9 Further Comments

9.1 Additional Licensor

▸▿ 10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1 Purpose

f12jvf is an iterative solver used to find some of the eigenvalues and the corresponding eigenvectors of a polynomial eigenvalue problem defined by complex nonsymmetric matrices. This is part of a suite of routines that also includes f12jaf, f12jbf, f12jff and f12jgf.

2 Specification

Fortran Interface

Subroutine f12jvf (

handle, irevcm, deg, ze, n, x, ldx, y, ldy, k, m0, nconv, d, z, ldz, eps, iter, resid, ifail)

Integer, Intent (In)	::	deg, n, ldx, ldy, ldz
Integer, Intent (Inout)	::	irevcm, m0, iter, ifail
Integer, Intent (Out)	::	k, nconv
Real (Kind=nag_wp), Intent (Inout)	::	eps, resid(*)
Complex (Kind=nag_wp), Intent (Inout)	::	ze, x(ldx,), y(ldy,), d(), z(ldz,)
Type (c_ptr), Intent (In)	::	handle

C Header Interface

#include <nag.h>

void

f12jvf_ (void **handle, Integer *irevcm, const Integer *deg, Complex *ze, const Integer *n, Complex x[], const Integer *ldx, Complex y[], const Integer *ldy, Integer *k, Integer *m0, Integer *nconv, Complex d[], Complex z[], const Integer *ldz, double *eps, Integer *iter, double resid[], Integer *ifail)

The routine may be called by the names f12jvf or nagf_sparseig_feast_poly_gen_solve.

3 Description

The suite of routines is designed to calculate some of the eigenvalues,

λ

, and the corresponding eigenvectors,

x

, of a polynomial eigenvalue problem

\sum_{i = 0}^{p} λ^{i} A_{i} x = 0

where the coefficient matrices

A_{i}

are sparse, complex and nonsymmetric. The suite can also be used to find selected eigenvalues/eigenvectors of smaller scale, dense problems.

f12jvf is a reverse communication routine, based on the FEAST eigensolver, described in Polizzi (2009), which finds eigenvalues using contour integration. Prior to calling f12jvf, one of the contour definition routines f12jff or f12jgf must be used to define nodes and weights for a contour around a region in the complex plane within which eigenvalues will be sought. f12jff or f12jgf uses this interval to define nodes and weights for the contour to be used by f12jvf.

The setup routine f12jaf and one of the contour definition routines f12jff or f12jgf must be called before f12jvf. Between the calls to f12jaf and f12jff or f12jgf, options may be set by calls to the option setting routine f12jbf.

f12jvf 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} {ze}^{i} A_{i}$ , where $ze$ is a point on the search contour;
–optionally, compute a factorization of the matrix ${(\sum_{i = 0}^{p} {ze}^{i} A_{i})}^{H}$ (this need only be done if the factorization $\sum_{i = 0}^{p} {ze}^{i} A_{i}$ does not allow linear systems involving ${(\sum_{i = 0}^{p} {ze}^{i} A_{i})}^{H}$ to be solved);
–solve a linear system involving $\sum_{i = 0}^{p} {ze}^{i} A_{i}$ or ${(\sum_{i = 0}^{p} {ze}^{i} A_{i})}^{H}$ , using the factorizations above;
–compute the matrix product $x = A_{k} z$ , for a given value of $k$ ;
–compute the matrix product $x = {A^{H}}_{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 routine f12jbf (see Section 11.1 in f12jbf for details on setting options and of the default settings). The search contour itself is defined by a call to f12jff or f12jgf.

4 References

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

5 Arguments

Note: this routine 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 than x and y must remain unchanged.

1: $handle$ – Type (c_ptr) Input

On entry: the handle to the internal data structure used by the NAG FEAST suite. It needs to be initialized by f12jaf. It must not be changed between calls to the NAG FEAST suite.

2: $irevcm$ – Integer Input/Output

On initial entry:

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.

$irevcm = 1$: The calling program must compute a factorization of the matrix $\sum_{i = 0}^{p} {ze}^{i} A_{i}$ suitable for solving a linear system, for example using f11dnf which computes an incomplete $L U$ factorization of a complex sparse matrix. All arguments to the routine must remain unchanged.
Note: the factorization can be computed in single precision.
$irevcm = 2$: The calling program must compute the solution to the linear system $(\sum_{i = 0}^{p} {ze}^{i} A_{i}) w = y$ , overwriting $y$ with the result $w$ . The matrix $\sum_{i = 0}^{p} {ze}^{i} A_{i}$ has previously been factorized (when $irevcm = 1$ was returned) and this factorization can be reused here.
Note: the solve can be performed in single precision.
$irevcm = 3$: Optionally, the calling program must compute a factorization of the matrix ${(\sum_{i = 0}^{p} {ze}^{i} A_{i})}^{H}$ . This need only be done if it is not possible to use the factorization computed when $irevcm = 1$ was returned to solve linear systems involving ${(\sum_{i = 0}^{p} {ze}^{i} A_{i})}^{H}$ . If this factorization is to be computed, then the factorization from $irevcm = 1$ must not be overwritten.
Note: the factorization can be performed in single precision.
$irevcm = 4$: The calling program must compute the solution to the linear system ${(\sum_{i = 0}^{p} {ze}^{i} A_{i})}^{H} w = y$ , overwriting $y$ with the result $w$ . If it is not possible to use the factorization of $\sum_{i = 0}^{p} {ze}^{i} A_{i}$ (computed when $irevcm = 1$ was returned) then the factorization of ${(\sum_{i = 0}^{p} {ze}^{i} A_{i})}^{H}$ (computed when $irevcm = 2$ was returned) should be used here.
Note: the factorization can be performed in single precision.
$irevcm = 5$: The calling program must compute $A_{k} z$ , storing the result in $x$ .
$irevcm = 6$: The calling program must compute $A_{k}^{H} z$ , storing the result in $x$ .

On final exit:

irevcm = 0

: f12jvf has completed its tasks. The value of ifail determines whether the iteration has been successfully completed, or whether errors have been detected.

Constraint: on initial entry,

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

n \times m0

and are stored in the arrays x, y and z, respectively.

Note: any values you return to f12jvf 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 f12jvf. If your code does inadvertently return any NaNs or infinities, f12jvf is likely to produce unexpected results.

3: $\deg$ – Integer Input

On entry: the degree,

p

, of the polynomial eigenvalue problem.

Constraint:

\deg > 0

4: $ze$ – Complex (Kind=nag_wp) Input/Output

On initial entry: need not be set.

On intermediate exit: contains the current point on the contour.

If $irevcm = 1$ , then this must be used by the calling program to form a factorization of the matrix $\sum_{i = 0}^{p} {ze}^{i} A_{i}$ .
If $irevcm = 3$ , then, optionally, this can be used to form a factorization of ${(\sum_{i = 0}^{p} {ze}^{i} A_{i})}^{H}$ .

5: $n$ – Integer Input

On entry: the order of the matrices

A_{i}

that define the eigenvalue problem.

Constraint:

n \geq 1

6: $x (ldx, *)$ – Complex (Kind=nag_wp) array Input/Output

Note: the second dimension of the array x must be at least

2 \times m0

On initial entry: need not be set.

On intermediate exit:

if $irevcm = 5$ , the calling program must compute $A_{k} z$ , storing the result in $x$ prior to re-entry.
If $irevcm = 6$ , the calling program must compute $A_{k}^{H} z$ , storing the result in $x$ prior to re-entry.

Note: the matrices

x

and

z

are stored in the first m0 columns of the arrays x and z, respectively.

7: $ldx$ – Integer Input

On entry: the first dimension of the array x as declared in the (sub)program from which f12jvf is called.

Constraint:

ldx \geq n

8: $y (ldy, *)$ – Complex (Kind=nag_wp) array Input/Output

Note: the second dimension of the array y must be at least

m0

On initial entry: need not be set.

On intermediate exit:

if $irevcm = 2$ , the calling program must compute the solution to the linear system $(\sum_{i = 0}^{p} {ze}^{i} A_{i}) w = y$ , overwriting $y$ with the result $w$ , prior to re-entry. The linear system has m0 right-hand sides.
If $irevcm = 4$ , the calling program must compute the solution to the linear system ${(\sum_{i = 0}^{p} {ze}^{i} A_{i})}^{H} w = y$ , overwriting $y$ with the result $w$ , prior to re-entry. The linear system has m0 right-hand sides.

9: $ldy$ – Integer Input

On entry: the first dimension of the array y as declared in the (sub)program from which f12jvf is called.

Constraint:

ldy \geq n

10: $k$ – Integer Output

On intermediate exit:

if $irevcm = 5$ , denotes which of the coefficient matrices $A_{k}$ needs to be used to form $A_{k} z$ .
If $irevcm = 6$ , denotes which of the coefficient matrices $A_{k}$ needs to be used to form $A_{k}^{H} z$ .

11: $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 f12jvf to be too large, then a new smaller suitable choice is returned.

Constraint:

0 < m0 \leq n

12: $nconv$ – Integer Output

On exit: the number of eigenvalues found within the search contour.

13: $d (*)$ – Complex (Kind=nag_wp) array Input/Output

Note: the dimension of the array d must be at least

m0

On initial entry: if the option

Subspace = Yes

was set using the option setting routine f12jbf, 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

Subspace = Yes

was set using the option setting routine f12jbf, then on final exit d contains an estimate of the eigenvalues after a single contour integral.

14: $z (ldz, *)$ – Complex (Kind=nag_wp) array Input/Output

Note: the second dimension of the array z must be at least

2 \times m0

On initial entry: if the option

Subspace = Yes

was set using the option setting routine f12jbf, then z should contain an initial guess at the eigenvector search subspace, otherwise z need not be set.

Note: if the option

Subspace = Yes

was set but only right eigenvectors are to be computed, then only the first m0 columns of z need to be set. However, if the optional parameter

Eigenvectors = Two-sided

was set using f12jbf, so that both left and right eigenvectors are returned, then the first m0 columns of z should contain the initial guess for the right eigenvector search subspace and the remaining m0 columns should contain an initial guess for the left eigenvector search subspace.

On intermediate exit: must not be changed.

On final exit: the first nconv columns of z contain the right eigenvectors corresponding to the eigenvalues found within the contour. If the option

Eigenvectors = Two-sided

was set in the option setting routine f12jbf, the columns

m0 + 1 : m0 + nconv

of z contain the left eigenvectors corresponding to the eigenvalues found within the contour.

Note: if the option

Execution Mode = Subspace

was set using the option setting routine f12jbf, then on final exit columns

1 : m0

of z contain the current search subspace after one contour integral.

15: $ldz$ – Integer Input

On entry: the first dimension of the array z as declared in the (sub)program from which f12jvf is called.

Constraint:

ldz \geq n

16: $eps$ – Real (Kind=nag_wp) 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

| {trace}_{k} - {trace}_{k - 1} | / (| E_{mid} | + r)

, where

{trace}_{k}

is the sum of the eigenvalues found at the

k

th iteration,

E_{mid}

is the centre of mass of the contour and

r

is the radius of the circle, centred at

E_{mid}

which just encloses the contour.

17: $iter$ – Integer Input/Output

On initial entry: need not be set.

On exit: the number of subspace iterations performed.

18: $resid (*)$ – Real (Kind=nag_wp) array Input/Output

Note: the dimension of the array resid must be at least

2 \times m0

On initial entry: need not be set.

On final exit: for

i = 1, \dots, nconv

resid (i)

contains the relative residual, in the

1

-norm, of the

i

th eigenpair found, that is

resid (i) = ‖ (\sum_{j = 0}^{p} λ_{i}^{j} A_{j}) z_{i} ‖ / (‖ λ^{p} A_{p} z_{i} ‖ \times (| E_{mid} | + r))

, where

E_{mid}

is the centre of mass of the contour and

r

is the radius of the circle, centred at

E_{mid}

which just encloses the contour.

If the optional parameter

Eigenvectors = Two-sided

was set in the option setting routine f12jbf, the entries

m0 + 1 : m0 + nconv

of resid contain corresponding residuals for the left eigenvectors,

resid (i) = ‖ {(\sum_{j = 0}^{p} λ_{i}^{j} A_{j})}^{H} z_{i} ‖ / (‖ λ^{p} A_{p}^{H} z_{i} ‖ \times (| E_{mid} | + r))

19: $ifail$ – Integer Input/Output

On initial entry: ifail must be set to

0

−1

1

to set behaviour on detection of an error; these values have no effect when no error is detected.

A value of

0

causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of

−1

means that an error message is printed while a value of

1

means that it is not.

If halting is not appropriate, the value

−1

1

is recommended. If message printing is undesirable, then the value

1

is recommended. Otherwise, the value

−1

is recommended since useful values can be provided in some output arguments even when

ifail \neq 0

on exit. When the value $- 1$ or $1$ is used it is essential to test the value of ifail on exit.

On final exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

ifail = 0

−1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

$ifail = 1$: Either one of the contour setting routines f12jff or f12jgf has not been called prior to the first call of this routine or the supplied handle has become corrupted.

$ifail = 2$: No eigenvalues were found within the search contour.

$ifail = 3$: The routine 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 f12jbf, 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.

$ifail = 4$: The size of the eigenvector search subspace, m0, is too small.

$ifail = 5$: The optional parameter $Execution Mode = Subspace$ was set using f12jbf. Columns $1 : m0$ of z contain the search subspace after one contour integral and d contains an estimate of the eigenvalues.

$ifail = 6$: The optional parameter $Execution Mode = Estimate$ was set using f12jbf. nconv contains a stochastic estimate of the number of eigenvalues within the contour.

$ifail = 7$: An internal error occurred in the reduced eigenvalue solver. Please contact NAG.

$ifail = 8$: On entry, $n = ⟨ value ⟩$ .
Constraint: $n \geq 1$ .

$ifail = 9$: On entry, $m0 = ⟨ value ⟩$ and $n = ⟨ value ⟩$ .
Constraint: $0 < m0 \leq n$ .

$ifail = 10$: On entry, $ldx = ⟨ value ⟩$ and $n = ⟨ value ⟩$ .
Constraint: $ldx \geq n$ .

$ifail = 11$: On entry, $ldy = ⟨ value ⟩$ and $n = ⟨ value ⟩$ .
Constraint: $ldy \geq n$ .

$ifail = 12$: On entry, $ldz = ⟨ value ⟩$ and $n = ⟨ value ⟩$ .
Constraint: $ldz \geq n$ .

$ifail = 13$: On initial entry, $irevcm = ⟨ value ⟩$ .
Constraint: $irevcm = 0$ .

On intermediate entry, $irevcm = ⟨ value ⟩$ .
Constraint: $irevcm = 1$ , $2$ , $3$ , $4$ , $5$ or $6$ .

$ifail = 14$: On entry, $\deg = ⟨ value ⟩$ .
Constraint: $\deg > 0$ .

$ifail = 15$: The option $Subspace = Yes$ was set using the option setting routine f12jbf but no nonzero elements were found in the supplied subspace.

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

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

irevcm = 1

2

3

4

can be performed in single precision, without any loss of accuracy in the final eigenvalues and eigenvectors.

8 Parallelism and Performance

f12jvf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.

f12jvf 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 routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9 Further Comments

Ideally, when using f12jvf 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 in 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 f12jvf is approximately

2 \times m0 \times m0 \times \deg \times \deg

9.1 Additional Licensor

Parts of the code for f12jvf are distributed under the BSD software License. Please refer to Library Licensors for further details.

10 Example

This example solves the quadratic eigenproblem

(λ^{2} A_{2} + λ A_{1} + A_{0}) x = 0

, where

A_{0} = (\begin{matrix} 0.0 + 0.0 i & 7.5 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i \\ - 1.7 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.4 i & 0.2 + 0.0 i \\ 0.0 + 0.0 i & - 1.6 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i \end{matrix}),

A_{1} = (\begin{matrix} 0.0 + 1.0 i & 0.0 + 0.0 i & 0.0 + 1.2 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & 0.0 + 2.9 i & 0.0 + 0.0 i & 0.0 - 3.4 i \\ 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.1 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 1.2 i \end{matrix}),

and

A_{2} = (\begin{matrix} 1.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & 1.7 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & 0.0 + 0.0 i & 0.2 + 0.0 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i & - 5.8 + 0.0 i \end{matrix}) .

The matrices are stored in coordinate storage format and appropriate sparse linear algebra routines are used throughout. The contour within which eigenvalues are sought is a semicircle connecting the points

0.5 + 1.5 i

and

- 1.5 - 0.5 i

f12jv: FL CL CPP AD

NAG FL Interfacef12jvf (feast_​poly_​gen_​solve)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

9.1 Additional Licensor

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
f12jvf (feast_poly_gen_solve)