Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.1.1
\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
Functions
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.

Function Documentation

◆ givens_rotation()

template<typename NumberType >
std::array<NumberType, 3> 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.
Author
Denis Davydov, 2017

◆ hyperbolic_rotation()

template<typename NumberType >
std::array<NumberType, 3> 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.
Author
Denis Davydov, 2017

◆ 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().
Author
Denis Davydov, 2017

◆ 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.

chebyshev8.png

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.
Author
Denis Davydov, 2017