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) |
|
C++ Header Interface
#include <nag.h> extern "C" {
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
eigenvalues of largest absolute value and the corresponding eigenvectors for the real eigenvalue problem
where
is an
by
matrix such that
for a given positive definite matrix
.
is said to be
-symmetric. Different specifications of
allow for the solution of a variety of eigenvalue problems. For example, when
the routine finds the
eigenvalues of largest absolute magnitude for the standard symmetric eigenvalue problem
The routine is intended for the case where
is sparse.
As a second example, when
where
the routine finds the
eigenvalues of largest absolute magnitude for the generalized symmetric eigenvalue problem
The routine is intended for the case where
and
are sparse.
The routine does not require
explicitly, but
is specified via
image which, given an
-element vector
, computes the image
given by
For instance, in the above example, where
,
image will need to solve the positive definite system of equations
for
.
To find the
eigenvalues of smallest absolute magnitude of
(3) we can choose
and hence find the reciprocals of the required eigenvalues, so that
image will need to solve
for
, and correspondingly for
(4) we can choose
and solve
for
.
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
is positive definite. Throughout
is assumed to be symmetric and, where necessary, nonsingularity is also assumed.
Eigenvalues Required |
Problem |
|
|
|
|
Largest |
Compute |
Solve |
Compute |
Smallest (Find ) |
Solve |
Solve |
Solve , |
Furthest from
(Find ) |
Compute
|
Solve |
Compute
|
Closest to
(Find ) |
Solve |
Solve |
Solve |
Table 1
The Requirement of
image for Various Problems.
The matrix
also need not be supplied explicitly, but is specified via
dot which, given
-element vectors
and
, computes the generalized dot product
.
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 vectors. Initial estimates to eigenvectors, corresponding to the eigenvalues of of largest absolute value, may be supplied to f02fjf. When possible should be chosen so that the th eigenvalue is not too close to the required eigenvalues, but if is initially chosen too small then f02fjf may be re-entered, supplying approximations to the eigenvectors found so far and with then increased.
At each major iteration f02fjf solves an by () 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:
– Integer
Input
-
On entry: , the order of the matrix .
Constraint:
.
-
2:
– Integer
Input/Output
-
On entry: , the number of eigenvalues required.
Constraint:
.
On exit:
, the number of eigenvalues actually found. It is equal to
if
on exit, and is less than
if
,
or
. See
Sections 6 and
9 for further information.
-
3:
– 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:
will often be a reasonable choice in the absence of better information.
Constraint:
.
-
4:
– Integer
Input/Output
-
On entry: the maximum number of major iterations (eigenvalue sub-problems) to be performed. If
, the value
is used in place of
noits.
On exit: the number of iterations actually performed.
-
5:
– 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
significant figures,
tol should be set to about
.
is accepted as an eigenvalue as soon as two successive approximations to
differ by less than
, where
is the latest approximation to
. Once an eigenvalue has been accepted, an eigenvector is accepted as soon as
, where
is the normalized residual of the current approximation to the eigenvector (see
Section 9 for further information). The values of the
and
can be printed from
monit. If
tol is supplied outside the range (
), where
is the
machine precision, the value
is used in place of
tol.
-
6:
– real (Kind=nag_wp) Function, supplied by the user.
External Procedure
-
dot must return the value
for given vectors
and
. For the standard eigenvalue problem, where
,
dot must return the dot product
.
The specification of
dot is:
Fortran Interface
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) |
|
C++ Header Interface
#include <nag.h> extern "C" {
double |
dot_ (Integer &iflag, const Integer &n, const double z[], const double w[], double ruser[], const Integer &lruser, Integer iuser[], const Integer &liuser) |
}
|
-
1:
– 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
. 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:
– Integer
Input
-
On entry: the number of elements in the vectors and and the order of the matrix .
-
3:
– Real (Kind=nag_wp) array
Input
-
On entry: the vector for which is required.
-
4:
– Real (Kind=nag_wp) array
Input
-
On entry: the vector for which is required.
-
5:
– 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:
– Integer
Input
-
On entry: the dimension of the array
ruser as declared in the (sub)program from which
f02fjf is called.
-
7:
– 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:
– 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:
– Subroutine, supplied by the user.
External Procedure
-
image must return the vector
for a given vector
.
The specification of
image is:
Fortran Interface
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) |
|
C++ Header Interface
#include <nag.h> extern "C" {
void |
image_ (Integer &iflag, const Integer &n, const double z[], double w[], double ruser[], const Integer &lruser, Integer iuser[], const Integer &liuser) |
}
|
-
1:
– 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
. If
iflag is negative on exit from
image,
f02fjf will exit immediately with
ifail set to
iflag.
-
2:
– Integer
Input
-
On entry: , the number of elements in the vectors and , and the order of the matrix .
-
3:
– Real (Kind=nag_wp) array
Input
-
On entry: the vector for which is required.
-
4:
– Real (Kind=nag_wp) array
Output
-
On exit: the vector .
-
5:
– 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:
– Integer
Input
-
On entry: the dimension of the array
ruser as declared in the (sub)program from which
f02fjf is called.
-
7:
– 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:
– 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:
– 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
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[]) |
|
C++ Header Interface
#include <nag.h> extern "C" {
void |
monit_ (const Integer &istate, const Integer &nextit, const Integer &nevals, const Integer &nevecs, const Integer &k, const double f[], const double d[]) |
}
|
-
1:
– Integer
Input
-
On entry: specifies the state of
f02fjf.
- No eigenvalue or eigenvector has just been accepted.
- One or more eigenvalues have been accepted since the last call to monit.
- One or more eigenvectors have been accepted since the last call to monit.
- One or more eigenvalues and eigenvectors have been accepted since the last call to monit.
- Return from f02fjf is about to occur.
-
2:
– Integer
Input
-
On entry: the number of the next iteration.
-
3:
– Integer
Input
-
On entry: the number of eigenvalues accepted so far.
-
4:
– Integer
Input
-
On entry: the number of eigenvectors accepted so far.
-
5:
– Integer
Input
-
On entry: , the number of simultaneous iteration vectors.
-
6:
– 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
and an element remains at
until the corresponding vector is tested.
-
7:
– Real (Kind=nag_wp) array
Input
-
On entry: contains the latest approximation to the absolute value of the th eigenvalue of .
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:
– Integer
Input
-
On entry: the number of approximate vectors that are being supplied in
x. If
novecs is outside the range
, the value
is used in place of
novecs.
-
10:
– Real (Kind=nag_wp) array
Input/Output
-
On entry: if
, the first
novecs columns of
x must contain approximations to the eigenvectors corresponding to the
novecs eigenvalues of largest absolute value of
. 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
,
,
or
, the first
columns contain the eigenvectors corresponding to the eigenvalues returned in the first
elements of
d; and the next
columns contain approximations to the eigenvectors corresponding to the approximate eigenvalues returned in the next
elements of
d. Here
is the value returned in
m, the number of eigenvalues actually found. The
th column is used as workspace.
-
11:
– Integer
Input
-
On entry: the first dimension of the array
x as declared in the (sub)program from which
f02fjf is called.
Constraint:
.
-
12:
– Real (Kind=nag_wp) array
Output
-
On exit: if
,
,
or
, the first
elements contain the first
eigenvalues in decreasing order of magnitude; and the next
elements contain approximations to the next
eigenvalues. Here
is the value returned in
m, the number of eigenvalues actually found.
contains the value
where
is the latest interval over which Chebyshev acceleration is performed.
-
13:
– Real (Kind=nag_wp) array
Workspace
-
14:
– Integer
Input
-
On entry: the dimension of the array
work as declared in the (sub)program from which
f02fjf is called.
Constraint:
.
-
15:
– 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:
– Integer
Input
-
On entry: the dimension of the array
ruser as declared in the (sub)program from which
f02fjf is called.
-
17:
– 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:
– Integer
Input
-
On entry: the dimension of the array
iuser as declared in the (sub)program from which
f02fjf is called.
-
19:
– Integer
Input/Output
-
On entry:
ifail must be set to
,
or
to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of means that an error message is printed while a value of means that it is not.
If halting is not appropriate, the value
or
is recommended. If message printing is undesirable, then the value
is recommended. Otherwise, the value
is recommended since useful values can be provided in some output arguments even when
on exit.
When the value or is used it is essential to test the value of ifail on exit.
On exit:
unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
or
, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
-
On entry, and .
Constraint: .
On entry, and .
Constraint: .
On entry,
lwork is too small. Minimum size required:
.
On entry, .
Constraint: .
On entry, and .
Constraint: .
On entry, .
Constraint: .
-
Not all requested eigenvalues and vectors have been obtained.
Approximations to the
th eigenvalue are oscillating rapidly indicating that severe cancellation is occurring in the
th eigenvector and so
m is returned as
. A restart with a larger value of
k may permit convergence.
-
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.
-
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.
-
Convergence of eigenvalue sub-problem occurred.
Not all requested eigenvalues and vectors have been obtained.
-
User set
iflag negative in
dot or
image.
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.
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.
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.
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
. 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
,
image needs to solve equations of the form
for
and routines from
Chapters F01 and
F04
will frequently be useful in this context. In particular, if
is a positive definite variable band matrix,
f04mcf may be used after
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
, to the
th eigenvalue, is accepted as soon as
and the previous approximation differ by less than
. 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
is the last eigenvalue in such a group and we define the residual
as
where
is the projection of
, with respect to
, onto the space spanned by
, and
is the current approximation to the
th eigenvector, then the value
returned in
monit is given by
and each vector in the group is accepted as an eigenvector if
where
is the current approximation to
. The values of the
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
of the Gram–Schmidt factorization of
, rather than forming
.
10
Example
This example finds the four eigenvalues of smallest absolute value and corresponding eigenvectors for the generalized symmetric eigenvalue problem
, where
and
are the
by
matrices
tol is taken as
and
iteration vectors are used.
f11jaf is used to factorize the matrix
, prior to calling
f02fjf, and
f11jcf is used within
image to solve the equations
for
.
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
10.2
Program Data
10.3
Program Results