spral_ssmfe
 Sparse Symmetric MatrixFree Eigensolver¶
Purpose¶
This package computes extreme (leftmost and/or rightmost) eigenpairs \(\{\lambda_i, x_i\}\) of the following eigenvalue problems:
the standard eigenvalue problem
\[A x = \lambda x,\]the generalized eigenvalue problem
\[A x = \lambda B x,\]the buckling problem
\[B x = \lambda A x,\]
where \(A\) and \(B\) are real symmetric (or Hermitian) matrices and \(B\) is positive definite.
This package provides a userfriendly wrapper around
spral_ssmfe_expert
, which in turn provides a wrapper around
spral_ssmfe_core
. If more finetuned control of the eigensolver is
required, use those modules instead.
Usage overview¶
The eigensolver subroutines behind spral_ssmfe
implement a block
iterative algorithm. The block nature of this algorithm allows the user
to benefit from highly optimized linear algebra subroutines and from the
ubiquitous multicore architecture of modern computers. It also makes
this algorithm more reliable than Krylovbased algorithms employed e.g.
by ARPACK in the presence of clustered eigenvalues. However, convergence
of the iterations may be slow if the density of the spectrum is high.
Thus, good performance (in terms of speed) is contingent on the following two factors:
 the number of desired eigenpairs must be substantial (e.g. not fewer than the number of CPU cores), and
 the employment of a convergence acceleration technique.
The acceleration techniques that can be used are shiftandinvert and preconditioning.
The former requires the direct solution of linear systems with the matrix \(A\) or its linear combination with \(B\), for which a sparse symmetric indefinite solver (such as HSL_MA97 or SPRAL_SSIDS) can be employed.
The latter applies to the case of positive definite \(A\) and requires a matrix or an operator \(T\), called a preconditioner, such that the vector \(v = T f\) is an approximation to the solution \(u\) of the system \(A u = f\) (see the simple example below). Note: This technique is only recommended for experienced users.
Subroutines¶

subroutine
ssmfe_standard
(rci, left, mep, lambda, n, x, ldx, keep, options, inform)¶ Computes the leftmost eigenpairs of the standard eigenvalue problem
\[Ax = \lambda x\]Optionally uses preconditioning.
Uses reversecommunication. Upon return the user must perform a task specified by the rci parameter and recall the routine. Possible values of rci and associated tasks are:
rci%job Task to be performed 3 None. Fatal error, see inform%flag. 2 Restart computation. Nonfatal error, see inform%flag. 1 None. Computation complete. 1 Calculate \(Y = AX\). 2 Apply preconditioner \(Y = TX\). The matrices \(X\) and \(Y\) are components of rci.
Parameters:  rci [ssmfe_rcid,inout] :: Reverse communication type. rci%job must be
set to 0 before the first call. (Type
ssmfe_rciz
in complex version).  left [integer,in] :: Number of left eigenpairs to find.
 mep [integer,in] :: Number of working eigenpairs. See method section for guidance on selecting a good value. Must be at least left.
 lambda (mep) [real,inout] :: Current eigenvalue estimates in ascending order.
 n [integer,in] :: Size of matrix \(A\).
 x (ldx, n) [real,inout] :: Current eigenvector estimates corresponding to eigenvalues in lambda. Used to supply initial estimates if options%user_x>0. (Type complex in complex version).
 ldx [integer,in] :: Leading dimension of x.
 keep [ssmfe_keepd,inout] :: Internal workspace used by routine. (Type
ssmfe_keepz
in complex version).  options [ssmfe_options,in] :: specifies algorithm options to be used.
 inform [ssmfe_inform,inout] :: returns information about the exection of the routine.
 rci [ssmfe_rcid,inout] :: Reverse communication type. rci%job must be
set to 0 before the first call. (Type

subroutine
ssmfe_standard_shift
(rci, sigma, left, right, mep, lambda, n, x, ldx, keep, options, inform)¶ Computes eigenpairs of the standard eigenvalue problem
\[Ax = \lambda x\]in the vicinity of a given value \(\sigma\).
Uses reversecommunication. Upon return the user must perform a task specified by the rci parameter and recall the routine. Possible values of rci and associated tasks are:
rci%job Task to be performed 3 None. Fatal error, see inform%flag. 2 Restart computation. Nonfatal error, see inform%flag. 1 None. Computation complete. 1 Calculate \(Y = AX\). 9 Solve \((A\sigma I)Y = X\) for Y. The matrices \(X\) and \(Y\) are components of rci.
Parameters:  rci [ssmfe_rcid,inout] :: Reverse communication type. rci%job must be
set to 0 before the first call. (Type
ssmfe_rciz
in complex version).  sigma [real,in] :: Shift value \(sigma\).
 left [integer,in] :: Number of left eigenpairs to find.
 right [integer,in] :: Number of right eigenpairs to find.
 mep [integer,in] :: Number of working eigenpairs. See method section for guidance on selecting a good value. Must be at least left+right.
 lambda (mep) [real,inout] :: Current eigenvalue estimates in ascending order.
 n [integer,in] :: Size of matrix \(A\).
 x (ldx, n) [real,inout] :: Current eigenvector estimates corresponding to eigenvalues in lambda. Used to supply initial estimates if options%user_x>0. (Type complex in complex version).
 ldx [integer,in] :: Leading dimension of x.
 keep [ssmfe_keepd,inout] :: Internal workspace used by routine. (Type
ssmfe_keepz
in complex version).  options [ssmfe_options,in] :: specifies algorithm options to be used.
 inform [ssmfe_inform,inout] :: returns information about the exection of the routine.
 rci [ssmfe_rcid,inout] :: Reverse communication type. rci%job must be
set to 0 before the first call. (Type

subroutine
ssmfe_generalized
(rci, left, mep, lambda, n, x, ldx, keep, options, inform)¶ Computes the leftmost eigenpairs of the generalized eigenvalue problem
\[Ax = \lambda B x\]Optionally uses preconditioning.
Uses reversecommunication. Upon return the user must perform a task specified by the rci parameter and recall the routine. Possible values of rci and associated tasks are:
rci%job Task to be performed 3 None. Fatal error, see inform%flag. 2 Restart computation. Nonfatal error, see inform%flag. 1 None. Computation complete. 1 Calculate \(Y = AX\). 2 Apply preconditioner \(Y = TX\). 3 Calculate \(Y = BX\). The matrices \(X\) and \(Y\) are components of rci.
Parameters:  rci [ssmfe_rcid,inout] :: Reverse communication type. rci%job must be
set to 0 before the first call. (Type
ssmfe_rciz
in complex version).  left [integer,in] :: Number of left eigenpairs to find.
 mep [integer,in] :: Number of working eigenpairs. See method section for guidance on selecting a good value. Must be at least left.
 lambda (mep) [real,inout] :: Current eigenvalue estimates in ascending order.
 n [integer,in] :: Size of matrix \(A\).
 x (ldx, n) [real,inout] :: Current eigenvector estimates corresponding to eigenvalues in lambda. Used to supply initial estimates if options%user_x>0. (Type complex in complex version).
 ldx [integer,in] :: Leading dimension of x.
 keep [ssmfe_keepd,inout] :: Internal workspace used by routine. (Type
ssmfe_keepz
in complex version).  options [ssmfe_options,in] :: specifies algorithm options to be used.
 inform [ssmfe_inform,inout] :: returns information about the exection of the routine.
 rci [ssmfe_rcid,inout] :: Reverse communication type. rci%job must be
set to 0 before the first call. (Type

subroutine
ssmfe_generalized_shift
(rci, sigma, left, right, mep, lambda, n, x, ldx, keep, options, inform)¶ Computes eigenpairs of the generalized eigenvalue problem
\[Ax = \lambda B x\]in the vicinity of a given value \(\sigma\).
Uses reversecommunication. Upon return the user must perform a task specified by the rci parameter and recall the routine. Possible values of rci and associated tasks are:
rci%job Task to be performed 3 None. Fatal error, see inform%flag. 2 Restart computation. Nonfatal error, see inform%flag. 1 None. Computation complete. 1 Calculate \(Y = AX\). 3 Calculate \(Y = BX\). 9 Solve \((A\sigma B)Y = X\) for Y. The matrices \(X\) and \(Y\) are components of rci.
Parameters:  rci [ssmfe_rcid,inout] :: Reverse communication type. rci%job must be
set to 0 before the first call. (Type
ssmfe_rciz
in complex version).  sigma [real,in] :: Shift value \(sigma\).
 left [integer,in] :: Number of left eigenpairs to find.
 right [integer,in] :: Number of right eigenpairs to find.
 mep [integer,in] :: Number of working eigenpairs. See method section for guidance on selecting a good value. Must be at least left+right.
 lambda (mep) [real,inout] :: Current eigenvalue estimates in ascending order.
 n [integer,in] :: Size of matrix \(A\).
 x (ldx, n) [real,inout] :: Current eigenvector estimates corresponding to eigenvalues in lambda. Used to supply initial estimates if options%user_x>0. (Type complex in complex version).
 ldx [integer,in] :: Leading dimension of x.
 keep [ssmfe_keepd,inout] :: Internal workspace used by routine. (Type
ssmfe_keepz
in complex version).  options [ssmfe_options,in] :: specifies algorithm options to be used.
 inform [ssmfe_inform,inout] :: returns information about the exection of the routine.
 rci [ssmfe_rcid,inout] :: Reverse communication type. rci%job must be
set to 0 before the first call. (Type

subroutine
ssmfe_buckling
(rci, sigma, left, right, mep, lambda, n, x, ldx, keep, options, inform)¶ Computes eigenpairs of the buckling problem
\[Bx = \lambda A x\]in the vicinity of a given value \(\sigma\).
Uses reversecommunication. Upon return the user must perform a task specified by the rci parameter and recall the routine. Possible values of rci and associated tasks are:
rci%job Task to be performed 3 None. Fatal error, see inform%flag. 2 Restart computation. Nonfatal error, see inform%flag. 1 None. Computation complete. 1 Calculate \(Y = AX\). 3 Calculate \(Y = BX\). 9 Solve \((B\sigma A)Y = X\) for Y. The matrices \(X\) and \(Y\) are components of rci.
Parameters:  rci [ssmfe_rcid,inout] :: Reverse communication type. rci%job must be
set to 0 before the first call. (Type
ssmfe_rciz
in complex version).  sigma [real,in] :: Shift value \(sigma\).
 left [integer,in] :: Number of left eigenpairs to find.
 right [integer,in] :: Number of right eigenpairs to find.
 mep [integer,in] :: Number of working eigenpairs. See method section for guidance on selecting a good value. Must be at least left+right.
 lambda (mep) [real,inout] :: Current eigenvalue estimates in ascending order.
 n [integer,in] :: Size of matrix \(A\).
 x (ldx, n) [real,inout] :: Current eigenvector estimates corresponding to eigenvalues in lambda. Used to supply initial estimates if options%user_x>0. (Type complex in complex version).
 ldx [integer,in] :: Leading dimension of x.
 keep [ssmfe_keepd,inout] :: Internal workspace used by routine. (Type
ssmfe_keepz
in complex version).  options [ssmfe_options,in] :: specifies algorithm options to be used.
 inform [ssmfe_inform,inout] :: returns information about the exection of the routine.
 rci [ssmfe_rcid,inout] :: Reverse communication type. rci%job must be
set to 0 before the first call. (Type

subroutine
ssmfe_free
(keep, inform)¶ Free memory allocated in keep and inform. Unnecessary if both are going out of scope.
Parameters:  keep [ssmfe_keepd,inout] :: Workspace to be freed.
(Type
ssmfe_keepz
in complex version).  inform [ssmfe_inform,inout] :: Information type to be freed.
 keep [ssmfe_keepd,inout] :: Workspace to be freed.
(Type
Derived types¶

type
ssmfe_rcid
¶ Realvalued reverse communication interface (RCI) type.
Type fields: %
job
[integer] :: Reversecommunication task to perform.%
nx
[integer] :: Number of columns in x and y.%
x
(n, nx) [real] :: Vector to be transformed by RCI task.%
y
(n, nx) [real] :: Vector to store result of RCI task.

type
ssmfe_rciz
¶ Complexvalued reverse communication interface (RCI) type.
Type fields: %
job
[integer] :: Reversecommunication task to perform.%
nx
[integer] :: Number of columns in x and y.%
x
(n, nx) [complex] :: Vector to be transformed by RCI task.%
y
(n, nx) [complex] :: Vector to store result of RCI task.

type
ssmfe_options
¶ Options that control the algorithm.
Type fields: %
abs_tol_lambda
[real,default=0.0] :: absolute tolerance for estimated eigenvalue convergence test, see Section [ssmfe:method]. Negative values are treated as the default.%
abs_tol_residual
[real,default=0.0] :: absolute tolerance for residual convergence test, see Section [ssmfe:method]. Negative values are treated as the default.%
max_iterations
[integer,default=100] :: maximum number of iterations.%
rel_tol_lambda
[real,default=0.0] :: relative tolerance for estimated eigenvalue error convergence test, see Section [ssmfe:method]. Negative values are treated as the default.%
rel_tol_residual
[real,default=0.0] :: relative tolerance for residual convergence test, see Section [ssmfe:method]. If both abs_tol_residual and rel_tol_residual are 0.0, then the residual norms are not taken into consideration by the convergence test, see Section [ssmfe:method]. Negative values are treated as the default.%
tol_x
[real,default=1.0] :: tolerance for estimated eigenvector error convergence test, see Section [ssmfe:method]. If tol_x is set to 0.0, the eigenvector error is not estimated. If a negative value is assigned, the tolerance is set to sqrt(epsilon(lambda)).%
print_level
[integer,default=0] ::amount of printing. Possible values are:
<0 no printing 0 error and warning messages only 1 the type (standard or generalized) and the size of the problem, the number of eigenpairs requested, the error tolerances and the size of the subspace are printed before the iterations start 2 as above but, for each eigenpair tested for convergence, the iteration number, the index of the eigenpair, the eigenvalue, whether it has converged, the residual norm, and the error estimates are also printed >2 as 1 but with all eigenvalues, whether converged, residual norms and eigenvalue/eigenvector error estimates printed on each iteration. Note that for eigenpairs that are far from convergence, ‘rough’ error estimates are printed (the estimates that are actually used by the stopping criteria, see Section [ssmfe:method], only become available on the last few iterations).
%
unit_error
[integer,default=6] :: unit number for error messages. Printing suppressed if negative.%
unit_diagnostic
[integer,default=6] :: unit number for diagnostic messages. Printing suppressed if negative.%
unit_warning
[integer,default=6] :: unit number for warning messages. Printing suppressed if negative.%
left_gap
[real,default=0.0] :: minimal acceptable distance between last computed left eigenvalue and rest of spectrum. Forssmfe_standard()
andssmfe_generalized()
the last computed left eigenvalue is the rightmost of those computed. For other routines it is the leftmost. If set to a negative value \(\delta\), the minimal distance is taken as \(\delta\) times the average distance between the computed eigenvalues. Note that for this option to have any effect, the value of mep must be larger than left+right. See Section [ssmfe:method] for further explanation.%
max_left
[integer,default=1] :: number of eigenvalues to left of \(\sigma\), or a negative value if not known.%
max_right
[integer,default=1] :: number of eigenvalues to right of \(\sigma\), or a negative value if not known.%
right_gap
[real,default=0.0] :: as left_gap, but for right eigenvalues.%
user_x
[integer,default=0] :: number of eigenvectors for which an initial guess is supplied in x(:,:) on the first call. Such eigenvectors must be lineraly independent.

type
ssmfe_inform
¶ Information on progress of the algorithm.
Type fields: %
flag
[integer] :: return status of algorithm. See table below.%
iteration
[integer] :: number of iterations.%
left
[integer] :: number of converged left eigenvalues.%
next_left
[real] :: upon completion, next left eigenvalue in spectrum (see options%left_gap).%
next_right
[real] :: upon completion, next right eigenvalue in spectrum (see options%right_gap).%
non_converged
[integer] :: number of nonconverged eigenpairs.%
right
[integer] :: number of converged right eigenvalues.%
stat
[integer] :: allocation status in event of failure
inform%flag 1 rci%job is outofrange. 9 n is outofrange. 10 ldx is outofrange. 11 left is outofrange. 12 right is outofrange. 13 mep is less than the number of desired eigenpairs. 100 Not enough memory; inform%stat contains the value of the Fortran stat parameter. 200 \(B\) is not positive definite or user_x>0 and linearly dependent initial guesses were supplied. +1 The iterations have been terminated because no further improvement in accuracy is possible (this may happen if \(B\) or the preconditioner is not positive definite, or if the components of the residual vectors are so small that the roundoff errors make them essentially random). The value of inform%non_converged is set to the number of nonconverged eigenpairs. +2 The maximum number of iterations max_iterations has been exceeded. The value of inform%non_converged is set to the number of nonconverged eigenpairs. +3 The solver had run out of storage space for the converged eigenpairs before the gap in the spectrum required by options%left_gap and/or options%right_gap was reached. The value of inform%non_converged is set to the number of nonconverged eigenpairs. If the computation is terminated with the error code 2 or 3, the computation is not complete, but may be restarted with larger values of max_iterations and/or mep. In this case the user should set options%user_x to info%left + info%right and restart the reverse communication loop. An alternative option is to use one of the advanced solver procedures from
spral_ssmfe_expert
orspral_ssmfe_core
that delegate the storage of computed eigenpairs and the termination of the computation to the user.
Examples¶
Preconditioning example¶
The following code computes the 5 leftmost eigenpairs of the matrix \(A\) of order 100 that approximates the twodimensional Laplacian operator on a 20by20 grid. One forward and one backward GaussSeidel update are used for preconditioning, which halves the number of iterations compared with solving the same problem without preconditioning. The module laplace2d (examples/Fortran/ssmfe/laplace2d.f90) supplies a subroutine apply_laplacian() that multiplies a block of vectors by \(A\), and a subroutine apply_gauss_seidel_step() that computes \(y = T x\) for a block of vectors \(x\) by applying one forward and one backward update of the GaussSeidel method to the system \(A y = x\).
! examples/Fortran/ssmfe/precond_ssmfe.f90
! Laplacian on a square grid (using SPRAL_SSMFE routines)
program ssmfe_precond_example
use spral_ssmfe
use laplace2d ! implement Lapalacian and preconditioners
implicit none
integer, parameter :: wp = kind(0d0) ! Working precision is double
integer, parameter :: m = 20 ! grid points along each side
integer, parameter :: n = m*m ! problem size
integer, parameter :: nep = 5 ! eigenpairs wanted
real(wp) :: lambda(2*nep) ! eigenvalues
real(wp) :: X(n, 2*nep) ! eigenvectors
type(ssmfe_rcid ) :: rci ! reverse communication data
type(ssmfe_options) :: options ! options
type(ssmfe_keepd ) :: keep ! private data
type(ssmfe_inform ) :: inform ! information
integer :: i ! loop index
! the gap between the last converged eigenvalue and the rest of the spectrum
! must be at least 0.1 times average gap between computed eigenvalues
options%left_gap = 0.1
rci%job = 0
do ! reverse communication loop
call ssmfe_standard &
( rci, nep, 2*nep, lambda, n, X, n, keep, options, inform )
select case ( rci%job )
case ( 1 )
call apply_laplacian( m, m, rci%nx, rci%x, rci%y )
case ( 2 )
call apply_gauss_seidel_step( m, m, rci%nx, rci%x, rci%y )
case ( :1 )
exit
end select
end do
print '(i3, 1x, a, i3, 1x, a)', inform%left, 'eigenpairs converged in', &
inform%iteration, 'iterations'
print '(1x, a, i2, a, es13.7)', &
('lambda(', i, ') = ', lambda(i), i = 1, inform%left)
call ssmfe_free( keep, inform )
end program ssmfe_precond_example
This code produces the following output:
6 eigenpairs converged in 19 iterations
lambda( 1) = 4.4676695E02
lambda( 2) = 1.1119274E01
lambda( 3) = 1.1119274E01
lambda( 4) = 1.7770878E01
lambda( 5) = 2.2040061E01
lambda( 6) = 2.2040061E01
Note that the code computed one extra eigenpair because of the insufficient gap between the 5th and 6th eigenvalues.
Shiftandinvert example¶
The following code computes the eigenpairs of the matrix of order 64 that approximates the twodimensional Laplacian operator on 8by8 grid with eigenvalues near the shift sigma=1.0. For the shifted solve, LAPACK subroutines DSYTRS and DSYTRF are used, which perform the LDLTfactorization and the solution of the factorized system respectively. The matrix of the discretized Laplacian is computed by the subroutine set_2d_laplacian_matrix() from the laplace2d module (examples/Fortran/ssmfe/laplace2d.f90). The module ldltf (examples/Fortran/ssmfe/ldltf.f90) supplies the function num_neg_D() that counts the number of negative eigenvalues of the Dfactor.
! examples/Fortran/ssmfe/shift_invert.f90
! Laplacian on a rectangular grid by shiftinvert via LDLT factorization
program ssmfe_shift_invert_example
use spral_ssmfe
use laplace2d ! implement Lapalacian and preconditioners
use ldltf ! implements LDLT support routines
implicit none
integer, parameter :: wp = kind(0d0) ! Working precision
integer, parameter :: nx = 8 ! grid points along x
integer, parameter :: ny = 8 ! grid points along y
integer, parameter :: n = nx*ny ! problem size
real(wp), parameter :: sigma = 1.0 ! shift
integer :: ipiv(n) ! LDLT pivot index
real(wp) :: lambda(n) ! eigenvalues
real(wp) :: X(n, n) ! eigenvectors
real(wp) :: A(n, n) ! matrix
real(wp) :: LDLT(n, n) ! factors
real(wp) :: work(n*n) ! work array for dsytrf
integer :: lwork = n*n ! size of work
integer :: left, right ! wanted eigenvalues left and right of sigma
integer :: i ! index
type(ssmfe_options) :: options ! eigensolver options
type(ssmfe_inform ) :: inform ! information
type(ssmfe_rcid ) :: rci ! reverse communication data
type(ssmfe_keepd ) :: keep ! private data
call set_laplacian_matrix( nx, ny, A, n )
! perform LDLT factorization of the shifted matrix
LDLT = A
forall ( i = 1 : n ) LDLT(i, i) = A(i, i)  sigma
lwork = n*n
call dsytrf( 'L', n, LDLT, n, ipiv, work, lwork, i )
left = num_neg_D(n, LDLT, n, ipiv) ! all eigenvalues to the left from sigma
right = 5 ! 5 eigenvalues to the right from sigma
rci%job = 0
do
call ssmfe_standard_shift &
( rci, sigma, left, right, n, lambda, n, X, n, keep, options, inform )
select case ( rci%job )
case ( 1 )
call dgemm &
( 'N', 'N', n, rci%nx, n, 1.0_wp, A, n, rci%x, n, 0.0_wp, rci%y, n )
case ( 9 )
call dcopy( n * rci%nx, rci%x, 1, rci%y, 1 )
call dsytrs( 'L', n, rci%nx, LDLT, n, ipiv, rci%y, n, i )
case ( :1 )
exit
end select
end do
print '(1x, a, es10.2, 1x, a, i3, 1x, a)', 'Eigenvalues near', sigma, &
'(took', inform%iteration, 'iterations)'
print '(1x, a, i2, a, es13.7)', &
('lambda(', i, ') = ', lambda(i), i = 1, inform%left + inform%right)
call ssmfe_free( keep, inform )
end program ssmfe_shift_invert_example
This code produces the following output:
Eigenvalues near 1.00E+00 (took 5 iterations)
lambda( 1) = 2.4122952E01
lambda( 2) = 5.8852587E01
lambda( 3) = 5.8852587E01
lambda( 4) = 9.3582223E01
lambda( 5) = 1.1206148E+00
lambda( 6) = 1.1206148E+00
lambda( 7) = 1.4679111E+00
lambda( 8) = 1.4679111E+00
lambda( 9) = 1.7733184E+00
Hermitian example¶
The following code computes the 5 leftmost eigenpairs of the differential operator \(i \frac{d}{dx}\) acting in the space of periodic functions discretized by central differences on a uniform mesh of 80 steps.
! examples/Fortran/ssmfe/hermitian.f90  Example code for SPRAL_SSMFE package
! Hermitian operator example
program ssmfe_hermitian_example
use spral_ssmfe
implicit none
integer, parameter :: wp = kind(0d0) ! working precision
integer, parameter :: n = 80 ! problem size
integer, parameter :: nep = 5 ! eigenpairs wanted
real(wp) :: lambda(nep) ! eigenvalues
complex(wp) :: X(n, nep) ! eigenvectors
type(ssmfe_rciz ) :: rci ! reverse communication data
type(ssmfe_options) :: options ! options
type(ssmfe_keepz ) :: keep ! private data
type(ssmfe_inform ) :: inform ! information
integer :: i ! loop index
rci%job = 0
do ! reverse communication loop
call ssmfe_standard( rci, nep, nep, lambda, n, X, n, keep, options, inform )
select case ( rci%job )
case ( 1 )
call apply_idx( n, rci%nx, rci%x, rci%y )
case ( :1 )
exit
end select
end do
print '(i3, 1x, a, i3, 1x, a)', inform%left, 'eigenpairs converged in', &
inform%iteration, 'iterations'
print '(1x, a, i2, a, es14.7)', &
('lambda(', i, ') = ', lambda(i), i = 1, inform%left)
call ssmfe_free( keep, inform )
contains
subroutine apply_idx( n, m, x, y ) ! central differences for i d/dx
implicit none
complex(wp), parameter :: IM_ONE = (0.0D0, 1.0D0)
integer, intent(in) :: n, m
complex(wp), intent(in) :: x(n, m)
complex(wp), intent(out) :: y(n, m)
integer :: i, j, il, ir
do j = 1, m
do i = 1, n
if ( i == 1 ) then
il = n
else
il = i  1
end if
if ( i == n ) then
ir = 1
else
ir = i + 1
end if
y(i, j) = IM_ONE*(x(ir, j)  x(il, j))
end do
end do
end subroutine apply_idx
end program ssmfe_hermitian_example
This code produces the following output:
5 eigenpairs converged in 25 iterations
lambda( 1) = 2.0000000E+00
lambda( 2) = 1.9938347E+00
lambda( 3) = 1.9938347E+00
lambda( 4) = 1.9753767E+00
lambda( 5) = 1.9753767E+00
Method¶
spral_ssmfe_core`
, upon which spral_ssmfe
is built,
implements a block iterative algorithm based on the Jacobiconjugate
preconditioned gradients (JCPG) method [2], [3]. This algorithm
simultaneously computes \(m < n\) approximate eigenpairs, where the block
size \(m\) exceeds the number \(n_e\) of desired eigenpairs for the
sake of better convergence, namely, \(m = n_e + \min(10, 0.1 n_e)\).
An approximate eigenpair \(\{x,\lambda\}\) is considered to have converged if the following three conditions are all satisfied:
 if options%abs_tol_lambda and options%rel_tol_lambda are not both equal to zero, then the estimated error in the approximate eigenvalue must be less than or equal to \(\max(\mathrm{options\%abs\_tol\_lambda}, \delta*\mathrm{options\%rel\_tol\_lambda})\), where \(\delta\) is the estimated average distance between eigenvalues.
 if options%tol_x is not zero, then the estimated sine of the angle between the approximate eigenvector and the invariant subspace corresponding to the eigenvalue approximated by \(\lambda\) must be less than or equal to options%tol_x.
 if options%abs_tol_residual and options%rel_tol_residual are not both equal to zero, then the Euclidean norm of the residual, \(\A x  \lambda B x\_2\), must be less than or equal to \(\max(\mathrm{options\%abs\_tol\_residual}, \mathrm{options\%rel\_tol\_residual}*\\lambda B x\_2)\).
The extra eigenpairs are not checked for convergence, as their role is purely auxiliary.
If the gap between the last computed eigenvalue and the rest of the
spectrum is small, then the accuracy of the corresponding eigenvector
may be very low. To prevent this from happening, the user should set the
eigenpairs storage size mep to a value that is larger than the number of
desired eigenpairs, and set the options options%left_gap and
options%right_gap to nonzero values \(\delta_l\) and
\(\delta_r\). These values determine the size of the minimal
acceptable gaps between the computed eigenvalues and the rest of the
spectrum, \(\delta_l\) referring to either leftmost eigenvalues (for
ssmfe_standard()
and ssmfe_generalized()
only) or those
to the left of the shift sigma, and \(\delta_r\) to those to the right of
the shift sigma. Positive values of \(\delta_l\) and \(\delta_r\) set
the gap explicitely, and negative values require the gap to be not less than
their absolute value times the average distance between the computed
eigenvalues. A recommended value of \(\delta_l\) and
\(\delta_r\) is 0.1. The value of mep has little effect on
the speed of computation, hence it might be set to any reasonably large
value. The larger the value of mep, the larger the size of an eigenvalue
cluster for which accurate eigenvectors can be computed, notably: to
safeguard against clusters of size up to \(k\), it is sufficient to
set mep to the number of desired eigenpairs plus \(k  1\).
When using the solver procedures that employ the shiftandinvert technique, it is very important to ensure that the numbers of desired eigenvalues each side of the shift do not exceed the actual numbers of these eigenvalues, as the eigenpairs ‘approximating’ nonexisting eigenpairs of the problem will not converge. It is therefore strongly recommended that the user employs a linear system solver that performs the \(LDL^T\) factorization of the shifted system, e.g. HSL_MA97 or SPRAL_SSIDS. The \(LDL^T\) factorization of the matrix \(A  \sigma B\) consists in finding a lower triangular matrix \(L\), a blockdiagonal matrix \(D\) with \(1\times 1\) and \(2\times 2\) blocks on the diagonal and a permutation matrix \(P\) such that \(P^T(A  \sigma B)P = L D L^T\). By the inertia theorem, the number of eigenvalues to the left and right from the shift \(\sigma\) is equal to the number of negative and positive eigenvalues of \(D\), which allows quick computation of the eigenvalue numbers each side of the shift.
References¶
[1]  E. E. Ovtchinnikov and J. Reid (2010). A preconditioned block conjugate gradient algorithm for computing extreme eigenpairs of symmetric and Hermitian problems. Technical Report RALTR201019. 
[2]  E. E. Ovtchinnikov (2008). Jacobi correction equation, line search and conjugate gradients in Hermitian eigenvalue computation I: Computing an extreme eigenvalue. SIAM J. Numer. Anal., 46:2567–2592. 
[3]  E. E. Ovtchinnikov (2008). Jacobi correction equation, line search and conjugate gradients in Hermitian eigenvalue computation II: Computing several extreme eigenvalues. SIAM J. Numer. Anal., 46:2593–2619. 