NAG Library Routine Document
F02EKF
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, you need only read Sections 1 to 9 of this document. If, however, you wish to reset some or all of the settings this must be done by calling the option setting routine F12ADF from the user-supplied subroutine OPTION. Please refer to Section 10 for a detailed description of the specification of the optional parameters.
1 Purpose
F02EKF computes selected eigenvalues and eigenvectors of a real sparse general matrix.
2 Specification
SUBROUTINE F02EKF ( |
N, NNZ, A, ICOLZP, IROWIX, NEV, NCV, SIGMA, MONIT, OPTION, NCONV, W, V, LDV, RESID, IUSER, RUSER, IFAIL) |
INTEGER |
N, NNZ, ICOLZP(N+1), IROWIX(NNZ), NEV, NCV, NCONV, LDV, IUSER(*), IFAIL |
REAL (KIND=nag_wp) |
A(NNZ), SIGMA, V(LDV,*), RESID(NEV+1), RUSER(*) |
COMPLEX (KIND=nag_wp) |
W(NCV) |
EXTERNAL |
MONIT, OPTION |
|
3 Description
F02EKF computes selected eigenvalues and the corresponding right eigenvectors of a real sparse general matrix
:
A specified number, , of eigenvalues , or the shifted inverses , may be selected either by largest or smallest modulus, largest or smallest real part, or, largest or smallest imaginary part. Convergence is generally faster when selecting larger eigenvalues, smaller eigenvalues can always be selected by choosing a zero inverse shift (). When eigenvalues closest to a given real value are required then the shifted inverses of largest magnitude should be selected with shift equal to the required real value.
Note that even though is real, and may be complex. If is an eigenvector corresponding to a complex eigenvalue , then the complex conjugate vector is the eigenvector corresponding to the complex conjugate eigenvalue . The eigenvalues in a complex conjugate pair and are either both selected or both not selected.
The sparse matrix
is stored in compressed column storage (CCS) format. See
Section 2.1.3 in the F11 Chapter Introduction.
F02EKF uses an implicitly restarted Arnoldi iterative method to converge approximations to a set of required eigenvalues and corresponding eigenvectors. Further algorithmic information is given in
Section 8 while a fuller discussion is provided in the
F12 Chapter Introduction. If shifts are to be performed then operations using shifted inverse matrices are performed using a direct sparse solver; further information on the solver used is provided in the
F11 Chapter Introduction.
4 References
Golub G H and Van Loan C F (1996) Matrix Computations (3rd Edition) Johns Hopkins University Press, Baltimore
Lehoucq R B, Sorensen D C and Yang C (1998) ARPACK Users' Guide: Solution of Large-scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods SIAM, Philidelphia
5 Parameters
- 1: N – INTEGERInput
On entry: , the order of the matrix .
Constraint:
.
- 2: NNZ – INTEGERInput
On entry: the dimension of the array
A as declared in the (sub)program from which F02EKF is called. The number of nonzero elements of the matrix
and, if a nonzero shifted inverse is to be applied, all diagonal elements. A nonzero diagonal element is counted once in the latter case.
Constraint:
.
- 3: A(NNZ) – REAL (KIND=nag_wp) arrayInput/Output
On entry: the array of nonzero elements (and diagonal elements if a nonzero inverse shift is to be applied) of the by general matrix .
On exit: if a nonzero shifted inverse is to be applied then the diagonal elements of
have the shift value, as supplied in
SIGMA, subtracted.
- 4: ICOLZP() – INTEGER arrayInput
On entry:
contains the index in
A of the start of column
, for
;
must contain the value
. Thus the number of nonzero elements in column
of
is
; when shifts are applied this includes diagonal elements irrespective of value. See
Section 2.1.3 in the F11 Chapter Introduction.
- 5: IROWIX() – INTEGER arrayInput
On entry:
contains the row index for each entry in
A. See
Section 2.1.3 in the F11 Chapter Introduction.
- 6: NEV – INTEGERInput
On entry: the dimension of the array
RESID as declared in the (sub)program from which F02EKF is called. The number of eigenvalues to be computed.
Constraint:
.
- 7: NCV – INTEGERInput
On entry: the dimension of the array
W as declared in the (sub)program from which F02EKF is called.
The number of Arnoldi basis vectors to use during the computation.
At present there is no
a priori analysis to guide the selection of
NCV relative to
NEV. However, it is recommended that
. If many problems of the same type are to be solved, you should experiment with increasing
NCV while keeping
NEV fixed for a given test problem. This will usually decrease the required number of matrix-vector operations but it also increases the work and storage required to maintain the orthogonal basis vectors. The optimal ‘cross-over’ with respect to CPU time is problem dependent and must be determined empirically.
Constraint:
.
- 8: SIGMA – REAL (KIND=nag_wp)Input
On entry: if the
Shifted Inverse Real mode has been selected then
SIGMA contains the real shift used; otherwise
SIGMA is not referenced. This mode can be selected by setting the appropriate options in the user-supplied subroutine
OPTION.
- 9: MONIT – SUBROUTINE, supplied by the NAG Library or the user.External Procedure
MONIT is used to monitor the progress of F02EKF.
MONIT may be the dummy subroutine F02EKZ if no monitoring is actually required. (F02EKZ is included in the NAG Library.)
MONIT is called after the solution of each eigenvalue sub-problem and also just prior to return from F02EKF.
The specification of
MONIT is:
INTEGER |
NCV, NITER, NCONV, ISTAT, IUSER(*) |
REAL (KIND=nag_wp) |
RZEST(NCV), RUSER(*) |
COMPLEX (KIND=nag_wp) |
W(NCV) |
|
- 1: NCV – INTEGERInput
On entry: the dimension of the arrays
W and
RZEST. The number of Arnoldi basis vectors used during the computation.
- 2: NITER – INTEGERInput
On entry: the number of the current Arnoldi iteration.
- 3: NCONV – INTEGERInput
On entry: the number of converged eigenvalues so far.
- 4: W(NCV) – COMPLEX (KIND=nag_wp) arrayInput
-
On entry: the first
NCONV elements of
W contain the converged approximate eigenvalues.
- 5: RZEST(NCV) – REAL (KIND=nag_wp) arrayInput
-
On entry: the first
NCONV elements of
RZEST contain the Ritz estimates (error bounds) on the converged approximate eigenvalues.
- 6: ISTAT – INTEGERInput/Output
On entry: set to zero.
On exit: if set to a nonzero value F02EKF returns immediately with .
- 7: IUSER() – INTEGER arrayUser Workspace
- 8: RUSER() – REAL (KIND=nag_wp) arrayUser Workspace
-
MONIT is called with the parameters
IUSER and
RUSER as supplied to F02EKF. You are free to use the arrays
IUSER and
RUSER to supply information to
MONIT as an alternative to using COMMON global variables.
MONIT must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which F02EKF is called. Parameters denoted as
Input must
not be changed by this procedure.
- 10: OPTION – SUBROUTINE, supplied by the NAG Library or the user.External Procedure
-
You can supply non-default options to the Arnoldi eigensolver by repeated calls to
F12ADF from within
OPTION. For example, you can set the mode to
Shifted Inverse Real, you can increase the
Iteration Limit beyond its default and you can print varying levels of detail on the iterative process using
Print Level.
If only the default options (including that the eigenvalues of largest magnitude are sought) are to be used then
OPTION may be the dummy subroutine F02EKY. (F02EKY is included in the NAG Library.) See
Section 9 for an example of using
OPTION to set some non-default options.
The specification of
OPTION is:
INTEGER |
ICOMM(*), ISTAT, IUSER(*) |
REAL (KIND=nag_wp) |
COMM(*), RUSER(*) |
|
- 1: ICOMM() – INTEGER arrayCommunication Array
On entry: contains details of the default option set. This array must be passed as parameter
ICOMM in any call to
F12ADF.
On exit: contains data on the current options set which may be altered from the default set via calls to
F12ADF.
- 2: COMM() – REAL (KIND=nag_wp) arrayCommunication Array
On entry: contains details of the default option set. This array must be passed as parameter
ICOMM in any call to
F12ADF.
On exit: contains data on the current options set which may be altered from the default set via calls to
F12ADF.
- 3: ISTAT – INTEGERInput/Output
On entry: set to zero.
On exit: if set to a nonzero value F02EKF returns immediately with .
- 4: IUSER() – INTEGER arrayUser Workspace
- 5: RUSER() – REAL (KIND=nag_wp) arrayUser Workspace
-
OPTION is called with the parameters
IUSER and
RUSER as supplied to F02EKF. You are free to use the arrays
IUSER and
RUSER to supply information to
OPTION as an alternative to using COMMON global variables.
OPTION must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which F02EKF is called.
- 11: NCONV – INTEGEROutput
On exit: the number of converged approximations to the selected eigenvalues. On successful exit, this will normally be either
NEV or
depending on the number of complex conjugate pairs of eigenvalues returned.
- 12: W(NCV) – COMPLEX (KIND=nag_wp) arrayOutput
On exit: the first
NCONV elements contain the converged approximations to the selected eigenvalues. Since complex conjugate pairs of eigenvalues appear together, it is possible (given an odd number of converged real eigenvalues) for F02EKF to return one more eigenvalue than requested.
- 13: V(LDV,) – REAL (KIND=nag_wp) arrayOutput
-
Note: the second dimension of the array
V
must be at least
.
On exit: contains the eigenvectors associated with the eigenvalue
, for
(stored in
W). For complex conjugate pairs of eigenvalues,
, the real and imaginary parts of the corresponding eigenvectors are stored, respectively, in
and
, for
; the second vector stored is for the first of the conjugate pair which must be negated for the second of the pair.
- 14: LDV – INTEGERInput
On entry: the first dimension of the array
V as declared in the (sub)program from which F02EKF is called.
Constraint:
.
- 15: RESID() – REAL (KIND=nag_wp) arrayOutput
-
On exit: the residual for the estimates to the eigenpair and is returned in , for .
- 16: IUSER() – INTEGER arrayUser Workspace
-
IUSER is not used by F02EKF, but is passed directly to
MONIT and
OPTION and may be used to pass information to these routines as an alternative to using COMMON global variables.
- 17: RUSER() – REAL (KIND=nag_wp) arrayUser Workspace
-
RUSER is not used by F02EKF, but is passed directly to
MONIT and
OPTION and may be used to pass information to these routines as an alternative to using COMMON global variables.
- 18: IFAIL – INTEGERInput/Output
-
On entry:
IFAIL must be set to
,
. If you are unfamiliar with this parameter you should refer to
Section 3.3 in the Essential Introduction 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
is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is
.
When the value 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, .
Constraint: .
-
On entry, .
Constraint: .
On entry, and .
Constraint: .
-
On entry, in shifted inverse mode, the th diagonal element of is not defined, for .
-
On entry, for , and .
Constraint: .
On entry, .
Constraint: .
On entry, and .
Constraint: .
-
On entry, in specification of column , and for
, and
.
Constraint: .
-
On entry, .
Constraint: .
-
On entry, and .
Constraint: .
On entry, and .
Constraint: .
-
On entry, the matrix is nearly numerically singular and could
not be inverted. Try perturbing the value of sigma.
Norm of matrix ,
Reciprocal condition number .
On entry, the matrix is numerically singular and could
not be inverted. Try perturbing the value of .
-
User requested termination in
MONIT,
.
-
User requested termination in
OPTION,
.
-
On entry, and .
Constraint: .
-
The maximum number of iterations
, the optional parameter
Iteration Limit has been set to
.
-
An internal call to
F12ABF returned with
.
This error should not occur. Please contact
NAG.
-
An internal call to
F12ABF returned with
.
-
The maximum number of iterations has been reached.
The maximum number of iterations .
The number of converged eigenvalues .
See the routine document for further details.
-
No shifts could be applied during a cycle of the
implicitly restarted Arnoldi iteration.
-
Could not build an Arnoldi factorization. The size
of the current Arnoldi factorization .
-
Error in internal call to compute eigenvalues and corresponding
error bounds of the current upper Hessenberg matrix.
Please contact
NAG.
-
An internal call to
F12ACF returned with
.
-
The number of eigenvalues found to sufficient accuracy,
is zero.
-
Internal inconsistency in the number of converged Ritz values number counted , number expected .
-
During calculation of a real Schur form, there was a failure to compute eigenvalues in a total of iterations.
-
The computed Schur form could not be reordered by an internal call.
This routine returned with
.
Please contact
NAG.
-
In calculating eigenvectors, an internal call returned with an error.
Please contact
NAG.
-
Dynamic memory allocation failed.
7 Accuracy
The relative accuracy of a Ritz value (eigenvalue approximation),
, is considered acceptable if its Ritz estimate
. The default value for
Tolerance is the
machine precision given by
X02AJF. The Ritz estimates are available via the
MONIT subroutine at each iteration, or can be printed by setting option
Print Level to a positive value.
F02EKF calls routines based on the ARPACK suite in
Chapter F12. These routines use an implicitly restarted Arnoldi iterative method to converge to approximations to a set of required eigenvalues (see the
F12 Chapter Introduction).
In the default
Regular mode, only matrix-vector multiplications are performed using the sparse matrix
during the Arnoldi process. Each iteration is therefore cheap computationally, relative to the alternative,
Shifted Inverse Real, mode described below. It is most efficient (i.e., the total number of iterations required is small) when the eigenvalues of largest magnitude are sought and these are distinct.
Although there is an option for returning the smallest eigenvalues using this mode (see
Smallest Magnitude option), the number of iterations required for convergence will be far greater or the method may not converge at all. However, where convergence is achieved,
Regular mode may still prove to be the most efficient since no inversions are required. Where smallest eigenvalues are sought and
Regular mode is not suitable, or eigenvalues close to a given real value are sought, the
Shifted Inverse Real mode should be used.
If the
Shifted Inverse Real mode is used (via a call to
F12ADF in
OPTION) then the matrix
is used in linear system solves by the Arnoldi process. This is first factorized internally using the direct
factorization routine
F11MEF. The condition number of
is then calculated by a call to
F11MGF. If the condition number is too big then the matrix is considered to be nearly singular, i.e.,
is an approximate eigenvalue of
, and the routine exits with an error. In this situation it is normally sufficient to perturb
by a small amount and call F02EKF again. After successful factorization, subsequent solves are performed by calls to
F11MFF.
Finally, F02EKF transforms the eigenvectors. Each eigenvector (real or complex) is normalized so that , and the element of largest absolute value is real and positive.
The monitoring routine
MONIT provides some basic information on the convergence of the Arnoldi iterations. Much greater levels of detail on the Arnoldi process are available via option
Print Level. If this is set to a positive value then information will be printed, by default, to standard output. The
Monitoring option may be used to select a channel number for the printing of monitoring information and
X04ACF may be used to associate that unit number with a file.
9 Example
This example computes the four eigenvalues of the matrix
which lie closest to the value
on the real line, and their corresponding eigenvectors, where
is the tridiagonal matrix with elements
9.1 Program Text
Program Text (f02ekfe.f90)
9.2 Program Data
Program Data (f02ekfe.d)
9.3 Program Results
Program Results (f02ekfe.r)
10 Optional Parameters
Internally F02EKF calls routines from the suite
F12AAF,
F12ABF,
F12ACF,
F12ADF and
F12AEF. Several optional parameters for these computational routines define choices in the problem specification or the algorithm logic. In order to reduce the number of formal parameters of F02EKF these optional parameters are also used here and have associated
default values that are usually appropriate. Therefore, you need only specify those optional parameters whose values are to be different from their default values.
Optional parameters may be specified via the user-supplied subroutine
OPTION in the call to F02EKF.
OPTION must be coded such that one call to
F12ADF is necessary to set each optional parameter. All optional parameters you do not specify are set to their default values.
The remainder of this section can be skipped if you wish to use the default values for all optional parameters.
The following is a list of the optional parameters available. A full description of each optional parameter is provided in
Section 10.1.
10.1 Description of the Optional Parameters
For each option, we give a summary line, a description of the optional parameter and details of constraints.
The summary line contains:
- the keywords, where the minimum abbreviation of each keyword is underlined;
- a parameter value,
where the letters , denote options that take character, integer and real values respectively;
- the default value, where the symbol is a generic notation for machine precision (see X02AJF).
Keywords and character values are case and white space insensitive.
Advisory | | Default the value returned by X04ABF
|
If the optional parameter
List is set then optional parameter specifications are printed on the advisory message channel number
.
This special keyword may be used to reset all optional parameters to their default values.
Iteration Limit | |
Default |
The limit on the number of Arnoldi iterations that can be performed before F02EKF exits with .
The Arnoldi iterative method converges on a number of eigenvalues with given properties. The default is to compute the eigenvalues of largest magnitude using
Largest Magnitude. Alternatively, eigenvalues may be chosen which have
Largest Real part,
Largest Imaginary part,
Smallest Magnitude,
Smallest Real part or
Smallest Imaginary part.
Note that these options select the eigenvalue properties for eigenvalues of the linear operator determined by the computational mode and problem type.
Normally each optional parameter specification is not printed to
the advisory channel
as it is supplied. Optional parameter
List may be used to enable printing and optional parameter
Nolist may be used to suppress the printing.
If
, monitoring information is output to channel number
during the solution of each problem; this may be the same as the
Advisory channel number. The type of information produced is dependent on the value of
Print Level, see the description of the optional parameter
Print Level for details of the information produced. Please see
X04ACF to associate a file with a given channel number.
This controls the amount of printing produced by F02EKF as follows.
| No output except error messages. |
| The set of selected options. |
| Problem and timing statistics when all calls to F12ABF have been completed. |
| A single line of summary output at each Arnoldi iteration. |
| If , then at each iteration, the length and additional steps of the current Arnoldi factorization and the number of converged Ritz values; during re-orthogonalization, the norm of initial/restarted starting vector. |
| Problem and timing statistics on final exit from F12ABF. If , then at each iteration, the number of shifts being applied, the eigenvalues and estimates of the Hessenberg matrix , the size of the Arnoldi basis, the wanted Ritz values and associated Ritz estimates and the shifts applied; vector norms prior to and following re-orthogonalization. |
| If , then on final iteration, the norm of the residual; when computing the Schur form, the eigenvalues and Ritz estimates both before and after sorting; for each iteration, the norm of residual for compressed factorization and the compressed upper Hessenberg matrix ; during re-orthogonalization, the initial/restarted starting vector; during the Arnoldi iteration loop, a restart is flagged and the number of the residual requiring iterative refinement; while applying shifts, the indices of the shifts being applied. |
| If , then during the Arnoldi iteration loop, the Arnoldi vector number and norm of the current residual; while applying shifts, key measures of progress and the order of ; while computing eigenvalues of , the last rows of the Schur and eigenvector matrices; when computing implicit shifts, the eigenvalues and Ritz estimates of . |
| If , then during Arnoldi iteration loop: norms of key components and the active column of , norms of residuals during iterative refinement, the final upper Hessenberg matrix ; while applying shifts: number of shifts, shift values, block indices, updated matrix ; while computing eigenvalues of : the matrix , the computed eigenvalues and Ritz estimates. |
These options define the computational mode which in turn defines the form of operation to be performed.
Given a standard eigenvalue problem in the form
then the following modes are available with the appropriate operator
.
An approximate eigenvalue has deemed to have converged when the corresponding Ritz estimate is within
Tolerance relative to the magnitude of the eigenvalue.