dg::mat::SqrtCauchyInt< Container > Struct Template Reference

Cauchy integral \( \sqrt{A} b= A\frac{ 2 K' \sqrt{m}}{\pi N} \sum_{j=1}^{N} ( w_j^2 I + A)^{-1} c_j d_j b \) More...

Public Types
using	container_type = Container
using	value_type = dg::get_value_type<Container>

Public Member Functions
	SqrtCauchyInt ()
	SqrtCauchyInt (const Container &copyable)
	Construct Rhs operator.
template<class ... Params>
void	construct (Params &&...ps)
	Perfect forward parameters to one of the constructors.
const double &	w () const
	The \( w\) in \( ( w^2I + A)^{-1}\).
template<class MatrixType >
auto	make_denominator (MatrixType &A) const
template<class MatrixType0 , class MatrixType1 , class ContainerType0 , class ContainerType1 >
void	operator() (MatrixType0 &&A, MatrixType1 &&wAinv, const ContainerType0 &b, ContainerType1 &x, std::array< value_type, 2 > EVs, unsigned steps, int exp=+1)
	Cauchy integral \( x = \sqrt{A}b = A \frac{ 2 K' \sqrt{m}}{\pi N} \sum_{j=1}^{N} (w_j^2 I + A)^{-1} c_j d_j b \)

Detailed Description

template<class Container>
struct dg::mat::SqrtCauchyInt< Container >

Cauchy integral \( \sqrt{A} b= A\frac{ 2 K' \sqrt{m}}{\pi N} \sum_{j=1}^{N} ( w_j^2 I + A)^{-1} c_j d_j b \)

A is the matrix, b is the vector, w is a scalar m is the smallest eigenvalue of A, K' is the conjuated complete elliptic integral and \(c_j\) and \(d_j\) are the jacobi functions

Note

If we leave away the first A on the right hand side we approximate the inverse square root.

This class is based on the approach (method 3) of the paper Computing A alpha log(A), and Related Matrix Functions by Contour Integrals by N. Hale et al

Member Typedef Documentation

container_type

template<class Container >

using dg::mat::SqrtCauchyInt< Container >::container_type = Container

value_type

template<class Container >

using dg::mat::SqrtCauchyInt< Container >::value_type = dg::get_value_type<Container>

Constructor & Destructor Documentation

SqrtCauchyInt() [1/2]

template<class Container >

dg::mat::SqrtCauchyInt< Container >::SqrtCauchyInt ( )

inline

SqrtCauchyInt() [2/2]

template<class Container >

dg::mat::SqrtCauchyInt< Container >::SqrtCauchyInt ( const Container & copyable )

inline

Construct Rhs operator.

Parameters

copyable

Member Function Documentation

construct()

template<class Container >

template<class ... Params>

void dg::mat::SqrtCauchyInt< Container >::construct ( Params &&... ps )

inline

Perfect forward parameters to one of the constructors.

Template Parameters

Params

deduced by the compiler

Parameters

ps	parameters forwarded to constructors

make_denominator()

template<class Container >

template<class MatrixType >

auto dg::mat::SqrtCauchyInt< Container >::make_denominator ( MatrixType & A ) const

inline

The functor returning \( (w^2I + A )\) A is stored by reference (use if needed)

operator()()

template<class Container >

template<class MatrixType0 , class MatrixType1 , class ContainerType0 , class ContainerType1 >

void dg::mat::SqrtCauchyInt< Container >::operator() ( MatrixType0 && A, MatrixType1 && wAinv, const ContainerType0 & b, ContainerType1 & x, std::array< value_type, 2 > EVs, unsigned steps, int exp = +1 )

inline

Cauchy integral \( x = \sqrt{A}b = A \frac{ 2 K' \sqrt{m}}{\pi N} \sum_{j=1}^{N} (w_j^2 I + A)^{-1} c_j d_j b \)

Note

The Eigenvalues can be estimated from a few lanczos iterations (which is at least more reliable than doing it semi-analytically)

dg::mat::UniversalLanczos lanczos( A.weights(), 20);
auto T = lanczos.tridiag( A, A.weights(), A.weights());
auto EVs = dg::mat::compute_extreme_EV( T);

Parameters

A	A self-adjoint or non-self-adjoint Matrix
wAinv	The operator \( (w^2 I + A)^{-1}\) (construct with the help of w() and make_denominator(A), we provide an initial guess)
b	is input vector
x	contains the result
EVs	{minimum Eigenvalue of A, maximum Eigenvalue of A}
steps	Number of steps to use in the integration

Note

The Jacobi elliptic functions are related to the Mathematica functions via jacobi_cn(k,u ) = JacobiCN_(u,k^2), ... and the complete elliptic integral of the first kind via comp_ellint_1(k) = EllipticK(k^2)

Parameters

exp	If +1 then the sqrt is computed else the inverse sqrt

w()

template<class Container >

const double & dg::mat::SqrtCauchyInt< Container >::w ( ) const

inline

The \( w\) in \( ( w^2I + A)^{-1}\).

---

The documentation for this struct was generated from the following file:

• sqrt_cauchy.h