spral_ssmfe_core - Sparse Symmetric Matrix-Free Eigensolver (Core Algorithm)¶

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 eigenvalue problem

  \[AB x = \lambda x\]

where \(A\) and \(B\) are real symmetric (or Hermitian) matrices and \(B\) is positive definite.

The modules spral_ssmfe and spral_ssmfe_expert provide a simplified interface to this routine, and should be used if the user does not require access to low level features provided in this package.

Usage overview¶

spral_ssmfe_core implements a block iterative algorithm for simultaneous computation of several eigenpairs for the problems above. 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 Krylov-based algorithms employed by e.g. 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:

1. the number of desired eigenpairs must be substantial (e.g. not fewer than the number of CPU cores), and
2. the employment of a convergence acceleration technique.

The acceleration techniques that can be used are shift-and-invert and preconditioning.

The former rewrites the eigenvalue problem for a matrix \(M\) as \(Ax = \lambda x\) with \(A = (M - \sigma I)^{-1}\), where \(I\) is the identity matrix and \(\sigma\) is a real value near eigenvalues of interest, and the generalized problem \(M x = \mu B x\) as the \(Ax = \lambda B x\) with \(A = (M - \sigma B)^{-1}\).

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.

For futher detail on the algorithm, see the outline in the method section below, and associated references.

The routine spral_ssmfe provides a user-friendly interface to this algorithm, whilst spral_ssmfe_expert provides an interface that allows users to manage their own memory. If this routine is used instead of spral_ssmfe_expert, the user is additionally responsible for deciding when a ufficient number of eigenpairs have been computed to sufficient accuracy. The amount of computation performed by the solver subroutines in this package and the memory they use are negligible. These features facilitate the use of these subroutines for shared-memory, out-of-core and hybrid computation.

Subroutines¶

To use the solver procedures, the user must maintain a workspace of (kw+1) blocks each containing m vectors of size n. For notational convienience we refer to this workspace as a Fortran array W(n,m,0:kw), but the user is free to store it as they wish. Note the block dimension is indexed from zero, not from one. The following table provides minimum values of kw for each setup:

	minAprod=T		minAprod=F
Problem	minBprod=T	minBprod=F	minBprod=T	minBprod=F
standard	7	5	3	3
standard_shift	7	5	N/A	N/A
generalized	7	5	5	3
generalized_shift	7	7	N/A	N/A
buckling	7	7	N/A	N/A

Further, the user must also store the converged eigenvectors \(X\), and (for generalised problems) their \(B\)-images \(BX\) using separate storage, e.g. X(n,mep), BX(n,mep). In addition to being output, the routine may need to reorthagonalise against these from time to time.

The first (zero-indexed) block holds the eigenvector approximations: the user must fill this block with \(m\) linearly independent vectors before the first call to a solver procedure.

The number of desired eigenpairs may exceed \(m\): whenever converged eigenpairs have been detected, a solver procedure reports the indices of these eigenpairs and they must be moved by the user to separate storage (X(:)).

subroutine ssmfe(rci, problem, left, right, m, lambda, rr, ind, keep, options, inform)¶

Computes specified number of leftmost and rightmost eigenvalues and corresponding eigenvectors.

Uses reverse-communication. 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	Failed to converge, see inform%flag.
-1	None. Computation complete.
1	Calculate \(\bar{V} = AU\).
2	Apply preconditioner \(\bar{V} = TU\). (Copy if T=I).
3	Compute \(\bar{V} = BU\)
4	Test convergence for each of \(m\) eigenvalues: If eigenpair i has converged, set inform%converged(i) to a positive number. Otherwise, leave at current value. Tests may use the estimated eigenvalue errors inform%err_lambda(i) and eigenvector errors inform%err_x(i).
5	Copy converged eigenvectors \(X\) to user storage: If rci%i>0: W(:,rci%jx:rci%jx+rci%nx-1,rci%kx). Else: W(:,rci%jx-rci%nx+1:rci%jx,rci%kx). Optionally save their \(B\)-images (if \(B\ne I\)) If rci%i>0: W(:,rci%jx:rci%jx+rci%nx-1,rci%ky). Else: W(:,rci%jx-rci%nx+1:rci%jx,rci%ky).
11	If rci%i.eq.0, copy \(\bar{V} = U\). Otherwise, reorder columns of block rci%kx such that column ind(j) becomes the new column j for j=1, ..., rci%nx Note: if rci%kx.eq.rci%ky, only reorder once.
12	Compute the dot products \[r_{ii} = U_i \cdot \bar{V}_i\]
13	Perform the scalings \[U_i = U_i/\sqrt{(U_i\cdot \bar{V}_i)}\] and \[\bar{V}_i = \bar{V}_i/\sqrt{(U_i\cdot \bar{V}_i)}\] for each column \(U_i\) and \(\bar{V}_i\) of \(U\) and \(\bar{V}\). Note: if rci%kx.eq.rci%ky, only scale once.
14	Perform the updates \[\bar{V}_i = \bar{V}_i + r_{ii} U_i\] for each column \(\bar{V}_i\) of \(\bar{V}\)
15	Perform the update \[R = \alpha U^T V + \beta R\]
16	Perform the update \[V = \alpha U R + \beta V\]
17	Perform the update \[U = \alpha U R\] Note: \(V\) may be used as a workspace
21	\(B\)-orthogonalize columns of \(V\) to all vectors \(X\) by solving \[(X^TBX) Q = X^T \bar{V}\] for \(Q\) and updating \[\begin{split}U & = & U - XQ \\\end{split}\] If \(B\ne I\), the :math:bar{V} must also be updated as \[\begin{split}\bar{V} & = & \bar{V} - BXQ,\end{split}\] or \[\bar{V} = BU\]
22	\(B\)-orthogonalize columns of \(U\) to all vectors \(X\) by solving \[(X^TBX) Q = X^T U\] for \(Q\) and updating \[U = U - BXQ\]
999	Restart: If rci%k>0: Restart suggested with block size m >= rci%nx + rci%i + rci%j, adjusting workspace size to match. Set rci%i=0 and rci%j=0 and recall the routine. If a restart is not desirable, routine may be recalled with no change to parameters. If rci%k=0: Restart required with the same block size. In both cases, the first block W(:,:,0) should retain vectors rci%jx:rci%jx+rci%nx-1, filling remaining vectors randomly such that the entire set of columns is linearly independent from each other and also from the converged eigenvectors.

The matrices are defined as follows:

• \(U\) = W(:, rci%jx:rci%jx+rci%nx-1, rci%kx)
• \(V\) = W(:, rci%jy:rci%jy+rci%ny-1, rci%ky)
• \(\bar{V}\) = W(:, rci%jy:rci%jy+rci%nx-1, rci%ky)
• \(R\) = rr(rci%i:rci%i+rci%nx-1, rci%j:rci%j+rci%ny-1, rci%k)

and \(\alpha\) and \(\beta\) are given by rci%alpha and rci%beta respectively. We use the notation \(r_{ii}\) to refer to the \(i\)-th diagonal element of \(R\), being rr(rci%i+i-1,rci%j+i-1,rci%k).

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). problem [integer,in] :: Problem to be solved, one of: 0 \(Ax=\lambda x\) >0 \(Ax=\lambda Bx\) <0 \(ABx=\lambda x\) left [integer,in] :: Number of left eigenpairs to find. On return with rci%job=5, can be set to zero if sufficient have been found. right [integer,in] :: Number of right eigenpairs to find. On return with rci%job=5, can be set to zero if sufficient have been found. lambda (m) [real,inout] :: Current eigenvalue estimates in ascending order. m [integer,in] :: Block size of workspace W. Must be at least 2. rr (2*m,2*m,3) [real,inout] :: reverse communication workspace. (Type complex in complex version). ind (m) [integer,inout] :: reverse communication workspace. keep [ssmfe_core_keep,inout] :: Internal workspace used by routine. options [ssmfe_core_options,in] :: specifies algorithm options to be used. inform [ssmfe_inform,inout] :: returns information about the exection of the routine.

subroutine ssmfe_largest(rci, problem, nep, m, lambda, rr, ind, keep, options, inform)¶

Computes specified number of eigenvalues of largest magnitude and corresponding eigenvectors.

Uses reverse-communication. 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 as for ssmfe() above.

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). problem [integer,in] :: Problem to be solved, one of: 0 \(Ax=\lambda x\) >0 \(Ax=\lambda Bx\) <0 \(ABx=\lambda x\) nep [integer,in] :: Number of eigenpairs to find. lambda (m) [real,inout] :: Current eigenvalue estimates in ascending order. m [integer,in] :: Block size of workspace W. Must be at least 2. rr (2*m,2*m,3) [real,inout] :: reverse communication workspace. (Type complex in complex version). ind (m) [integer,inout] :: reverse communication workspace. keep [ssmfe_core_keep,inout] :: Internal workspace used by routine. options [ssmfe_core_options,in] :: specifies algorithm options to be used. inform [ssmfe_inform,inout] :: returns information about the exection of the routine.

subroutine ssmfe_free([keep, inform])¶

Free memory allocated in keep and inform. Unnecessary if both are going out of scope.

Options:	keep [ssmfe_expert_keep,inout] :: Workspace to be freed. inform [ssmfe_inform,inout] :: Information type to be freed.

Derived types¶

type ssmfe_core_options¶

Options that control the algorithm.

Type fields:

% err_est [integer,default=2] :: error estimation scheme, one of: 1 Residual error bounds: modified Davis-Kahan estimate for eigenvector error and Lehmann bounds for eigenvale error (see method section). 2 (default) Convergence curve-based estimate. % extra_left [integer,default=0] :: number of extra approximate eigenvectors corresponding to leftmost eigenvalues used to enhance convergence. % extra_right [integer,default=0] :: number of extra approximate eigenvectors corresponding to rightmost eigenvalues used to enhance convergence. % minAprod [logical,default=.true.] :: If true, minimize number of multiplications with \(A\) by requiring 2 additional blocks of memory for the workspace W(:,:,:). If false, three returns with rci%job=1 occur per iteration instead of one. Must be true if problem < 0. % minBprod [logical,default=.true.] :: If true, minimize number of multiplications with \(B\) by requiring 2 additional blocks of memory for the workspace W(:,:,:). If false, at least three returns with rci%job=3 occur per iteration instead of one. % min_gap [real,default=0.0] :: Restart sensitivity: if the relative distance between the last eigenvalue of interest on either margin of the spectrum and the rest of the spectrum is smaller than min_gap, the solver procedure suggests restart (rci%job=999). The values rci%i and rci%j are set to the numbers of eigenvalues on the left and right margin of the spectrum that are too close to the eigenvalues of interest, causing slow convergence. The default value of 0.0 means no restart is ever suggested. Must be in the range \([0.0,1.0]\). % cf_max [real,default=1.0] :: Stagnation sensitivity: if the value \(q_{ij}\) (see method section) is greater than cf_max for \(i > 5\), the eigenpair is marked as stagnated by setting inform%converged(j) to a negative value. The default value of 1.0 indicates that the estimated asymptotic convergence factor is not used for stagnation detection. Must be in the range \([0.5, 1.0]\).

type ssmfe_infrom¶

Information on progress of the algorithm.

Type fields:

% converged (mep) [integer,allocatable] :: Convergence status. If converged(j)>0, the eigenpair (lambda(j), X(j)) converged on iteration converged(j). If converged(j)=0, the eigenpair (lambda(j), X(j)) is still converging. If converged(j)<0, the eigenpair (lambda(j), X(j)) stagnated at iteration converged(j). % err_lambda (mep) [real,allocatable] :: estimated eigenvalue errors for converged and stagnated eigenvalues. % err_x (mep) [real,allocatable] :: estimated eigenvector errors for converged and stagnated eigenvectors. % flag [integer] :: return status of algorithm. See table below. % iteration [integer] :: number of iterations. % residual_norms (mep) [real,allocatable] :: Euclidean norms of residuals for (lambda(:), X(:)) on return with rci%job=4, 5. % stat [integer] :: allocation status in event of failure

inform%flag
-1	m is out-of-range.
-2	rci%job is out-of-range.
-3	options%err_est is out-of-range.
-4	options%minAprod is incompatible with selected routine.
-5	options%extra_left or options%extra_right is out-of-range.
-6	options%min_gap is out-of-range.
-7	options%cf_max is out-of-range.
-11	left is out-of-range (ssmfe()) or nep is out-of-range (ssmfe_largest()).
-12	right is out-of-range.
-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 round-off errors make them essentially random). The value of inform%non_converged is set to the number of non-converged eigenpairs.

Examples¶

! examples/Fortran/ssmfe/precond_core.f90
! Laplacian on a square grid (using SPRAL_SSMFE_CORE routines)
program ssmfe_core_precond_example
  use spral_random
  use spral_ssmfe_core
  use laplace2d ! implement Laplacian and preconditioners
  implicit none

  integer, parameter :: wp = kind(0d0) ! working precision is double

  integer, parameter :: l   = 20    ! grid points along each side
  integer, parameter :: n   = l*l   ! problem size
  integer, parameter :: nep = 5     ! eigenpairs wanted
  integer, parameter :: m = 3       ! dimension of the iterated subspace
  real(wp), parameter :: tol = 1e-6 ! eigenvector tolerance

  type(random_state) :: state       ! PRNG state

  real(wp), external :: dnrm2, ddot ! BLAS functions

  integer :: ncon                 ! number of converged eigenpairs
  integer :: i, j, k
  integer :: ind(m)               ! permutation index
  real(wp) :: s
  real(wp) :: lambda(n)           ! eigenvalues
  real(wp) :: X(n, n)             ! eigenvectors
  real(wp) :: lmd(m)              ! work array
  real(wp) :: rr(2*m, 2*m, 3)     ! work array
  real(wp) :: W(n, m, 0:5)        ! work array
  real(wp) :: U(n, m)             ! work array
  type(ssmfe_rcid        ) :: rci       ! reverse communication data
  type(ssmfe_core_options) :: options   ! options
  type(ssmfe_core_keep   ) :: keep      ! private data
  type(ssmfe_inform      ) :: inform    ! information

  ! Initialize W to lin indep vectors by randomizing
  do i=1,n
    do j=1,m
      W(i,j,0) = random_real(state, positive=.true.)
    end do
  end do

  ncon = 0
  rci%job = 0
  do ! reverse communication loop
!    call ssmfe( rci, 0, nep, 0, m, lmd, rr, ind, keep, options, inform )
!    call ssmfe( rci, 0, 0, nep, m, lmd, rr, ind, keep, options, inform )
    call ssmfe_largest( rci, 0, nep, m, lmd, rr, ind, keep, options, inform )
    select case ( rci%job )
    case ( 1 )
      call apply_laplacian &
        ( l, l, rci%nx, W(1 : n, rci%jx : rci%jx + rci%nx - 1, rci%kx), &
          W(1 : n, rci%jy : rci%jy + rci%nx - 1, rci%ky ) )
      call dscal( n*rci%nx, -1.0D0, W(1, rci%jy, rci%ky), 1 )
    case ( 2 )
      call dcopy( n*rci%nx, W(1, rci%jx, rci%kx), 1, W(1, rci%jy, rci%ky), 1 ) 
!      call apply_gauss_seidel_step &
!        ( l, l, rci%nx, W(1 : n, rci%jx : rci%jx + rci%nx - 1, rci%kx), &
!          W(1 : n, rci%jy : rci%jy + rci%nx - 1, rci%ky ) )
    case ( 999 )
      print *, rci%i, rci%j, rci%k
      exit
    case ( 4 )
      print *, 'iteration', inform%iteration
      do j = 1, m
        print '(e14.7, 2e10.1)', lmd(j), inform%residual_norms(j), inform%err_X(j)
        if ( inform%converged(j) /= 0 ) then
          cycle
        else if ( inform%err_X(j) > 0 .and. inform%err_X(j) < tol ) then
          inform%converged(j) = 1
        end if
      end do
    case ( 5 )
      if ( rci%i < 0 ) cycle
      do k = 0, rci%nx - 1
        j = ncon + k + 1
        lambda(j) = lmd(rci%jx + k)
        call dcopy( n, W(1, rci%jx + k, 0), 1, X(1, j), 1 )
      end do
      ncon = ncon + rci%nx
!      if ( ncon > 0 .or. inform%iteration > 250 ) exit
      if ( ncon >= nep .or. inform%iteration > 10 ) exit
    case ( 11 )
      if ( rci%i == 0 ) then
        if ( rci%kx /= rci%ky .or. rci%jx > rci%jy ) then
          call dcopy &
            ( n * rci%nx, W(1, rci%jx, rci%kx), 1, W(1, rci%jy, rci%ky), 1 )
        else if ( rci%jx < rci%jy ) then
          do j = rci%nx - 1, 0, -1
            call dcopy &
              ( n, W(1, rci%jx + j, rci%kx), 1, W(1, rci%jy + j, rci%ky), 1 )
          end do
        end if
      else
        do i = 1, n
          do j = 1, rci%nx
            U(i, j) = W(i, ind(j), rci%kx)
          end do
          do j = 1, rci%nx
            W(i, j, rci%kx) = U(i, j)
          end do
          if ( rci%ky /= rci%kx ) then
            do j = 1, rci%nx
              U(i, j) = W(i, ind(j), rci%ky)
            end do
            do j = 1, rci%nx
              W(i, j, rci%ky) = U(i, j)
            end do
          end if
        end do
      end if
    case ( 12 )
      do i = 0, rci%nx - 1
        rr(rci%i + i, rci%j + i, rci%k) = &
          ddot(n, W(1, rci%jx + i, rci%kx), 1, W(1, rci%jy + i, rci%ky), 1)
      end do
    case ( 13 )
      do i = 0, rci%nx - 1
        if ( rci%kx == rci%ky ) then
          s = dnrm2(n, W(1, rci%jx + i, rci%kx), 1)
          if ( s > 0 ) &
            call dscal( n, 1/s, W(1, rci%jx + i, rci%kx), 1 )
        else
          s = sqrt(abs(ddot &
            (n, W(1, rci%jx + i, rci%kx), 1, W(1, rci%jy + i, rci%ky), 1)))
          if ( s > 0 ) then
            call dscal( n, 1/s, W(1, rci%jx + i, rci%kx), 1 )
            call dscal( n, 1/s, W(1, rci%jy + i, rci%ky), 1 )
          else
            call dcopy( n, 0.0D0, 0, W(1, rci%jy + i, rci%ky), 1 )
          end if
        end if
      end do
    case ( 14 )
      do i = 0, rci%nx - 1
        s = -rr(rci%i + i, rci%j + i, rci%k)
        call daxpy&
          ( n, s, W(1, rci%jx + i, rci%kx), 1, W(1, rci%jy + i, rci%ky), 1 )
      end do    
    case ( 15 )
      if ( rci%nx > 0 .and. rci%ny > 0 ) &
        call dgemm &
          ( 'T', 'N', rci%nx, rci%ny, n, &
            rci%alpha, W(1, rci%jx, rci%kx), n, W(1, rci%jy, rci%ky), n, &
            rci%beta, rr(rci%i, rci%j, rci%k), 2*m )
    case ( 16, 17 )    
      if ( rci%ny < 1 ) cycle
      if ( rci%nx < 1 ) then
        if ( rci%job == 17 ) cycle
        if ( rci%beta == 1.0D0 ) cycle
        do j = rci%jy, rci%jy + rci%ny - 1
          call dscal( n, rci%beta, W(1, j, rci%ky), 1 )
        end do
        cycle
      end if
      if ( rci%job == 17 ) then
        call dgemm &
          ( 'N', 'N', n, rci%ny, rci%nx, &
            1.0D0, W(1, rci%jx, rci%kx), n, rr(rci%i, rci%j, rci%k), 2*m, &
            0.0D0, W(1, rci%jy, rci%ky), n )
        call dcopy &
          ( n * rci%ny, W(1, rci%jy, rci%ky), 1, W(1, rci%jx, rci%kx), 1 )
      else
        call dgemm &
          ( 'N', 'N', n, rci%ny, rci%nx, &
            rci%alpha, W(1, rci%jx, rci%kx), n, rr(rci%i, rci%j, rci%k), 2*m, &
            rci%beta, W(1, rci%jy, rci%ky), n )
      end if
    case ( 21, 22 )
      if ( ncon > 0 ) then
        call dgemm &
          ( 'T', 'N', ncon, rci%nx, n, &
            1.0D0, X, n, W(1, rci%jy, rci%ky), n, 0.0D0, U, n )
        call dgemm &
          ( 'N', 'N', n, rci%nx, ncon, &
            -1.0D0, X, n, U, n, 1.0D0, W(1, rci%jx, rci%kx), n )
      end if
    case ( :-1 )
      exit
    end select
  end do
  print '(i3, 1x, a, 1x, i3, 1x, a)', ncon, 'eigenpairs converged in', inform%iteration, 'iterations'
  print '(1x, a, i1, a, es14.7)', &
    ('lambda(', i, ') = ', lambda(i), i = 1, ncon)
  call ssmfe_free( keep, inform )
end program ssmfe_core_precond_example

 5 eigenpairs converged in  72 iterations
lambda(1) = 4.4676695E-02
lambda(2) = 1.1119274E-01
lambda(3) = 1.1119274E-01
lambda(4) = 1.7770878E-01
lambda(5) = 2.2040061E-01

Method¶