NAG FL Interface
f02fjf (real_​symm_​sparse_​eigsys)

1 Purpose

f02fjf finds eigenvalues and eigenvectors of a real sparse symmetric or generalized symmetric eigenvalue problem.

2 Specification

Fortran Interface
Subroutine f02fjf ( n, m, k, noits, tol, dot, image, monit, novecs, x, ldx, d, work, lwork, ruser, lruser, iuser, liuser, ifail)
Integer, Intent (In) :: n, k, novecs, ldx, lwork, lruser, liuser
Integer, Intent (Inout) :: m, noits, iuser(liuser), ifail
Real (Kind=nag_wp), External :: dot
Real (Kind=nag_wp), Intent (In) :: tol
Real (Kind=nag_wp), Intent (Inout) :: x(ldx,k), ruser(lruser)
Real (Kind=nag_wp), Intent (Out) :: d(k), work(lwork)
External :: image, monit
C Header Interface
#include <nag.h>
void  f02fjf_ (const Integer *n, Integer *m, const Integer *k, Integer *noits, const double *tol,
double (NAG_CALL *dot)(Integer *iflag, const Integer *n, const double z[], const double w[], double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser),
void (NAG_CALL *image)(Integer *iflag, const Integer *n, const double z[], double w[], double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser),
void (NAG_CALL *monit)(const Integer *istate, const Integer *nextit, const Integer *nevals, const Integer *nevecs, const Integer *k, const double f[], const double d[]),
const Integer *novecs, double x[], const Integer *ldx, double d[], double work[], const Integer *lwork, double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser, Integer *ifail)
The routine may be called by the names f02fjf or nagf_eigen_real_symm_sparse_eigsys.

3 Description

f02fjf finds the m eigenvalues of largest absolute value and the corresponding eigenvectors for the real eigenvalue problem
Cx=λ x (1)
where C is an n by n matrix such that
BC=CTB (2)
for a given positive definite matrix B. C is said to be B-symmetric. Different specifications of C allow for the solution of a variety of eigenvalue problems. For example, when
C=A  and  B=I  where  A=AT  
the routine finds the m eigenvalues of largest absolute magnitude for the standard symmetric eigenvalue problem
Ax=λx. (3)
The routine is intended for the case where A is sparse.
As a second example, when
C=B-1A  
where
A=AT  
the routine finds the m eigenvalues of largest absolute magnitude for the generalized symmetric eigenvalue problem
Ax=λBx. (4)
The routine is intended for the case where A and B are sparse.
The routine does not require C explicitly, but C is specified via image which, given an n-element vector z, computes the image w given by
w=Cz.  
For instance, in the above example, where C=B-1A, image will need to solve the positive definite system of equations Bw=Az for w.
To find the m eigenvalues of smallest absolute magnitude of (3) we can choose C=A-1 and hence find the reciprocals of the required eigenvalues, so that image will need to solve Aw=z for w, and correspondingly for (4) we can choose C=A-1B and solve Aw=Bz for w.
A table of examples of choice of image is given in Table 1. It should be remembered that the routine also returns the corresponding eigenvectors and that B is positive definite. Throughout A is assumed to be symmetric and, where necessary, nonsingularity is also assumed.
Eigenvalues
Required
Problem
  Ax=λx B=I Ax=λBx ABx=λx
Largest Compute w=Az Solve Bw=Az Compute w=ABz
Smallest (Find 1/λ) Solve Aw=z Solve Aw=Bz Solve Av=z, Bw=v
Furthest from σ
(Find λ-σ)
Compute
w=A-σIz
Solve Bw=A-σBz Compute
w=AB-σIz
Closest to σ
(Find 1/λ-σ)
Solve A-σIw=z Solve A-σBw=Bz Solve AB-σIw=z
Table 1
The Requirement of image for Various Problems.
The matrix B also need not be supplied explicitly, but is specified via dot which, given n-element vectors z and w, computes the generalized dot product wTBz.
f02fjf is based upon routine SIMITZ (see Nikolai (1979)), which is itself a derivative of the Algol procedure ritzit (see Rutishauser (1970)), and uses the method of simultaneous (subspace) iteration. (See Parlett (1998) for a description, analysis and advice on the use of the method.)
The routine performs simultaneous iteration on k>m vectors. Initial estimates to pk eigenvectors, corresponding to the p eigenvalues of C of largest absolute value, may be supplied to f02fjf. When possible k should be chosen so that the kth eigenvalue is not too close to the m required eigenvalues, but if k is initially chosen too small then f02fjf may be re-entered, supplying approximations to the k eigenvectors found so far and with k then increased.
At each major iteration f02fjf solves an r by r (rk) eigenvalue sub-problem in order to obtain an approximation to the eigenvalues for which convergence has not yet occurred. This approximation is refined by Chebyshev acceleration.

4 References

Nikolai P J (1979) Algorithm 538: Eigenvectors and eigenvalues of real generalized symmetric matrices by simultaneous iteration ACM Trans. Math. Software 5 118–125
Parlett B N (1998) The Symmetric Eigenvalue Problem SIAM, Philadelphia
Rutishauser H (1969) Computational aspects of F L Bauer's simultaneous iteration method Numer. Math. 13 4–13
Rutishauser H (1970) Simultaneous iteration method for symmetric matrices Numer. Math. 16 205–223

5 Arguments

1: n Integer Input
On entry: n, the order of the matrix C.
Constraint: n1.
2: m Integer Input/Output
On entry: m, the number of eigenvalues required.
Constraint: m1.
On exit: m, the number of eigenvalues actually found. It is equal to m if ifail=0 on exit, and is less than m if ifail=2, 3 or 4. See Sections 6 and 9 for further information.
3: k Integer Input
On entry: the number of simultaneous iteration vectors to be used. Too small a value of k may inhibit convergence, while a larger value of k incurs additional storage and additional work per iteration.
Suggested value: k=m+4 will often be a reasonable choice in the absence of better information.
Constraint: m<kn.
4: noits Integer Input/Output
On entry: the maximum number of major iterations (eigenvalue sub-problems) to be performed. If noits0, the value 100 is used in place of noits.
On exit: the number of iterations actually performed.
5: tol Real (Kind=nag_wp) Input
On entry: a relative tolerance to be used in accepting eigenvalues and eigenvectors. If the eigenvalues are required to about t significant figures, tol should be set to about 10-t. di is accepted as an eigenvalue as soon as two successive approximations to di differ by less than d~i×tol/10, where d~i is the latest approximation to di. Once an eigenvalue has been accepted, an eigenvector is accepted as soon as difi/di-dk<tol, where fi is the normalized residual of the current approximation to the eigenvector (see Section 9 for further information). The values of the fi and di can be printed from monit. If tol is supplied outside the range (ε,1.0), where ε is the machine precision, the value ε is used in place of tol.
6: dot real (Kind=nag_wp) Function, supplied by the user. External Procedure
dot must return the value wTBz for given vectors w and z. For the standard eigenvalue problem, where B=I, dot must return the dot product wTz.
The specification of dot is:
Fortran Interface
Function dot ( iflag, n, z, w, ruser, lruser, iuser, liuser)
Real (Kind=nag_wp) :: dot
Integer, Intent (In) :: n, lruser, liuser
Integer, Intent (Inout) :: iflag, iuser(liuser)
Real (Kind=nag_wp), Intent (In) :: z(n), w(n)
Real (Kind=nag_wp), Intent (Inout) :: ruser(lruser)
C Header Interface
double  dot_ (Integer *iflag, const Integer *n, const double z[], const double w[], double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser)
1: iflag Integer Input/Output
On entry: is always non-negative.
On exit: may be used as a flag to indicate a failure in the computation of wTBz. If iflag is negative on exit from dot, f02fjf will exit immediately with ifail set to iflag. Note that in this case dot must still be assigned a value.
2: n Integer Input
On entry: the number of elements in the vectors z and w and the order of the matrix B.
3: zn Real (Kind=nag_wp) array Input
On entry: the vector z for which wTBz is required.
4: wn Real (Kind=nag_wp) array Input
On entry: the vector w for which wTBz is required.
5: ruserlruser Real (Kind=nag_wp) array User Workspace
dot is called with the argument ruser as supplied to f02fjf. You should use the array ruser to supply information to dot.
6: lruser Integer Input
On entry: the dimension of the array ruser as declared in the (sub)program from which f02fjf is called.
7: iuserliuser Integer array User Workspace
dot is called with the argument iuser as supplied to f02fjf. You should use the array iuser to supply information to dot.
8: liuser Integer Input
On entry: the dimension of the array iuser as declared in the (sub)program from which f02fjf is called.
dot must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which f02fjf is called. Arguments denoted as Input must not be changed by this procedure.
Note: dot should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by f02fjf. If your code inadvertently does return any NaNs or infinities, f02fjf is likely to produce unexpected results.
7: image Subroutine, supplied by the user. External Procedure
image must return the vector w=Cz for a given vector z.
The specification of image is:
Fortran Interface
Subroutine image ( iflag, n, z, w, ruser, lruser, iuser, liuser)
Integer, Intent (In) :: n, lruser, liuser
Integer, Intent (Inout) :: iflag, iuser(liuser)
Real (Kind=nag_wp), Intent (In) :: z(n)
Real (Kind=nag_wp), Intent (Inout) :: ruser(lruser)
Real (Kind=nag_wp), Intent (Out) :: w(n)
C Header Interface
void  image_ (Integer *iflag, const Integer *n, const double z[], double w[], double ruser[], const Integer *lruser, Integer iuser[], const Integer *liuser)
1: iflag Integer Input/Output
On entry: is always non-negative.
On exit: may be used as a flag to indicate a failure in the computation of w. If iflag is negative on exit from image, f02fjf will exit immediately with ifail set to iflag.
2: n Integer Input
On entry: n, the number of elements in the vectors w and z, and the order of the matrix C.
3: zn Real (Kind=nag_wp) array Input
On entry: the vector z for which Cz is required.
4: wn Real (Kind=nag_wp) array Output
On exit: the vector w=Cz.
5: ruserlruser Real (Kind=nag_wp) array User Workspace
image is called with the argument ruser as supplied to f02fjf. You should use the array ruser to supply information to image.
6: lruser Integer Input
On entry: the dimension of the array ruser as declared in the (sub)program from which f02fjf is called.
7: iuserliuser Integer array User Workspace
image is called with the argument iuser as supplied to f02fjf. You should use the array iuser to supply information to image.
8: liuser Integer Input
On entry: the dimension of the array iuser as declared in the (sub)program from which f02fjf is called.
image must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which f02fjf is called. Arguments denoted as Input must not be changed by this procedure.
Note: image should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by f02fjf. If your code inadvertently does return any NaNs or infinities, f02fjf is likely to produce unexpected results.
8: monit Subroutine, supplied by the NAG Library or the user. External Procedure
monit is used to monitor the progress of f02fjf. monit may be the dummy subroutine f02fjz if no monitoring is actually required. (f02fjz is included in the NAG Library.) monit is called after the solution of each eigenvalue sub-problem and also just prior to return from f02fjf. The arguments istate and nextit allow selective printing by monit.
The specification of monit is:
Fortran Interface
Subroutine monit ( istate, nextit, nevals, nevecs, k, f, d)
Integer, Intent (In) :: istate, nextit, nevals, nevecs, k
Real (Kind=nag_wp), Intent (In) :: f(k), d(k)
C Header Interface
void  monit_ (const Integer *istate, const Integer *nextit, const Integer *nevals, const Integer *nevecs, const Integer *k, const double f[], const double d[])
1: istate Integer Input
On entry: specifies the state of f02fjf.
istate=0
No eigenvalue or eigenvector has just been accepted.
istate=1
One or more eigenvalues have been accepted since the last call to monit.
istate=2
One or more eigenvectors have been accepted since the last call to monit.
istate=3
One or more eigenvalues and eigenvectors have been accepted since the last call to monit.
istate=4
Return from f02fjf is about to occur.
2: nextit Integer Input
On entry: the number of the next iteration.
3: nevals Integer Input
On entry: the number of eigenvalues accepted so far.
4: nevecs Integer Input
On entry: the number of eigenvectors accepted so far.
5: k Integer Input
On entry: k, the number of simultaneous iteration vectors.
6: fk Real (Kind=nag_wp) array Input
On entry: a vector of error quantities measuring the state of convergence of the simultaneous iteration vectors. See tol and Section 9 for further details. Each element of f is initially set to the value 4.0 and an element remains at 4.0 until the corresponding vector is tested.
7: dk Real (Kind=nag_wp) array Input
On entry: di contains the latest approximation to the absolute value of the ith eigenvalue of C.
monit must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which f02fjf is called. Arguments denoted as Input must not be changed by this procedure.
9: novecs Integer Input
On entry: the number of approximate vectors that are being supplied in x. If novecs is outside the range 0,k, the value 0 is used in place of novecs.
10: xldxk Real (Kind=nag_wp) array Input/Output
On entry: if 0<novecsk, the first novecs columns of x must contain approximations to the eigenvectors corresponding to the novecs eigenvalues of largest absolute value of C. Supplying approximate eigenvectors can be useful when reasonable approximations are known, or when f02fjf is being restarted with a larger value of k. Otherwise it is not necessary to supply approximate vectors, as simultaneous iteration vectors will be generated randomly by f02fjf.
On exit: if ifail=0, 2, 3 or 4, the first m columns contain the eigenvectors corresponding to the eigenvalues returned in the first m elements of d; and the next k-m-1 columns contain approximations to the eigenvectors corresponding to the approximate eigenvalues returned in the next k-m-1 elements of d. Here m is the value returned in m, the number of eigenvalues actually found. The kth column is used as workspace.
11: ldx Integer Input
On entry: the first dimension of the array x as declared in the (sub)program from which f02fjf is called.
Constraint: ldxn.
12: dk Real (Kind=nag_wp) array Output
On exit: if ifail=0, 2, 3 or 4, the first m elements contain the first m eigenvalues in decreasing order of magnitude; and the next k-m-1 elements contain approximations to the next k-m-1 eigenvalues. Here m is the value returned in m, the number of eigenvalues actually found. dk contains the value e where -e,e is the latest interval over which Chebyshev acceleration is performed.
13: worklwork Real (Kind=nag_wp) array Workspace
14: lwork Integer Input
On entry: the dimension of the array work as declared in the (sub)program from which f02fjf is called.
Constraint: lwork3×k+maxk×k,2×n.
15: ruserlruser Real (Kind=nag_wp) array User Workspace
ruser is not used by f02fjf, but is passed directly to dot and image and may be used to pass information to these routines.
16: lruser Integer Input
On entry: the dimension of the array ruser as declared in the (sub)program from which f02fjf is called.
17: iuserliuser Integer array User Workspace
iuser is not used by f02fjf, but is passed directly to dot and image and may be used to pass information to these routines.
18: liuser Integer Input
On entry: the dimension of the array iuser as declared in the (sub)program from which f02fjf is called.
19: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 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 or 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 ifail0 on exit. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On 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 or -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
On entry, k=value and n=value.
Constraint: kn.
On entry, ldx=value and n=value.
Constraint: ldxn.
On entry, lwork is too small. Minimum size required: value.
On entry, m=value.
Constraint: m1.
On entry, m=value and k=value.
Constraint: m<k.
On entry, n=value.
Constraint: n1.
ifail=2
Not all requested eigenvalues and vectors have been obtained.
Approximations to the rth eigenvalue are oscillating rapidly indicating that severe cancellation is occurring in the rth eigenvector and so m is returned as r-1. A restart with a larger value of k may permit convergence.
ifail=3
Not all requested eigenvalues and vectors have been obtained.
The rate of convergence of the remaining eigenvectors suggests that more than noits iterations would be required and so the input value of m has been reduced. A restart with a larger value of k may permit convergence.
ifail=4
Not all requested eigenvalues and vectors have been obtained.
noits iterations have been performed. A restart, possibly with a larger value of k, may permit convergence.
ifail=5
Convergence of eigenvalue sub-problem occurred.
Not all requested eigenvalues and vectors have been obtained.
ifail<0
User set iflag negative in dot or image.
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

Eigenvalues and eigenvectors will normally be computed to the accuracy requested by the argument tol, but eigenvectors corresponding to small or to close eigenvalues may not always be computed to the accuracy requested by the argument tol. Use of the monit to monitor acceptance of eigenvalues and eigenvectors is recommended.

8 Parallelism and Performance

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

The time taken by f02fjf will be principally determined by the time taken to solve the eigenvalue sub-problem and the time taken by dot and image. The time taken to solve an eigenvalue sub-problem is approximately proportional to nk2. It is important to be aware that several calls to dot and image may occur on each major iteration.
As can be seen from Table 1, many applications of f02fjf will require the image to solve a system of linear equations. For example, to find the smallest eigenvalues of Ax=λ Bx, image needs to solve equations of the form Aw=Bz for w and routines from Chapters F01 and F04 will frequently be useful in this context. In particular, if A is a positive definite variable band matrix, f04mcf may be used after A has been factorized by f01mcf. Thus factorization need be performed only once prior to calling f02fjf. An illustration of this type of use is given in the example program.
An approximation d~h, to the ith eigenvalue, is accepted as soon as d~h and the previous approximation differ by less than d~h×tol/ 10. Eigenvectors are accepted in groups corresponding to clusters of eigenvalues that are equal, or nearly equal, in absolute value and that have already been accepted. If dr is the last eigenvalue in such a group and we define the residual rj as
rj=Cxj-yr  
where yr is the projection of Cxj, with respect to B, onto the space spanned by x1,x2,, xr, and xj is the current approximation to the jth eigenvector, then the value fi returned in monit is given by
fi = maxrjB / CxjB xB2 = xTBx  
and each vector in the group is accepted as an eigenvector if
dr fr/dr-e<tol,  
where e is the current approximation to d~k. The values of the fi are systematically increased if the convergence criteria appear to be too strict. See Rutishauser (1970) for further details.
The algorithm implemented by f02fjf differs slightly from SIMITZ (see Nikolai (1979)) in that the eigenvalue sub-problem is solved using the singular value decomposition of the upper triangular matrix R of the Gram–Schmidt factorization of Cxr, rather than forming RTR.

10 Example

This example finds the four eigenvalues of smallest absolute value and corresponding eigenvectors for the generalized symmetric eigenvalue problem Ax=λBx, where A and B are the 16 by 16 matrices
A=-14 -4 1 1 1 -4 1 1 1 -4 1 1 1 -4 1 1 1 1 -4 1 1 1 1 -4 1 1 1 1 -4 1 1 1 1 -4 1 1 1 1 -4 1 1 1 1 -4 1 1 1 1 -4 1 1 1 1 -4 1 1 1 1 -4 1 1 1 -4 1 1 1 -4 1 1 1 -4  
B=-12 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2 1 1 -2  
tol is taken as 0.0001 and 6 iteration vectors are used. f11jaf is used to factorize the matrix A, prior to calling f02fjf, and f11jcf is used within image to solve the equations Aw=Bz for w.
Output from monit occurs each time istate is nonzero. Note that the required eigenvalues are the reciprocals of the eigenvalues returned by f02fjf.

10.1 Program Text

Program Text (f02fjfe.f90)

10.2 Program Data

Program Data (f02fjfe.d)

10.3 Program Results

Program Results (f02fjfe.r)