Class xlifepp::SolverUtils#

template<class ScalarType, class MV, class OP> class SolverUtils#

Collaboration diagram for xlifepp::SolverUtils:

xlifepp’s templated, static class providing utilities for the solvers.

This class provides concrete, templated implementations of utilities necessary for the solvers. These utilities include sorting, orthogonalization, projecting/solving local eigensystems, and sanity checking. These are internal utilties, so the user should not alter this class.

Constructor/Destructor

SolverUtils()#: Constructor.

inline virtual ~SolverUtils()#: Destructor.

Sorting Methods

static void permuteVectors(const int n, const std::vector<int> &perm, MV &Q, std::vector<MagnitudeType> *resids = nullptr)#: Permute the vectors in a multivector according to the permutation vector perm, and optionally the residual vector resids.

static void permuteVectors(const std::vector<int> &perm, MatrixDense &Q)#: Permute the columns of a Teuchos::SerialDenseMatrix according to the permutation vector perm.

Basis update methods

static void applyHouse(int k, MV &V, const MatrixDense &H, const std::vector<ScalarType> &tau, SmartPtr<MV> workMV = _smPtrNull)#

apply a sequence of Householder reflectors (from GEQRF) to a multivector, using minimal workspace.

This routine applies a sequence of Householder reflectors, $H_{1} H_{2} \dots H_{k}$ , to a multivector $V$ . The reflectors are applied individually, as rank-one updates to the multivector. The benefit of this is that the only required workspace is a one-column multivector. This workspace can be provided by the user. If it is not, it will be allocated locally on each call to applyHouse.

Each $H_{i}$ ( $i = 1, \dots, k \leq n$ ) has the form $H_{i} = I - τ_{i} v_{i} v_{i}^{T}$ where $τ_{i}$ is a scalar and $v_{i}$ is a vector with $v_{i} (1 : i - 1) = 0$ and $e_{i}^{T} v_{i} = 1$ ; $v (i + 1 : n)$ is stored below H(i,i) and $τ_{i}$ in tau[i-1]

.

If the multivector is

m \times n

and we apply

k

Householder reflectors, the total cost of the method is

4 m n k - 2 m (k^{2} - k)

flops. For

k = n

, this becomes

2 m n^{2}

, the same as for a matrix-matrix multiplication by the accumulated Householder reflectors.

Parameters:

k – [in] the number of Householder reflectors composing the product
V – [in/out] the multivector to be modified, with $n$ columns
H – [in] a $n \times k$ matrix containing the encoded Householder vectors, as returned from GEQRF (see below)
tau – [in] the $n$ coefficients for the Householder reflects, as returned from GEQRF
workMV – [work] (optional) a multivector used for workspace. it need contain only a single vector; it if contains more, only the first vector will be modified.

Note

zero-based indexing used for data structures H and tau, while one-based indexing used for mathematic object $v_{i}$ .

Eigensolver Projection Methods

static int directSolver(int size, const MatrixDense &KK, SmartPtr<const MatrixDense> MM, MatrixDense &EV, std::vector<MagnitudeType> &theta, int &nev, int esType = 0)#

Routine for computing the first NEV generalized eigenpairs of the Hermitian pencil (KK, MM)

  @param size [in] Dimension of the eigenproblem (KK, MM)
  @param KK [in] Hermitian "stiffness" matrix 
  @param MM [in] Hermitian positive-definite "mass" matrix
  @param EV [in] Dense matrix to store the nev eigenvectors 
  @param theta [in] Array to store the eigenvalues (Size = nev)
  @param nev [in/out] Number of the smallest eigenvalues requested (in) / computed (out)
  @param esType [in] Flag to select the algorithm
  <ul>
  <li> esType =  0  (default) (Cholesky factorization of MM)
                    with deflation of MM to get orthonormality of 
                    eigenvectors (\f$S^TMMS = I\f$)
  <li> esType =  1  (Cholesky factorization of MM)
                    (no check of orthonormality)
  <li> esType = 10  Simple eigenproblem on KK
                    (MM is not referenced in this case)
  </ul>

  \note The code accesses only the upper triangular part of KK and MM.
  \return Integer \c info on the status of the computation

Return the integer info on the status of the computation

info = 0 >> Success
info = - 20 >> Failure in LAPACK routine

Sanity Checking Methods

static NumTraits<ScalarType>::RealScalar errorEquality(const MV &X, const MV &MX, SmartPtr<const OP> M = _smPtrNull)#: Return the maximum coefficient of the matrix $M * X - M X$ scaled by the maximum coefficient of MX.

Note

When M is not specified, the identity is used.