f12jr: FL CL CPP AD

NAG CL Interface
f12jrc (feast_complex_herm_solve)

Note: this function 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 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.

Keyword Search:

NAG Library Manual, Mark 27.3

Interfaces: FL CL CPP AD

NAG CL Interface Introduction

F12 (Sparseig) Chapter Contents

F12 (Sparseig) Chapter Introduction

f12jr: 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

Settings help

CL Name Style:
Short (a00aac)
Long (impl_details)
Full (nag_info_impl_details)

1 Purpose

f12jrc is an iterative solver used to find some of the eigenvalues and the corresponding eigenvectors of a standard or generalized eigenvalue problem defined by complex Hermitian matrices. This is part of a suite of functions that also includes f12jac, f12jbc and f12jec.

2 Specification

#include <nag.h>

void	f12jrc (void handle, Integer irevcm, Complex ze, Integer n, Complex x[], Integer pdx, Complex y[], Integer pdy, Integer m0, Integer nconv, double d[], Complex z[], Integer pdz, double eps, Integer iter, double resid[], NagError fail)

The function may be called by the names: f12jrc or nag_sparseig_feast_complex_herm_solve.

3 Description

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

λ

, and the corresponding eigenvectors,

x

, of a standard eigenvalue problem

A x = λ x

, or of a generalized eigenvalue problem

A x = λ B x

, where the coefficient matrices

A

and

B

are sparse Hermitian, and

B

is positive definite. The suite can also be used to find selected eigenvalues/eigenvectors of smaller scale, dense, Hermitian problems.

f12jrc is a reverse communication function, based on the FEAST eigensolver, described in Polizzi (2009), which finds eigenvalues using contour integration. Prior to calling f12jrc, the contour definition function f12jec is used to specify a search interval on the real line,

[E_{\min}, E_{\max}]

, within which eigenvalues will be sought (note that the eigenvalues of complex Hermitian eigenproblems are themselves real). f12jec uses this interval to define nodes and weights for an elliptical contour to be used by f12jrc.

The setup function f12jac and the contour definition function f12jec must be called before f12jrc. Between the calls to f12jac and f12jec, options may be set by calls to the option setting function f12jbc.

f12jrc 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 $ze B - A$ , where $ze$ is a point on the search contour;
–optionally, compute a factorization of the matrix ${(ze B - A)}^{H}$ (this need only be done if the factorization $ze B - A$ does not allow linear systems involving ${(ze B - A)}^{H}$ to be solved);
–solve a linear system involving $ze B - A$ or ${(ze B - A)}^{H}$ , using the factorizations above;
–compute the matrix product $x = A z$ ;
–compute the matrix product $x = B z$ ;
–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 f12jec.

4 References

Polizzi E (2009) Density-Matrix-Based Algorithms for Solving Eigenvalue Problems Phys. Rev. B. 79 115112

5 Arguments

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

1: $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: $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 $ze B - A$ suitable for solving a linear system, for example using f11dnc 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 $(ze B - A) w = y$ , overwriting $y$ with the result $w$ . The matrix $ze B - A$ 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 ${(ze B - A)}^{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 ${(ze B - A)}^{H}$ .
Note: the solve can be performed in single precision.
$irevcm = 4$: The calling program must compute the solution to the linear system ${(ze B - A)}^{H} w = y$ , overwriting $y$ with the result $w$ . If it is not possible to use the factorization of $ze B - A$ (computed when $irevcm = 1$ was returned) then the factorization of ${(ze B - A)}^{H}$ (computed when $irevcm = 2$ was returned) should be used here.
$irevcm = 5$: The calling program must compute $A z$ , storing the result in $x$ .
$irevcm = 6$: The calling program must compute $B z$ , storing the result in $x$ . If a standard eigenproblem is being solved (so that $B = I$ ) then the calling program should set $x = z$ .

On final exit:

irevcm = 0

: f12jrc 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,

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

3: $ze$ – Complex * 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 $ze B - A$ .
If $irevcm = 3$ , then, optionally, this can be used to form a factorization of ${(ze B - A)}^{H}$ .

4: $n$ – Integer Input

On entry: the order of the matrix

A

(and the order of the matrix

B

for the generalized problem) that defines the eigenvalue problem.

Constraint:

n \geq 1

5: $x [\dim]$ – Complex Input/Output

Note: the dimension, dim, of the array x must be at least

pdx \times m0

The

(i, j)

th element of the matrix

X

is stored in

x [(j - 1) \times pdx + i - 1]

On initial entry: need not be set.

On intermediate exit:

if $irevcm = 5$ , the calling program must compute $A z$ , storing the result in $x$ prior to re-entry.
If $irevcm = 6$ , the calling program must compute $B 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.

6: $pdx$ – Integer Input

On entry: the stride separating matrix row elements in the array x.

Constraint:

pdx \geq n

7: $y [\dim]$ – Complex Input/Output

Note: the dimension, dim, of the array y must be at least

pdy \times m0

The

(i, j)

th element of the matrix

Y

is stored in

y [(j - 1) \times pdy + i - 1]

On initial entry: need not be set.

On intermediate exit:

if $irevcm = 2$ , the calling program must compute the solution to the linear system $(ze B - A) 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 ${(ze B - A)}^{H} w = y$ , overwriting $y$ with the result $w$ , prior to re-entry. The linear system has m0 right-hand sides.

8: $pdy$ – Integer Input

On entry: the stride separating matrix row elements in the array y.

Constraint:

pdy \geq n

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

Constraint:

0 < m0 \leq n

10: $nconv$ – Integer * Output

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

Note: if the optional parameter

Execution Mode = Estimate

was set in the option setting function f12jbc, then nconv contains a stochastic estimate of the number of eigenvalues within the search contour.

11: $d [\dim]$ – double Input/Output

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

m0

On initial entry: if the option

Subspace = 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

Subspace = 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.

12: $z [\dim]$ – Complex Input/Output

Note: the dimension, dim, of the array z must be at least

pdz \times m0

The

(i, j)

th element of the matrix

Z

is stored in

z [(j - 1) \times pdz + i - 1]

On initial entry: if the option

Subspace = 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

Execution Mode = Subspace

was set using the option setting function f12jbc, then on final exit columns

1 : m0

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

13: $pdz$ – Integer Input

On entry: the stride separating matrix row elements in the array z.

Constraint:

pdz \geq n

14: $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

| {trace}_{k} - {trace}_{k - 1} | / \max (| E_{\min} |, | E_{\max} |)

, where

{trace}_{k}

is the sum of the eigenvalues found at the

k

th iteration, and

E_{\min}

and

E_{\max}

are the lower and upper bounds respectively of the eigenvalue search interval.

15: $iter$ – Integer * Input/Output

On initial entry: need not be set.

On exit: the number of subspace iterations performed.

16: $resid [\dim]$ – double Input/Output

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

m0

On initial entry: need not be set.

On final exit: for

i = 1, \dots, nconv

resid [i - 1]

contains the relative residual, in the

1

-norm, of the

i

th eigenpair found, that is

resid [i - 1] = ‖ A z_{i} - λ_{i} B z_{i} ‖ / (‖ B z_{i} ‖ \times \max (| E_{\min} |, | E_{\max} |))

, where

E_{\min}

and

E_{\max}

are the lower and upper bounds respectively of the eigenvalue search interval.

17: $fail$ – NagError * Input/Output

The NAG error argument (see Section 7 in the Introduction to the NAG Library CL Interface).

6 Error 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 $⟨ value ⟩$ 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, $n = ⟨ value ⟩$ .
Constraint: $n \geq 1$ .

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

On intermediate entry, $irevcm = ⟨ value ⟩$ .
Constraint: $irevcm = 1$ , $2$ , $3$ , $4$ , $5$ or $6$ .
NE_INT_2: On entry, $m0 = ⟨ value ⟩$ and $n = ⟨ value ⟩$ .
Constraint: $0 < m0 \leq n$ .

On entry, $pdx = ⟨ value ⟩$ and $n = ⟨ value ⟩$ .
Constraint: $pdx \geq n$ .

On entry, $pdy = ⟨ value ⟩$ and $n = ⟨ value ⟩$ .
Constraint: $pdy \geq n$ .

On entry, $pdz = ⟨ value ⟩$ and $n = ⟨ value ⟩$ .
Constraint: $pdz \geq n$ .
NE_INTERNAL_EIGVAL_FAIL: An internal error occurred in the reduced eigenvalue solver. A possible cause is that the matrix $B$ is not positive definite.
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 $Execution Mode = Estimate$ was set using f12jbc so only a stochastic estimate of the number of eigenvalues within the contour has been returned.

The optional parameter $Execution Mode = Subspace$ was set using f12jbc. Columns $1 : 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 $Subspace = Yes$ was set using the option setting function f12jbc but no nonzero elements were found in the supplied subspace.

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

f12jrc is not threaded in any implementation.

9 Further Comments

Ideally, when using f12jrc you should have an idea of the distribution of the eigenvalue spectrum to allow good choices of the search interval and m0 to be made. For best performance, m0 should exceed the number of eigenvalues in the search interval by a factor of approximately

1.5

A stochastic estimate of the number of eigenvalues within the search contour can be obtained by setting

Execution Mode = Estimate

in the option setting function f12jbc. In this case, f12jrc can be called with a small value of m0 (for example

~ 10

). On final output nconv will contain an estimate of the number of eigenvalues, which can then be used to guide the choice of m0.

The complex allocatable memory required by f12jrc is approximately

4 \times m0 \times m0

f12jjc can be used to solve real symmetric eigenvalue problems.

9.1 Additional Licensor

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

10 Example

This example solves the generalized eigenvalue problem

A x = λ B x

, where

A = (\begin{matrix} 0.5 + 0.0 i & 0.0 + 0.0 i & 0.0 - 0.4 i & 0.0 + 0.0 i & 0.2 - 1.1 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & 1.7 + 0.0 i & 0.0 + 0.0 i & 0.1 - 0.1 i & 0.7 - 0.5 i & 0.0 + 0.0 i \\ 0.0 + 0.4 i & 0.0 + 0.0 i & 0.2 + 0.0 i & - 0.3 + 0.5 i & 0.0 + 0.0 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & 0.0 + 0.0 i & - 0.3 - 0.5 i & - 5.8 + 0.0 i & 0.0 + 0.0 i & - 0.1 + 2.1 i \\ 0.0 + 0.0 i & 0.7 + 0.5 i & 0.0 + 0.0 i & 0.0 + 0.0 i & - 0.1 + 0.0 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i & - 0.1 + 2.1 i & 0.0 + 0.0 i & 1.6 + 0.0 i \end{matrix}),

and

B = (\begin{matrix} 8.3 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 - 0.3 i & 0.0 + 0.0 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & 6.8 + 0.0 i & 0.0 - 0.1 i & 0.0 + 0.0 i & 0.0 + 0.0 i & - 0.1 - 0.1 i \\ 0.0 + 0.0 i & 0.0 + 0.1 i & 4.7 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i \\ 0.0 + 0.3 i & 0.0 + 0.0 i & 0.0 + 0.0 i & 9.1 + 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.0 i & 8.9 + 0.0 i & 0.0 + 0.0 i \\ 0.0 + 0.0 i & - 0.1 + 0.1 i & 0.0 + 0.0 i & 0.0 + 0.0 i & 0.0 + 0.0 i & 7.6 + 0.0 i \end{matrix}) .

Only those eigenvalues lying within the interval

[- 1, 0]

are returned (together with the associated right eigenvectors). The matrices

A

and

B

are stored in symmetric coordinate storage format, with appropriate sparse matrix multiplication routines and sparse linear system solvers used.

f12jr: FL CL CPP AD

NAG CL Interfacef12jrc (feast_​complex_​herm_​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 CL Interface
f12jrc (feast_complex_herm_solve)