Extension: Matrix functions
#include "dg/matrix/matrix.h"
dg::mat::UniversalLanczos< ContainerType > Class Template Reference

Tridiagonalize \(A\) and approximate \(f(A)b \approx |b|_W V f(T) e_1\) via Lanczos algorithm. A is self-adjoint in the weights \( W\). More...

Public Types

using value_type = get_value_type< ContainerType >
 value type of the ContainerType class More...
 
using HCooMatrix = cusp::coo_matrix< int, value_type, cusp::host_memory >
 
using HDiaMatrix = cusp::dia_matrix< int, value_type, cusp::host_memory >
 
using HVec = dg::HVec
 

Public Member Functions

 UniversalLanczos ()
 Allocate nothing, Call construct method before usage. More...
 
 UniversalLanczos (const ContainerType &copyable, unsigned max_iterations)
 Allocate memory for the method. More...
 
template<class ... Params>
void construct (Params &&...ps)
 Perfect forward parameters to one of the constructors. More...
 
void set_max (unsigned new_max)
 Set the maximum number of iterations. More...
 
unsigned get_max () const
 Get the current maximum number of iterations. More...
 
void set_verbose (bool verbose)
 Set or unset debugging output during iterations. More...
 
value_type get_bnorm () const
 Norm of b from last call to operator() More...
 
template<class MatrixType , class ContainerType0 , class ContainerType1 , class ContainerType2 , class FuncTe1 >
unsigned solve (ContainerType0 &x, FuncTe1 f, MatrixType &&A, const ContainerType1 &b, const ContainerType2 &weights, value_type eps, value_type nrmb_correction=1., std::string error_norm="universal", value_type res_fac=1., unsigned q=1)
 \( x = f(A)b \approx ||b||_W V f(T) e_1 \) via Lanczos and matrix function computation. A is self-adjoint in the weights \( W\). More...
 
template<class MatrixType , class ContainerType0 , class ContainerType1 >
const HDiaMatrixtridiag (MatrixType &&A, const ContainerType0 &b, const ContainerType1 &weights, value_type eps=1e-12, value_type nrmb_correction=1., std::string error_norm="universal", value_type res_fac=1., unsigned q=1)
 Tridiagonalization of A using Lanczos method. More...
 
unsigned get_iter () const
 Get the number of iterations in the last call to tridiag or solve (same as T.num_rows) More...
 

Detailed Description

template<class ContainerType>
class dg::mat::UniversalLanczos< ContainerType >

Tridiagonalize \(A\) and approximate \(f(A)b \approx |b|_W V f(T) e_1\) via Lanczos algorithm. A is self-adjoint in the weights \( W\).

The M-Lanczos method is based on the paper Novel Numerical Methods for Solving the Time-Space Fractional Diffusion Equation in Two Dimensions by Q. Yang et al, but adopts a more efficient implementation similar to that in the PCG method. Further also the conventional Lanczos method can be found there and also in text books such as Iteratvie Methods for Sparse Linear Systems 2nd edition by Yousef Saad.

Note
We have two stopping criteria. The residual criterion stops when \( \tau ||r_i||_W = \tau ||\vec b||_W \beta_i (T^{-1})_{1m} \leq \epsilon_{rel} ||\vec b||_W + \epsilon_{aps} \) where \( \tau \) is a residual factor accounting for the condition of various matrix functions
The universal stopping criterion is based on the paper Estimating the error in matrix function approximations by Q. Eshghi, N. and Reichel L., The iteration stops when

\[ ||\vec{e}_{f,m}||_W = ||\vec{b}||_W ||\left(\check{T} - f(\bar{T})\right)\vec{e}_1||_2 \leq \epsilon_{rel} ||\vec{b}||_W ||f(T)\vec e_1||_2 + \epsilon_{abs} \]

with

\[ \bar{T} = \begin{pmatrix} T & \beta_m \vec{e}_m & & \\ \beta_m \vec{e}_m^T & \alpha_{m-1} & \beta_{m-2} & \\ & \beta_{m-2} & \alpha_{m-2} & \beta_{m-1-q} \\ & & \beta_{n-1-q} & \alpha_{n-q} \end{pmatrix} \]

\[ \check{T} = \begin{pmatrix} f(T) & 0 \\ 0 & 0 \end{pmatrix} \]

The common lanczos method (and M-Lanczos) method are prone to loss of orthogonality for finite precision. Here, only the basic Paige fix is used. Thus the iterations should be kept as small as possible. Could be fixed via full, partial or selective reorthogonalization strategies, but so far no problems occured due to this.

Template Parameters
ContainerTypeAny class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others
  • dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp)
  • std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors)
  • double (scalar) and other primitive types ...
If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them
See also
See The dg dispatch system for a detailed explanation of our type dispatch system

Member Typedef Documentation

◆ HCooMatrix

template<class ContainerType >
using dg::mat::UniversalLanczos< ContainerType >::HCooMatrix = cusp::coo_matrix<int, value_type, cusp::host_memory>

◆ HDiaMatrix

template<class ContainerType >
using dg::mat::UniversalLanczos< ContainerType >::HDiaMatrix = cusp::dia_matrix<int, value_type, cusp::host_memory>

◆ HVec

template<class ContainerType >
using dg::mat::UniversalLanczos< ContainerType >::HVec = dg::HVec

◆ value_type

template<class ContainerType >
using dg::mat::UniversalLanczos< ContainerType >::value_type = get_value_type<ContainerType>

value type of the ContainerType class

Constructor & Destructor Documentation

◆ UniversalLanczos() [1/2]

template<class ContainerType >
dg::mat::UniversalLanczos< ContainerType >::UniversalLanczos ( )
inline

Allocate nothing, Call construct method before usage.

◆ UniversalLanczos() [2/2]

template<class ContainerType >
dg::mat::UniversalLanczos< ContainerType >::UniversalLanczos ( const ContainerType &  copyable,
unsigned  max_iterations 
)
inline

Allocate memory for the method.

Parameters
copyableA ContainerType must be copy-constructible from this
max_iterationsMaximum number of iterations to be used

Member Function Documentation

◆ construct()

template<class ContainerType >
template<class ... Params>
void dg::mat::UniversalLanczos< ContainerType >::construct ( Params &&...  ps)
inline

Perfect forward parameters to one of the constructors.

Template Parameters
Paramsdeduced by the compiler
Parameters
psparameters forwarded to constructors

◆ get_bnorm()

template<class ContainerType >
value_type dg::mat::UniversalLanczos< ContainerType >::get_bnorm ( ) const
inline

Norm of b from last call to operator()

Returns
bnorm

◆ get_iter()

template<class ContainerType >
unsigned dg::mat::UniversalLanczos< ContainerType >::get_iter ( ) const
inline

Get the number of iterations in the last call to tridiag or solve (same as T.num_rows)

Returns
the number of iterations in the last call to tridiag or solve

◆ get_max()

template<class ContainerType >
unsigned dg::mat::UniversalLanczos< ContainerType >::get_max ( ) const
inline

Get the current maximum number of iterations.

Returns
the current maximum

◆ set_max()

template<class ContainerType >
void dg::mat::UniversalLanczos< ContainerType >::set_max ( unsigned  new_max)
inline

Set the maximum number of iterations.

Parameters
new_maxNew maximum number

◆ set_verbose()

template<class ContainerType >
void dg::mat::UniversalLanczos< ContainerType >::set_verbose ( bool  verbose)
inline

Set or unset debugging output during iterations.

Parameters
verboseIf true, additional output will be written to std::cout during solution

◆ solve()

template<class ContainerType >
template<class MatrixType , class ContainerType0 , class ContainerType1 , class ContainerType2 , class FuncTe1 >
unsigned dg::mat::UniversalLanczos< ContainerType >::solve ( ContainerType0 &  x,
FuncTe1  f,
MatrixType &&  A,
const ContainerType1 &  b,
const ContainerType2 &  weights,
value_type  eps,
value_type  nrmb_correction = 1.,
std::string  error_norm = "universal",
value_type  res_fac = 1.,
unsigned  q = 1 
)
inline

\( x = f(A)b \approx ||b||_W V f(T) e_1 \) via Lanczos and matrix function computation. A is self-adjoint in the weights \( W\).

Tridiagonalize \(A\) using Lanczos algorithm with a residual or universal stopping criterion

Parameters
xoutput vector
fthe matrix function that is called like yH = f( T) where T is the tridiagonal matrix and returns the result of \( f(T)\vec e_1\) (for example dg::mat::make_FuncEigen_Te1( dg::SQRT<double>()) )
See also
Matrix-functions
Parameters
Aself-adjoint and semi-positive definit matrix
binput vector
weightsweights in which A is self-adjoint
epsrelative accuracy of M-Lanczos method
nrmb_correctionthe absolute error in units of eps to be respected
error_normEither "residual" or "universal"
res_facfactor \( \tau\) that is multiplied to the norm of the residual. Used to account for specific matrix function and operator in the convergence criterium
qThe q-number in the "universal stopping criterion
Returns
number of iterations of M-Lanczos routine

◆ tridiag()

template<class ContainerType >
template<class MatrixType , class ContainerType0 , class ContainerType1 >
const HDiaMatrix & dg::mat::UniversalLanczos< ContainerType >::tridiag ( MatrixType &&  A,
const ContainerType0 &  b,
const ContainerType1 &  weights,
value_type  eps = 1e-12,
value_type  nrmb_correction = 1.,
std::string  error_norm = "universal",
value_type  res_fac = 1.,
unsigned  q = 1 
)
inline

Tridiagonalization of A using Lanczos method.

Tridiagonalize \(A\) using Lanczos algorithm with a residual or universal stopping criterion

Parameters
AA self-adjoint, positive definit matrix
bThe initial vector that starts orthogonalization
weightsWeights that define the scalar product in which A is self-adjoint and in which the error norm is computed.
epsrelative accuracy of residual
nrmb_correctionthe absolute error C in units of eps to be respected
error_normEither "residual" or "universal"
res_facfactor \( \tau\) that is multiplied to the norm of the residual. Used to account for specific matrix function and operator in the convergence criterium
qThe q-number in the "universal stopping criterion
Returns
returns the tridiagonal matrix T. Note that \( T = (MV)^T A V \). The number of iterations is given by T.num_rows

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