 Reference documentation for deal.II version Git c1d6723b63 2020-01-25 23:51:38 -0500
Utilities::LinearAlgebra Namespace Reference

## Functions

template<typename NumberType >
std::array< NumberType, 3 > givens_rotation (const NumberType &x, const NumberType &y)

template<typename NumberType >
std::array< NumberType, 3 > hyperbolic_rotation (const NumberType &x, const NumberType &y)

template<typename OperatorType , typename VectorType >
double lanczos_largest_eigenvalue (const OperatorType &H, const VectorType &v0, const unsigned int k, VectorMemory< VectorType > &vector_memory, std::vector< double > *eigenvalues=nullptr)

template<typename OperatorType , typename VectorType >
void chebyshev_filter (VectorType &x, const OperatorType &H, const unsigned int n, const std::pair< double, double > unwanted_spectrum, const double tau, VectorMemory< VectorType > &vector_memory)

## Detailed Description

A collection of linear-algebra utilities.

## ◆ givens_rotation()

template<typename NumberType >
 std::array Utilities::LinearAlgebra::givens_rotation ( const NumberType & x, const NumberType & y )

Return the elements of a continuous Givens rotation matrix and the norm of the input vector.

That is for a given pair x and y, return $$c$$ , $$s$$ and $$\sqrt{x^2+y^2}$$ such that

$\begin{bmatrix} c & s \\ -s & c \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} \sqrt{x^2+y^2} \\ 0 \end{bmatrix}$

Note
The function is implemented for real valued numbers only.

## ◆ hyperbolic_rotation()

template<typename NumberType >
 std::array Utilities::LinearAlgebra::hyperbolic_rotation ( const NumberType & x, const NumberType & y )

Return the elements of a hyperbolic rotation matrix.

That is for a given pair x and y, return $$c$$ , $$s$$ and $$r$$ such that

$\begin{bmatrix} c & -s \\ -s & c \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} = \begin{bmatrix} r \\ 0 \end{bmatrix}$

Real valued solution only exists if $$|x|>|g|$$, the function will throw an error otherwise.

Note
The function is implemented for real valued numbers only.

## ◆ lanczos_largest_eigenvalue()

template<typename OperatorType , typename VectorType >
 double Utilities::LinearAlgebra::lanczos_largest_eigenvalue ( const OperatorType & H, const VectorType & v0, const unsigned int k, VectorMemory< VectorType > & vector_memory, std::vector< double > * eigenvalues = nullptr )

Estimate an upper bound for the largest eigenvalue of H by a k -step Lanczos process starting from the initial vector v0. Typical values of k are below 10. This estimator computes a k-step Lanczos decomposition $$H V_k=V_k T_k+f_k e_k^T$$ where $$V_k$$ contains k Lanczos basis, $$V_k^TV_k=I_k$$, $$T_k$$ is the tridiagonal Lanczos matrix, $$f_k$$ is a residual vector $$f_k^TV_k=0$$, and $$e_k$$ is the k-th canonical basis of $$R^k$$. The returned value is $$||T_k||_2 + ||f_k||_2$$. If eigenvalues is not nullptr, the eigenvalues of $$T_k$$ will be written there.

vector_memory is used to allocate memory for temporary vectors. OperatorType has to provide vmult operation with VectorType.

This function implements the algorithm from

@article{Zhou2006,
Title = {Self-consistent-field Calculations Using Chebyshev-filtered
Subspace Iteration},
Author = {Zhou, Yunkai and Saad, Yousef and Tiago, Murilo L. and
Chelikowsky, James R.},
Journal = {Journal of Computational Physics},
Year = {2006},
Volume = {219},
Pages = {172--184},
}
Note
This function uses Lapack routines to compute the largest eigenvalue of $$T_k$$.
This function provides an alternate estimate to that obtained from several steps of SolverCG with SolverCG<VectorType>::connect_eigenvalues_slot().

## ◆ chebyshev_filter()

template<typename OperatorType , typename VectorType >
 void Utilities::LinearAlgebra::chebyshev_filter ( VectorType & x, const OperatorType & H, const unsigned int n, const std::pair< double, double > unwanted_spectrum, const double tau, VectorMemory< VectorType > & vector_memory )

Apply Chebyshev polynomial of the operator H to x. For a non-defective operator $$H$$ with a complete set of eigenpairs $$H \psi_i = \lambda_i \psi_i$$, the action of a polynomial filter $$p$$ is given by $$p(H)x =\sum_i a_i p(\lambda_i) \psi_i$$, where $$x=: \sum_i a_i \psi_i$$. Thus by appropriately choosing the polynomial filter, one can alter the eigenmodes contained in $$x$$.

This function uses Chebyshev polynomials of first kind. Below is an example of polynomial $$T_n(x)$$ of degree $$n=8$$ normalized to unity at $$-1.2$$. By introducing a linear mapping $$L$$ from unwanted_spectrum to $$[-1,1]$$, we can dump the corresponding modes in x. The higher the polynomial degree $$n$$, the more rapid it grows outside of the $$[-1,1]$$. In order to avoid numerical overflow, we normalize polynomial filter to unity at tau. Thus, the filtered operator is $$p(H) = T_n(L(H))/T_n(L(\tau))$$.

The action of the Chebyshev filter only requires evaluation of vmult() of H and is based on the recursion equation for Chebyshev polynomial of degree $$n$$: $$T_{n}(x) = 2x T_{n-1}(x) - T_{n-2}(x)$$ with $$T_0(x)=1$$ and $$T_1(x)=x$$.

vector_memory is used to allocate memory for temporary objects.

This function implements the algorithm (with a minor fix of sign of $$\sigma_1$$) from

@article{Zhou2014,
Title = {Chebyshev-filtered subspace iteration method free of sparse
diagonalization for solving the Kohn--Sham equation},
Author = {Zhou, Yunkai and Chelikowsky, James R and Saad, Yousef},
Journal = {Journal of Computational Physics},
Year = {2014},
Volume = {274},
Pages = {770--782},
}
Note
If tau is equal to std::numeric_limits<double>::infinity(), no normalization will be performed.