ezarpack/mpi/solver_symmetric.hpp - symmetric real eigenproblems

template<typename Backend>
class arpack_solver<Symmetric, Backend>

Main solver class wrapping the Implicitly Restarted Lanczos Method (IRLM) for real symmetric eigenproblems (MPI version).

This specialization of arpack_solver calls ARPACK-NG routines pdsaupd() and pdseupd() to compute approximations to a few eigenpairs of a linear operator \( \hat O \) that is real and symmetric with respect to a real positive semi-definite symmetric matrix \( \hat B \). In other words,

\[ \langle \mathbf{x},\hat O \mathbf{y} \rangle = \langle \hat O \mathbf{x}, \mathbf{y} \rangle \]
for all vectors \( \mathbf{x} \), \( \mathbf{y} \) and with the scalar product defined as
\[ \langle \mathbf{x}, \mathbf{y} \rangle = \mathbf{x}^T \hat B \mathbf{y}. \]
A variant of the Lanczos algorithm is internally used instead of the Arnoldi iteration for this class of problems.

Template Parameters:

Backend – Tag type specifying what storage backend (matrix/vector algebra library) must be used by arpack_solver. The storage backend determines types of internally stored data arrays and input/output view objects returned by methods of the class.

Backend-specific array and view types

using real_vector_t = typename storage::real_vector_type

One-dimensional data array (vector) of real numbers.

using real_matrix_t = typename storage::real_matrix_type

Two-dimensional data array (matrix) of real numbers.

using int_vector_t = typename storage::int_vector_type

One-dimensional data array (vector) of integers.

using real_vector_view_t = typename storage::real_vector_view_type

Partial view (slice) of a real vector.

using real_vector_const_view_t = typename storage::real_vector_const_view_type

Partial constant view (slice) of a real vector.

using real_matrix_const_view_t = typename storage::real_matrix_const_view_type

Partial constant view (slice) of a real matrix.

using vector_const_view_t = real_vector_const_view_t

Storage-specific view type to expose real input vectors \( \mathbf{x} \). An argument of this type is passed as input to callable objects representing linear operators \( \hat O \) and \( \hat B \).

using vector_view_t = real_vector_view_t

Storage-specific view type to expose real output vectors \( \mathbf{y} \). An argument of this type receives output from callable objects representing linear operators \( \hat O \) and \( \hat B \).

Public Types

enum Mode

Computational modes for generalized eigenproblems.

Values:

enumerator Inverse

Regular inverse mode.

Solve a generalized eigenproblem \( \hat A\mathbf{x} = \lambda \hat M\mathbf{x} \) by reduction to the canonical form with \( \hat O = \hat M^{-1} \hat A \) and \( \hat B = \hat M \), where \( \hat A \) is symmetric and \( \hat M \) is symmetric positive-definite.

enumerator ShiftAndInvert

Shift-and-Invert mode.

Solve a generalized eigenproblem \( \hat A\mathbf{x} = \lambda \hat M\mathbf{x} \) by reduction to the canonical form with \( \hat O = (\hat A - \sigma\hat M)^{-1} \hat M \) and \( \hat B = \hat M \), where \( \hat A \) is symmetric and \( \hat M \) is symmetric positive semi-definite.

enumerator Buckling

Buckling mode.

Solve a generalized eigenproblem \( \hat K\mathbf{x} = \lambda \hat K_G\mathbf{x} \) by reduction to the canonical form with \( \hat O = (\hat K - \sigma\hat K_G)^{-1} \hat K \) and \( \hat B = \hat K \), where \( \hat K \) is symmetric positive semi-definite and \( \hat K_G \) is symmetric indefinite.

enumerator Cayley

Cayley mode.

Solve a generalized eigenproblem \( \hat A\mathbf{x} = \lambda \hat M\mathbf{x} \) by reduction to the canonical form with \(\hat O = (\hat A -\sigma\hat M)^{-1}(\hat A +\sigma M)\) and \( \hat B = \hat M \), where \( \hat A \) is symmetric and \( \hat M \) is symmetric positive semi-definite.

Public Functions

inline arpack_solver(unsigned int N, MPI_Comm const &comm)

Constructs a solver object and allocates internal data buffers to be used by ARPACK-NG.

This constructor partitions N-dimensional vectors into blocks and distributes those blocks among all MPI ranks in a communicator in the most even way.

Parameters:
  • N – Dimension of the eigenproblem.

  • comm – MPI communicator.

inline arpack_solver(std::vector<unsigned int> const &block_sizes, MPI_Comm const &comm)

Constructs a solver object and allocates internal data buffers to be used by ARPACK-NG.

This constructor accepts a list of sizes of MPI rank-local vector blocks.

Parameters:
  • block_sizes – Sizes of MPI rank-local vector blocks, one element per MPI rank.

  • comm – MPI communicator.

template<typename A, typename ShiftsF = exact_shifts_f>
inline void operator()(A &&a, params_t const &params, ShiftsF shifts_f = {})

Solve a standard eigenproblem \( \hat A\mathbf{x} = \lambda\mathbf{x}\).

Parameters:
  • a – A callable object representing the linear operator \( \hat A \). It must take two arguments,

    a(vector_const_view_t in, vector_view_t out)
    
    a is expected to act on the vector view in and write the result into the vector view out, out = a*in. Both in and out are MPI rank-local blocks of their respective N-dimensional vectors. If a has matrix elements connecting blocks stored on different MPI ranks, it is user’s duty to pass data between those ranks and to appropriately update out on all ranks. Given an instance as of the arpack_solver< Symmetric, Backend > class, in is also indirectly accessible as
    as.workspace_vector(as.in_vector_n())
    
    and out is accessible as
    as.workspace_vector(as.out_vector_n())
    

  • params – Set of input parameters for the Implicitly Restarted Lanczos Method.

  • shifts_f – Functor that implements a shift selection strategy for the implicit restarts. For the expected signature of the functor, see exact_shifts_f::operator()(). When this argument is omitted, the default “Exact Shift Strategy” is used, which is the right choice in most cases.

Throws:
  • ezarpack::ncv_insufficient – No shifts could be applied during a cycle of the IRL iteration.

  • ezarpack::maxiter_reached – Maximum number of IRL iterations has been reached. All possible eigenvalues of \( \hat O \) have been found.

  • std::runtime_error – Invalid input parameters and other errors reported by ARPACK-NG routines pdsaupd() and pdseupd().

template<typename OP, typename B, typename ShiftsF = exact_shifts_f>
inline void operator()(OP &&op, B &&b, Mode mode, params_t const &params, ShiftsF shifts_f = {})

Solve a generalized eigenproblem \( \hat A\mathbf{x} = \lambda\hat M\mathbf{x}\).

Operators \( \hat O \) and \( \hat B \) mentioned below are related to \( \hat A \) and \( \hat M \) via a Mode -dependent spectral transformation.

Parameters:
  • op – A callable object representing the linear operator \( \hat O \). It must take two arguments,

    op(vector_view_t in, vector_view_t out)
    
    In all computational modes except for Inverse, op is expected to act on the vector view in and write the result into the vector view out, out = op*in. In the Inverse mode, however, op must do the following, in = op*in, out = M^{-1}*in. Both in and out are MPI rank-local blocks of their respective N-dimensional vectors. If op has matrix elements connecting blocks stored on different MPI ranks, it is user’s duty to pass data between those ranks and to appropriately update in and out on all ranks. Given an instance as of the arpack_solver< Symmetric, Backend > class, in is also indirectly accessible as
    as.workspace_vector(as.in_vector_n())
    
    and out is accessible as
    as.workspace_vector(as.out_vector_n())
    

  • b – A callable object representing the linear operator \( \hat B \). It must take two arguments,

    b(vector_const_view_t in, vector_view_t out)
    
    b is expected to act on the vector view in and write the result into the vector view out, out = b*in. Both in and out are MPI rank-local blocks of their respective N-dimensional vectors. If b has matrix elements connecting blocks stored on different MPI ranks, it is user’s duty to pass data between those ranks and to appropriately update out on all ranks.

  • modeComputational mode to be used.

  • params – Set of input parameters for the Implicitly Restarted Lanczos Method.

  • shifts_f – Functor that implements a shift selection strategy for the implicit restarts. For the expected signature of the functor, see exact_shifts_f::operator()(). When this argument is omitted, the default “Exact Shift Strategy” is used, which is the right choice in most cases.

Throws:
  • ezarpack::ncv_insufficient – No shifts could be applied during a cycle of the IRL iteration.

  • ezarpack::maxiter_reached – Maximum number of IRL iterations has been reached. All possible eigenvalues of \( \hat O \) have been found.

  • std::runtime_error – Invalid input parameters and other errors reported by ARPACK-NG routines pdsaupd() and pdseupd().

inline int dim() const

Returns dimension of the eigenproblem.

inline MPI_Comm mpi_comm() const

Returns MPI communicator used to construct this solver.

inline int local_block_start() const

Returns the index of the first vector element in the local block.

inline int local_block_size() const

Returns the size of the local block.

inline int in_vector_n() const

Returns the index of the workspace vector, which is currently expected to be acted upon by linear operator \( \hat O \) or \( \hat B \).

inline int out_vector_n() const

Returns the index of the workspace vector, which is currently expected to receive result from application of linear operator \( \hat O \) or \( \hat B \).

inline real_vector_view_t workspace_vector(int n)

Returns a view of an MPI rank-local block of a vector within ARPACK-NG’s workspace array.

Parameters:

n – Index of the workspace vector. Valid values are 0, 1 and 2.

Throws:

std::runtime_error – Invalid index value.

inline unsigned int nconv() const

Number of “converged” Ritz values.

inline real_vector_const_view_t eigenvalues() const

Returns a constant view of a list of nconv() eigenvalues.

The values in the list are in ascending order.

Note

In the generalized eigenproblem modes, this method always returns eigenvalues of the original problem.

inline real_matrix_const_view_t eigenvectors() const

Returns a constant view of a matrix, whose nconv() columns are MPI rank-local blocks of converged Ritz basis vectors (eigenvectors).

Throws:

std::runtime_error – Ritz vectors have not been computed in the last IRLM run.

inline real_vector_view_t residual_vector()

Returns a view of the MPI rank-local block of the current residual vector.

When params_t::random_residual_vector is set to false, the view returned by this accessor can be used to set the initial residual vector.

inline bool Bx_available() const

Has \( \hat B\mathbf{x} \) already been computed at the current IRLM iteration?

inline real_vector_const_view_t Bx_vector() const

Returns a constant view of the MPI rank-local block of the most recently computed vector \( \hat B\mathbf{x} \).

inline stats_t stats() const

Returns computation statistics from the last IRLM run.

struct exact_shifts_f

If this functor is used to provide shifts for implicit restart, then the default ARPACK-NG’s shift strategy (Exact Shift Strategy) will be employed.

See also

Paragraph 4.4.1 of ARPACK Users’ Guide: Solution of Large Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods (R. B. Lehoucq, D. C. Sorensen, C. Yang, SIAM, 1998), http://li.mit.edu/Archive/Activities/Archive/CourseWork/Ju_Li/MITCourses/18.335/Doc/ARPACK/Lehoucq97.pdf#page=64

Public Functions

inline void operator()(real_vector_const_view_t ritz_values, real_vector_const_view_t ritz_bounds, real_vector_view_t shifts)

Trivial call operator. The actual shifts will be internally computed by ARPACK-NG.

Parameters:
  • ritz_values[in] View of a real vector with current params_t::ncv Ritz values.

  • ritz_bounds[in] View of a real vector with current estimated error bounds of the Ritz values.

  • shifts[out] Real vector view to receive the computed shifts.

struct params_t

Input parameters of the Implicitly Restarted Lanczos Method (IRLM).

Public Types

enum eigenvalues_select_t

Categories of eigenvalues to compute.

Values:

enumerator Largest

Largest (algebraic) eigenvalues.

enumerator Smallest

Smallest (algebraic) eigenvalues.

enumerator LargestMagnitude

Largest eigenvalues in magnitude.

enumerator SmallestMagnitude

Smallest eigenvalues in magnitude.

enumerator BothEnds

Eigenvalues at both ends of the spectrum. If n_eigenvalues is odd, compute one more from the high end than from the low end.

Public Functions

inline params_t(unsigned int n_eigenvalues, eigenvalues_select_t eigenvalues_select, bool compute_eigenvectors)

Constructs an IRLM parameter object with given n_eigenvalues, eigenvalues_select and compute_eigenvectors. The rest of the parameters are set to their defaults.

Public Members

unsigned int n_eigenvalues

Number of eigenvalues (Ritz values) to compute.

eigenvalues_select_t eigenvalues_select

Which of the eigenvalues to compute?

int ncv = -1

Number of Lanczos vectors to be generated. -1 stands for the default value min(2*n_eigenvalues + 2, N), where N is the dimension of the problem.

bool compute_eigenvectors

Compute Ritz vectors in addition to the eigenvalues?

bool random_residual_vector = true

Use a randomly generated initial residual vector?

double sigma = 0

Eigenvalue shift \( \sigma \) used if a spectral transformation is employed.

double tolerance = 0

Relative tolerance for Ritz value convergence. The default setting is machine precision.

unsigned int max_iter = INT_MAX

Maximum number of IRLM iterations allowed.

struct stats_t

Statistics regarding a completed IRLM run.

Public Members

unsigned int n_iter

Number of Arnoldi update iterations taken.

unsigned int n_op_x_operations

Total number of \( \hat O \mathbf{x} \) operations.

unsigned int n_b_x_operations

Total number of \( \hat B \mathbf{x} \) operations.

unsigned int n_reorth_steps

Total number of steps of re-orthogonalization.