NAG Library Routine Document

1Purpose

f02jcf solves the quadratic eigenvalue problem
 $λ2 ⁢ A + λ ⁢ B + C ⁢ x = 0 ,$
where $A$, $B$ and $C$ are real $n$ by $n$ matrices.
The routine returns the $2n$ eigenvalues, ${\lambda }_{\mathit{j}}$, for $\mathit{j}=1,2,\dots ,2n$, and can optionally return the corresponding right eigenvectors, ${x}_{j}$ and/or left eigenvectors, ${y}_{j}$ as well as estimates of the condition numbers of the computed eigenvalues and backward errors of the computed right and left eigenvectors. A left eigenvector satisfies the equation
 $yH ⁢ λ2⁢A+ λ⁢B+ C = 0 ,$
where ${y}^{\mathrm{H}}$ is the complex conjugate transpose of $y$.
$\lambda$ is represented as the pair $\left(\alpha ,\beta \right)$, such that $\lambda =\alpha /\beta$. Note that the computation of $\alpha /\beta$ may overflow and indeed $\beta$ may be zero.

2Specification

Fortran Interface
 Subroutine f02jcf ( scal, tol, n, a, lda, b, ldb, c, ldc, beta, vl, ldvl, vr, ldvr, s, bevl, bevr,
 Integer, Intent (In) :: scal, sense, n, lda, ldb, ldc, ldvl, ldvr Integer, Intent (Inout) :: ifail Integer, Intent (Out) :: iwarn Real (Kind=nag_wp), Intent (In) :: tol Real (Kind=nag_wp), Intent (Inout) :: a(lda,*), b(ldb,*), c(ldc,*), vl(ldvl,*), vr(ldvr,*), s(*), bevl(*), bevr(*) Real (Kind=nag_wp), Intent (Out) :: alphar(2*n), alphai(2*n), beta(2*n) Character (1), Intent (In) :: jobvl, jobvr
#include <nagmk26.h>
 void f02jcf_ (const Integer *scal, const char *jobvl, const char *jobvr, const Integer *sense, const double *tol, const Integer *n, double a[], const Integer *lda, double b[], const Integer *ldb, double c[], const Integer *ldc, double alphar[], double alphai[], double beta[], double vl[], const Integer *ldvl, double vr[], const Integer *ldvr, double s[], double bevl[], double bevr[], Integer *iwarn, Integer *ifail, const Charlen length_jobvl, const Charlen length_jobvr)

3Description

The quadratic eigenvalue problem is solved by linearizing the problem and solving the resulting $2n$ by $2n$ generalized eigenvalue problem. The linearization is chosen to have favourable conditioning and backward stability properties. An initial preprocessing step is performed that reveals and deflates the zero and infinite eigenvalues contributed by singular leading and trailing matrices.
The algorithm is backward stable for problems that are not too heavily damped, that is $‖B‖\le 10\sqrt{‖A‖·‖C‖}$.
Further details on the algorithm are given in Hammarling et al. (2013).

4References

Fan H -Y, Lin W.-W and Van Dooren P. (2004) Normwise scaling of second order polynomial matrices. SIAM J. Matrix Anal. Appl. 26, 1 252–256
Gaubert S and Sharify M (2009) Tropical scaling of polynomial matrices Lecture Notes in Control and Information Sciences Series 389 291–303 Springer–Verlag
Hammarling S, Munro C J and Tisseur F (2013) An algorithm for the complete solution of quadratic eigenvalue problems. ACM Trans. Math. Software. 39(3):18:1–18:119 http://eprints.maths.manchester.ac.uk/2061/

5Arguments

1:     $\mathbf{scal}$ – IntegerInput
On entry: determines the form of scaling to be performed on $A$, $B$ and $C$.
${\mathbf{scal}}=0$
No scaling.
${\mathbf{scal}}=1$ (the recommended value)
Fan, Lin and Van Dooren scaling if $\frac{‖B‖}{\sqrt{‖A‖×‖C‖}}<10$ and no scaling otherwise where $‖Z‖$ is the Frobenius norm of $Z$.
${\mathbf{scal}}=2$
Fan, Lin and Van Dooren scaling.
${\mathbf{scal}}=3$
Tropical scaling with largest root.
${\mathbf{scal}}=4$
Tropical scaling with smallest root.
Constraint: ${\mathbf{scal}}=0$, $1$, $2$, $3$ or $4$.
2:     $\mathbf{jobvl}$ – Character(1)Input
On entry: if ${\mathbf{jobvl}}=\text{'N'}$, do not compute left eigenvectors.
If ${\mathbf{jobvl}}=\text{'V'}$, compute the left eigenvectors.
If ${\mathbf{sense}}=1$, $2$, $4$, $5$, $6$ or $7$, jobvl must be set to 'V'.
Constraint: ${\mathbf{jobvl}}=\text{'N'}$ or $\text{'V'}$.
3:     $\mathbf{jobvr}$ – Character(1)Input
On entry: if ${\mathbf{jobvr}}=\text{'N'}$, do not compute right eigenvectors.
If ${\mathbf{jobvr}}=\text{'V'}$, compute the right eigenvectors.
If ${\mathbf{sense}}=1$, $3$, $4$, $5$, $6$ or $7$, jobvr must be set to 'V'.
Constraint: ${\mathbf{jobvr}}=\text{'N'}$ or $\text{'V'}$.
4:     $\mathbf{sense}$ – IntegerInput
On entry: determines whether, or not, condition numbers and backward errors are computed.
${\mathbf{sense}}=0$
Do not compute condition numbers, or backward errors.
${\mathbf{sense}}=1$
Just compute condition numbers for the eigenvalues.
${\mathbf{sense}}=2$
Just compute backward errors for the left eigenpairs.
${\mathbf{sense}}=3$
Just compute backward errors for the right eigenpairs.
${\mathbf{sense}}=4$
Compute backward errors for the left and right eigenpairs.
${\mathbf{sense}}=5$
Compute condition numbers for the eigenvalues and backward errors for the left eigenpairs.
${\mathbf{sense}}=6$
Compute condition numbers for the eigenvalues and backward errors for the right eigenpairs.
${\mathbf{sense}}=7$
Compute condition numbers for the eigenvalues and backward errors for the left and right eigenpairs.
Constraint: ${\mathbf{sense}}=0$, $1$, $2$, $3$, $4$, $5$, $6$ or $7$.
5:     $\mathbf{tol}$ – Real (Kind=nag_wp)Input
On entry: tol is used as the tolerance for making decisions on rank in the deflation procedure. If tol is zero on entry then  is used in place of tol, where machine precision is as returned by routine x02ajf. A diagonal element of a triangular matrix, $R$, is regarded as zero if $\left|{r}_{jj}\right|\le {\mathbf{tol}}×\mathrm{size}\left(X\right)$, or  when tol is zero, where $\mathrm{size}\left(X\right)$ is based on the size of the absolute values of the elements of the matrix $X$ containing the matrix $R$. See Hammarling et al. (2013) for the motivation. If tol is $-1.0$ on entry then no deflation is attempted. The recommended value for tol is zero.
6:     $\mathbf{n}$ – IntegerInput
On entry: the order of the matrices $A$, $B$ and $C$.
Constraint: ${\mathbf{n}}\ge 0$.
7:     $\mathbf{a}\left({\mathbf{lda}},*\right)$ – Real (Kind=nag_wp) arrayInput/Output
Note: the second dimension of the array a must be at least ${\mathbf{n}}$.
On entry: the $n$ by $n$ matrix $A$.
On exit: a is used as internal workspace, but if ${\mathbf{jobvl}}=\text{'V'}$ or ${\mathbf{jobvr}}=\text{'V'}$, a is restored on exit.
8:     $\mathbf{lda}$ – IntegerInput
On entry: the first dimension of the array a as declared in the (sub)program from which f02jcf is called.
Constraint: ${\mathbf{lda}}\ge {\mathbf{n}}$.
9:     $\mathbf{b}\left({\mathbf{ldb}},*\right)$ – Real (Kind=nag_wp) arrayInput/Output
Note: the second dimension of the array b must be at least ${\mathbf{n}}$.
On entry: the $n$ by $n$ matrix $B$.
On exit: b is used as internal workspace, but is restored on exit.
10:   $\mathbf{ldb}$ – IntegerInput
On entry: the first dimension of the array b as declared in the (sub)program from which f02jcf is called.
Constraint: ${\mathbf{ldb}}\ge {\mathbf{n}}$.
11:   $\mathbf{c}\left({\mathbf{ldc}},*\right)$ – Real (Kind=nag_wp) arrayInput/Output
Note: the second dimension of the array c must be at least ${\mathbf{n}}$.
On entry: the $n$ by $n$ matrix $C$.
On exit: c is used as internal workspace, but if ${\mathbf{jobvl}}=\text{'V'}$ or ${\mathbf{jobvr}}=\text{'V'}$, c is restored on exit.
12:   $\mathbf{ldc}$ – IntegerInput
On entry: the first dimension of the array c as declared in the (sub)program from which f02jcf is called.
Constraint: ${\mathbf{ldc}}\ge {\mathbf{n}}$.
13:   $\mathbf{alphar}\left(2×{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput
On exit: ${\mathbf{alphar}}\left(\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,2n$, contains the real part of ${\alpha }_{j}$ for the $j$th eigenvalue pair $\left({\alpha }_{j},{\beta }_{j}\right)$ of the quadratic eigenvalue problem.
14:   $\mathbf{alphai}\left(2×{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput
On exit: ${\mathbf{alphai}}\left(\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,2n$, contains the imaginary part of ${\alpha }_{j}$ for the $j$th eigenvalue pair $\left({\alpha }_{j},{\beta }_{j}\right)$ of the quadratic eigenvalue problem. If ${\mathbf{alphai}}\left(j\right)$ is zero then the $j$th eigenvalue is real; if ${\mathbf{alphai}}\left(j\right)$ is positive then the $j$th and $\left(j+1\right)$th eigenvalues are a complex conjugate pair, with ${\mathbf{alphai}}\left(j+1\right)$ negative.
15:   $\mathbf{beta}\left(2×{\mathbf{n}}\right)$ – Real (Kind=nag_wp) arrayOutput
On exit: ${\mathbf{beta}}\left(\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,2n$, contains the second part of the $j$th eigenvalue pair $\left({\alpha }_{j},{\beta }_{j}\right)$ of the quadratic eigenvalue problem, with ${\beta }_{j}\ge 0$. Infinite eigenvalues have ${\beta }_{j}$ set to zero.
16:   $\mathbf{vl}\left({\mathbf{ldvl}},*\right)$ – Real (Kind=nag_wp) arrayOutput
Note: the second dimension of the array vl must be at least $2×{\mathbf{n}}$ if ${\mathbf{jobvl}}=\text{'V'}$.
On exit: if ${\mathbf{jobvl}}=\text{'V'}$, the left eigenvectors ${y}_{j}$ are stored one after another in the columns of vl, in the same order as the corresponding eigenvalues. If the $j$th eigenvalue is real, then ${y}_{j}={\mathbf{vl}}\left(:,j\right)$, the $j$th column of vl. If the $j$th and $\left(j+1\right)$th eigenvalues form a complex conjugate pair, then ${y}_{j}={\mathbf{vl}}\left(:,j\right)+i×{\mathbf{vl}}\left(:,j+1\right)$ and ${y}_{j+1}={\mathbf{vl}}\left(:,j\right)-i×{\mathbf{vl}}\left(:,j+1\right)$. Each eigenvector will be normalized with length unity and with the element of largest modulus real and positive.
If ${\mathbf{jobvl}}=\text{'N'}$, vl is not referenced.
17:   $\mathbf{ldvl}$ – IntegerInput
On entry: the first dimension of the array vl as declared in the (sub)program from which f02jcf is called.
Constraint: ${\mathbf{ldvl}}\ge {\mathbf{n}}$.
18:   $\mathbf{vr}\left({\mathbf{ldvr}},*\right)$ – Real (Kind=nag_wp) arrayOutput
Note: the second dimension of the array vr must be at least $2×{\mathbf{n}}$ if ${\mathbf{jobvr}}=\text{'V'}$.
On exit: if ${\mathbf{jobvr}}=\text{'V'}$, the right eigenvectors ${x}_{j}$ are stored one after another in the columns of vr, in the same order as the corresponding eigenvalues. If the $j$th eigenvalue is real, then ${x}_{j}={\mathbf{vr}}\left(:,j\right)$, the $j$th column of vr. If the $j$th and $\left(j+1\right)$th eigenvalues form a complex conjugate pair, then ${x}_{j}={\mathbf{vr}}\left(:,j\right)+i×{\mathbf{vr}}\left(:,j+1\right)$ and ${x}_{j+1}={\mathbf{vr}}\left(:,j\right)-i×{\mathbf{vr}}\left(:,j+1\right)$. Each eigenvector will be normalized with length unity and with the element of largest modulus real and positive.
If ${\mathbf{jobvr}}=\text{'N'}$, vr is not referenced.
19:   $\mathbf{ldvr}$ – IntegerInput
On entry: the first dimension of the array vr as declared in the (sub)program from which f02jcf is called.
Constraint: ${\mathbf{ldvr}}\ge {\mathbf{n}}$.
20:   $\mathbf{s}\left(*\right)$ – Real (Kind=nag_wp) arrayOutput
Note: the dimension of the array s must be at least $2×{\mathbf{n}}$ if ${\mathbf{sense}}=1$, $5$, $6$ or $7$.
Also: computing the condition numbers of the eigenvalues requires that both the left and right eigenvectors be computed.
On exit: if ${\mathbf{sense}}=1$, $5$, $6$ or $7$, ${\mathbf{s}}\left(j\right)$ contains the condition number estimate for the $j$th eigenvalue (large condition numbers imply that the problem is near one with multiple eigenvalues). Infinite condition numbers are returned as the largest model real number (x02alf).
If ${\mathbf{sense}}=0$, $2$, $3$ or $4$, s is not referenced.
21:   $\mathbf{bevl}\left(*\right)$ – Real (Kind=nag_wp) arrayOutput
Note: the dimension of the array bevl must be at least $2×{\mathbf{n}}$ if ${\mathbf{sense}}=2$, $4$, $5$ or $7$.
On exit: if ${\mathbf{sense}}=2$, $4$, $5$ or $7$, ${\mathbf{bevl}}\left(j\right)$ contains the backward error estimate for the computed left eigenpair $\left({\lambda }_{j},{y}_{j}\right)$.
If ${\mathbf{sense}}=0$, $1$, $3$ or $6$, bevl is not referenced.
22:   $\mathbf{bevr}\left(*\right)$ – Real (Kind=nag_wp) arrayOutput
Note: the dimension of the array bevr must be at least $2×{\mathbf{n}}$ if ${\mathbf{sense}}=3$, $4$, $6$ or $7$.
On exit: if ${\mathbf{sense}}=3$, $4$, $6$ or $7$, ${\mathbf{bevr}}\left(j\right)$ contains the backward error estimate for the computed right eigenpair $\left({\lambda }_{j},{x}_{j}\right)$.
If ${\mathbf{sense}}=0$, $1$, $2$ or $5$, bevr is not referenced.
23:   $\mathbf{iwarn}$ – IntegerOutput
On exit: iwarn will be positive if there are warnings, otherwise iwarn will be $0$.
If ${\mathbf{ifail}}={\mathbf{0}}$ then:
• if ${\mathbf{iwarn}}=1$ then one, or both, of the matrices $A$ and $C$ is zero. In this case no scaling is performed, even if ${\mathbf{scal}}>0$;
• if ${\mathbf{iwarn}}=2$ then the matrices $A$ and $C$ are singular, or nearly singular, so the problem is potentially ill-posed;
• if ${\mathbf{iwarn}}=3$ then both the conditions for ${\mathbf{iwarn}}=1$ and ${\mathbf{iwarn}}=2$ above, apply. If ${\mathbf{iwarn}}=4$, $‖{\mathbf{b}}‖\ge 10\sqrt{‖A‖·‖C‖}$ and backward stability cannot be guaranteed.
If ${\mathbf{ifail}}={\mathbf{2}}$, iwarn returns the value of info from f08xaf (dgges).
If ${\mathbf{ifail}}={\mathbf{3}}$, iwarn returns the value of info from f08waf (dggev).
24:   $\mathbf{ifail}$ – IntegerInput/Output
On entry: ifail must be set to $0$, . If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value  is recommended. If the output of error messages is undesirable, then the value $1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is $0$. When the value  is used it is essential to test the value of ifail on exit.
On exit: ${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see Section 6).

6Error Indicators and Warnings

If on entry ${\mathbf{ifail}}=0$ or $-1$, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
${\mathbf{ifail}}=1$
The quadratic matrix polynomial is nonregular (singular).
${\mathbf{ifail}}=2$
The $QZ$ iteration failed in f08xaf (dgges).
iwarn returns the value of info returned by f08xaf (dgges). This failure is unlikely to happen, but if it does, please contact NAG.
${\mathbf{ifail}}=3$
The $QZ$ iteration failed in f08waf (dggev).
iwarn returns the value of info returned by f08waf (dggev). This failure is unlikely to happen, but if it does, please contact NAG.
${\mathbf{ifail}}=-1$
On entry, ${\mathbf{scal}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{scal}}=0$, $1$, $2$, $3$ or $4$.
${\mathbf{ifail}}=-2$
On entry, ${\mathbf{jobvl}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{jobvl}}=\text{'N'}$ or $\text{'V'}$.
On entry, ${\mathbf{sense}}=〈\mathit{\text{value}}〉$ and ${\mathbf{jobvl}}=〈\mathit{\text{value}}〉$.
Constraint: when ${\mathbf{jobvl}}=\text{'N'}$, ${\mathbf{sense}}=0$ or $3$,
when ${\mathbf{jobvl}}=\text{'V'}$, ${\mathbf{sense}}=1$, $2$, $4$, $5$, $6$ or $7$.
${\mathbf{ifail}}=-3$
On entry, ${\mathbf{jobvr}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{jobvr}}=\text{'N'}$ or $\text{'V'}$.
On entry, ${\mathbf{sense}}=〈\mathit{\text{value}}〉$ and ${\mathbf{jobvr}}=〈\mathit{\text{value}}〉$.
Constraint: when ${\mathbf{jobvr}}=\text{'N'}$, ${\mathbf{sense}}=0$ or $2$,
when ${\mathbf{jobvr}}=\text{'V'}$, ${\mathbf{sense}}=1$, $3$, $4$, $5$, $6$ or $7$.
${\mathbf{ifail}}=-4$
On entry, ${\mathbf{sense}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{sense}}=0$, $1$, $2$, $3$, $4$, $5$, $6$ or $7$.
${\mathbf{ifail}}=-6$
On entry, ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{n}}\ge 0$.
${\mathbf{ifail}}=-8$
On entry, ${\mathbf{lda}}=〈\mathit{\text{value}}〉$ and ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{lda}}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=-10$
On entry, ${\mathbf{ldb}}=〈\mathit{\text{value}}〉$ and ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{ldb}}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=-12$
On entry, ${\mathbf{ldc}}=〈\mathit{\text{value}}〉$ and ${\mathbf{n}}=〈\mathit{\text{value}}〉$.
Constraint: ${\mathbf{ldc}}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=-17$
On entry, ${\mathbf{ldvl}}=〈\mathit{\text{value}}〉$, ${\mathbf{n}}=〈\mathit{\text{value}}〉$ and ${\mathbf{jobvl}}=〈\mathit{\text{value}}〉$.
Constraint: when ${\mathbf{jobvl}}=\text{'V'}$, ${\mathbf{ldvl}}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=-19$
On entry, ${\mathbf{ldvr}}=〈\mathit{\text{value}}〉$, ${\mathbf{n}}=〈\mathit{\text{value}}〉$ and ${\mathbf{jobvr}}=〈\mathit{\text{value}}〉$.
Constraint: when ${\mathbf{jobvr}}=\text{'V'}$, ${\mathbf{ldvr}}\ge {\mathbf{n}}$.
${\mathbf{ifail}}=-99$
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-399$
Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.
${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7Accuracy

The algorithm is backward stable for problems that are not too heavily damped, that is $‖B‖\le 10\sqrt{‖A‖·‖C‖}$.

8Parallelism and Performance

f02jcf is threaded by NAG for parallel execution in multithreaded implementations of the NAG Library.
f02jcf 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.

None.

10Example

To solve the quadratic eigenvalue problem
 $λ2 A + λB + C x = 0$
where
 $A = 1 2 2 3 1 1 3 2 1 , B = 3 2 1 2 1 3 1 3 2 and C = 1 1 2 2 3 1 3 1 2 .$
The example also returns the left eigenvectors, condition numbers for the computed eigenvalues and backward errors of the computed right and left eigenpairs.

10.1Program Text

Program Text (f02jcfe.f90)

10.2Program Data

Program Data (f02jcfe.d)

10.3Program Results

Program Results (f02jcfe.r)