BLAS level 2 routines: Matrix-Vector

\( \alpha M \cdot x + \beta y\) and \( x^T M y \) More...

Collaboration diagram for BLAS level 2 routines: Matrix-Vector:

Namespaces
namespace	dg::blas2
	BLAS Level 2 routines.

Functions
template<class ContainerType1 , class MatrixType , class ContainerType2 >
get_value_type< MatrixType >	dg::blas2::dot (const ContainerType1 &x, const MatrixType &m, const ContainerType2 &y)
	\( x^T M y\); Binary reproducible general dot product More...
template<class MatrixType , class ContainerType >
get_value_type< MatrixType >	dg::blas2::dot (const MatrixType &m, const ContainerType &x)
	\( x^T M x\); Binary reproducible general dot product More...
template<class MatrixType , class ContainerType1 , class ContainerType2 >
void	dg::blas2::symv (get_value_type< ContainerType1 > alpha, MatrixType &&M, const ContainerType1 &x, get_value_type< ContainerType1 > beta, ContainerType2 &y)
	\( y = \alpha M x + \beta y\) More...
template<class MatrixType , class ContainerType1 , class ContainerType2 >
void	dg::blas2::symv (MatrixType &&M, const ContainerType1 &x, ContainerType2 &y)
	\( y = M x\) More...
template<class MatrixType , class ContainerType1 , class ContainerType2 >
void	dg::blas2::gemv (get_value_type< ContainerType1 > alpha, MatrixType &&M, const ContainerType1 &x, get_value_type< ContainerType1 > beta, ContainerType2 &y)
	\( y = \alpha M x + \beta y \); (alias for symv) More...
template<class MatrixType , class ContainerType1 , class ContainerType2 >
void	dg::blas2::gemv (MatrixType &&M, const ContainerType1 &x, ContainerType2 &y)
	\( y = M x\); (alias for symv) More...
template<class Stencil , class ContainerType , class ... ContainerTypes>
void	dg::blas2::parallel_for (Stencil f, unsigned N, ContainerType &&x, ContainerTypes &&... xs)
	\( f(i, x_0, x_1, ...)\ \forall i\); Customizable and generic for loop More...
template<class FunctorType , class MatrixType , class ContainerType1 , class ContainerType2 >
void	dg::blas2::stencil (FunctorType f, MatrixType &&M, const ContainerType1 &x, ContainerType2 &y)
	\( F(M, x, y)\) More...
template<class MatrixType , class AnotherMatrixType >
void	dg::blas2::transfer (const MatrixType &x, AnotherMatrixType &y)
	\( y = x\); Generic way to copy and/or convert a Matrix type to a different Matrix type More...

Detailed Description

\( \alpha M \cdot x + \beta y\) and \( x^T M y \)

Note

Successive calls to blas routines are executed sequentially

A manual synchronization of threads or devices is never needed in an application using these functions. All functions returning a value block until the value is ready.

Function Documentation

dot() [1/2]

template<class ContainerType1 , class MatrixType , class ContainerType2 >

get_value_type< MatrixType > dg::blas2::dot ( const ContainerType1 &  x, const MatrixType &  m, const ContainerType2 &  y  )

inline

\( x^T M y\); Binary reproducible general dot product

This routine computes the scalar product defined by the symmetric positive definite matrix M

\[ x^T M y = \sum_{i,j=0}^{N-1} x_i M_{ij} y_j \]

This code snippet demonstrates how to discretize and compute the norm of a function on a shared memory system

// define the function to integrate
double function(double x, double y, double amp){
    return amp*exp(x)*exp(y);
}
 
// create a grid of the domain [0,2]x[0,2] with 20 cells in x and y and 3 polynomial coefficients
dg::Grid2d g2d( 0, 2, 0, 2, 3, 20, 20);
 
// create the Gaussian weights (volume form) for the integration
const dg::HVec w2d = dg::create::weights( g2d);
 
// our function needs to depend only on x and y so let's fix the amplitude to 2
using namespace std::placeholders; //for _1, _2
auto functor = std::bind( function, _1, _2, 2.);
 
// discretize the function on the grid
const dg::HVec vec = dg::evaluate( functor, g2d);
 
// now compute the scalar product (the L2 norm)
double norm = dg::blas2::dot( vec, w2d, vec);
 
// norm is now: (e^4-1)^2

Attention

if one of the input vectors contains Inf or NaN or the product of the input numbers reaches Inf or Nan then the behaviour is undefined and the function may throw. See dg::ISNFINITE and dg::ISNSANE in that case

Note

Our implementation guarantees binary reproducible results up to and excluding the last mantissa bit of the result. Furthermore, the sum is computed with infinite precision and the result is then rounded to the nearest double precision number. Although the products are not computed with infinite precision, the order of multiplication is guaranteed. This is possible with the help of an adapted version of the dg::exblas library and works for single and double precision.

Parameters

x	Left input
m	The diagonal Matrix.
y	Right input (may alias x)

Returns

Generalized scalar product. If x and y are vectors of containers and m is not, then we sum the results of dg::blas2::dot( x[i], m, y[i])

Note

This routine is always executed synchronously due to the implicit memcpy of the result. With mpi the result is broadcasted to all processes. Also note that the behaviour is undefined when one of the containers contains nan

Template Parameters

MatrixType

MatrixType has to have a category derived from AnyVectorTag and must be compatible with the ContainerTypes

ContainerType

Any 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

dot() [2/2]

template<class MatrixType , class ContainerType >

get_value_type< MatrixType > dg::blas2::dot ( const MatrixType &  m, const ContainerType &  x  )

inline

\( x^T M x\); Binary reproducible general dot product

Equivalent to dg::blas2::dot( x,m,x)

\[ x^T M x = \sum_{i,j=0}^{N-1} x_i M_{ij} x_j \]

Parameters

m	The diagonal Matrix
x	Right input

Returns

Generalized scalar product. If x is a vector of containers and m is not, then we sum the results of dg::blas2::dot( m, x[i])

Note

This routine is always executed synchronously due to the implicit memcpy of the result.

This routine is equivalent to the call dg::blas2::dot( x, m, x); which should be prefered because it looks more explicit

Template Parameters

MatrixType

MatrixType has to have a category derived from AnyVectorTag and must be compatible with the ContainerTypes

ContainerType

Any 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

gemv() [1/2]

template<class MatrixType , class ContainerType1 , class ContainerType2 >

void dg::blas2::gemv ( get_value_type< ContainerType1 >  alpha, MatrixType &&  M, const ContainerType1 &  x, get_value_type< ContainerType1 >  beta, ContainerType2 &  y  )

inline

\( y = \alpha M x + \beta y \); (alias for symv)

This routine computes

\[ y = \alpha M x + \beta y \]

where \( M\) is a matrix (or functor that is called like M(alpha,x,beta,y)). There is nothing that prevents you from making the matrix M non-symmetric or even non-linear. In this sense the term "symv" (symmetrix-Matrix-Vector multiplication) is misleading. For better code readability we introduce aliases: dg::blas2::gemv (general Matrix-Vector multiplication) and dg::apply (general, possibly non-linear functor application).

This code snippet demonstrates how to derive a function on a device

//define a function to derive
double function(double x, double y){
    return sin(x)*sin(y);
}
 
//create a grid of the domain [0,2]x[0,2] with 20 cells in x and y and 3 polynomial coefficients
dg::Grid2d g2d( 0, 2, 0, 2, 3, 20, 20, dg::DIR);
 
//declare a device matrix
dg::DMatrix dx;
 
//create the x derivative on the grid and transfer the matrix to the device
dg::blas2::transfer( dg::create::dx(g2d), dx);
 
//discretize the function on the grid and transfer the result to the device
const dg::DVec x = dg::construct<dg::DVec>( dg::evaluate( function, g2d));
 
//allocate memory for the result
dg::DVec y=x;
 
//apply the derivative to x and store result in y
dg::blas2::symv(dx, x,y);
 
//or equivalently
dg::blas2::symv(1., dx, x, 0., y);

Parameters

alpha	A Scalar
M	The Matrix.
x	input vector
beta	A Scalar
y	contains the solution on output (may not alias x)

Attention

y may not alias x, the only exception is if MatrixType has the AnyVectorTag

If y on input contains a NaN or Inf, it may contain NaN or Inf on output as well even if beta is zero.

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example scalar or Container: Scalars and containers act as diagonal matrices. std::function<void( const ContainerType0&, ContainerType1&)> dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec> dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec> dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec> dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec> In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

ContainerType

Any 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

gemv() [2/2]

template<class MatrixType , class ContainerType1 , class ContainerType2 >

void dg::blas2::gemv ( MatrixType &&  M, const ContainerType1 &  x, ContainerType2 &  y  )

inline

\( y = M x\); (alias for symv)

This routine computes

\[ y = M x \]

where \( M\) is a matrix (or functor that is called like M(x,y)). There is nothing that prevents you from making the matrix M non-symmetric or even non-linear. In this sense the term "symv" (symmetrix-Matrix-Vector multiplication) is misleading. For better code readability we introduce aliases: dg::blas2::gemv (general Matrix-Vector multiplication) and dg::apply (general, possibly non-linear functor application)

This code snippet demonstrates how to derive a function on a device

//define a function to derive
double function(double x, double y){
    return sin(x)*sin(y);
}
 
//create a grid of the domain [0,2]x[0,2] with 20 cells in x and y and 3 polynomial coefficients
dg::Grid2d g2d( 0, 2, 0, 2, 3, 20, 20, dg::DIR);
 
//declare a device matrix
dg::DMatrix dx;
 
//create the x derivative on the grid and transfer the matrix to the device
dg::blas2::transfer( dg::create::dx(g2d), dx);
 
//discretize the function on the grid and transfer the result to the device
const dg::DVec x = dg::construct<dg::DVec>( dg::evaluate( function, g2d));
 
//allocate memory for the result
dg::DVec y=x;
 
//apply the derivative to x and store result in y
dg::blas2::symv(dx, x,y);
 
//or equivalently
dg::blas2::symv(1., dx, x, 0., y);

Parameters

M	The Matrix.
x	input vector
y	contains the solution on output (may not alias x)

Attention

y may not alias x, the only exception is if MatrixType has the AnyVectorTag and ContainerType1 ==ContainerType2

If y on input contains a NaN or Inf it is not a prioriy clear if it will contain a NaN or Inf on output as well. Our own matrix formats overwrite y and work correctly but for third-party libraries it is worth double-checking.

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example scalar or Container: Scalars and containers act as diagonal matrices. std::function<void( const ContainerType0&, ContainerType1&)> dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec> dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec> dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec> dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec> In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

ContainerType

Any 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

parallel_for()

template<class Stencil , class ContainerType , class ... ContainerTypes>

void dg::blas2::parallel_for ( Stencil  f, unsigned  N, ContainerType &&  x, ContainerTypes &&...  xs  )

inline

\( f(i, x_0, x_1, ...)\ \forall i\); Customizable and generic for loop

Attention

Only works for shared memory vectors (or scalars): no MPI, no Recursive (find reasons below).

For trivially parallel operations (no neighboring points involved) use dg::blas1::subroutine

This routine loops over an arbitrary user-defined "loop body" functor f with an arbitrary number of arguments \( x_s\) elementwise

\[ f(i, x_{0}, x_{1}, ...) \]

where i iterates from 0 to a given size N. The order of iterations is undefined. It is equivalent to the following

for(unsigned i=0; i<N; i++)
    f( i, *x_0[0], *x_1[0], ...);

With this function very general for-loops can be parallelized like for example a forward finite difference:

dg::DVec x( 100,2), y(100,4);
unsigned N = 100;
double hx = 1.;
// implement forward difference with periodic boundaries
dg::blas1::parallel_for( [&]DG_DEVICE( unsigned i, const double* x, double* y){
    unsigned ip = (i+1)%N;
    y[i] = (x[ip] - x[i])/hx;
}, N, x, y);
// y[i] now has the value 0

Note

In a way this function is a generalization of dg::blas1::subroutine to non-trivial parallelization tasks. However, this comes at a price: this function only works for containers with the dg::SharedVectorTag and sclar types. The reason it cannot work for MPI is that the stencil (and thus the communication pattern) is unkown. However, it can serve as an important building block for other parallel functions like dg::blas2::parallel_for.

This is the closest function we have to kokkos::parallel_for of the Kokkos library.

Parameters

f	the loop body
N	the total number of iterations in the for loop
x	the first argument
xs	other arguments

Attention

The user has to decide whether or not it is safe to alias input or output vectors. If in doubt, do not alias output vectors.

Template Parameters

Stencil

a function or functor with an arbitrary number of arguments and no return type; The first argument is an unsigned (the loop iteration), afterwards takes a const_pointer_type argument (const pointer to first element in vector) for each input argument in the call and a pointer_type argument (pointer to first element in vector) for each output argument. Scalars are forwarded "as is" scalar_type . Stencil must be callable on the device in use. In particular, with CUDA it must be a functor (not a function) and its signature must contain the __device__ specifier. (s.a. DG_DEVICE)

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the SharedVectorTag and AnyPolicyTag. Among others dg::HVec (serial), dg::DVec (cuda / omp) int, double and other primitive types ...

stencil()

template<class FunctorType , class MatrixType , class ContainerType1 , class ContainerType2 >

void dg::blas2::stencil ( FunctorType  f, MatrixType &&  M, const ContainerType1 &  x, ContainerType2 &  y  )

inline

\( F(M, x, y)\)

This routine calls

\[ F(i, [M], x, y) \]

for all \( i \in [0,N[\), where N is the number of rows in M, using dg::blas2::parallel_for, where [M] depends on the matrix type:

• for a csr matrix it is [M] = m.row_offsets, m.column_indices, m.values

Possible shared memory implementation

dg::blas2::parallel_for( F, m.num_rows, m.row_offsets, m.column_indices, m.values, x, y);

Other matrix types have not yet been implemented.

Note

Since the matrix is known, a communication pattern is available and thus the function works in parallel for MPI (unlike dg::blas2::parallel_for).

In a way this function is a generalization of dg::blas2::parallel_for to MPI vectors at the cost of having to encode the communication stencil in the matrix M and only one vector argument

Parameters

f	The filter function is called like f(i, m.row_offsets_ptr, m.column_indices_ptr, m.values_ptr, x_ptr, y_ptr)
M	The Matrix.
x	input vector
y	contains the solution on output (may not alias x)

Template Parameters

FunctorType	A type that is callable void operator()( unsigned, pointer, [m_pointers], const_pointer) For GPU vector the functor must be callable on the device.
MatrixType	So far only one of the cusp::csr_matrix types and their MPI variants dg::MPIDistMat<cusp::csr_matrix, Comm> are allowed

See also

dg::convert, dg::CSRMedianFilter, dg::create::window_stencil

Template Parameters

ContainerType

Any 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

symv() [1/2]

template<class MatrixType , class ContainerType1 , class ContainerType2 >

void dg::blas2::symv ( get_value_type< ContainerType1 >  alpha, MatrixType &&  M, const ContainerType1 &  x, get_value_type< ContainerType1 >  beta, ContainerType2 &  y  )

inline

\( y = \alpha M x + \beta y\)

This routine computes

\[ y = \alpha M x + \beta y \]

where \( M\) is a matrix (or functor that is called like M(alpha,x,beta,y)). There is nothing that prevents you from making the matrix M non-symmetric or even non-linear. In this sense the term "symv" (symmetrix-Matrix-Vector multiplication) is misleading. For better code readability we introduce aliases: dg::blas2::gemv (general Matrix-Vector multiplication) and dg::apply (general, possibly non-linear functor application).

This code snippet demonstrates how to derive a function on a device

//define a function to derive
double function(double x, double y){
    return sin(x)*sin(y);
}
 
//create a grid of the domain [0,2]x[0,2] with 20 cells in x and y and 3 polynomial coefficients
dg::Grid2d g2d( 0, 2, 0, 2, 3, 20, 20, dg::DIR);
 
//declare a device matrix
dg::DMatrix dx;
 
//create the x derivative on the grid and transfer the matrix to the device
dg::blas2::transfer( dg::create::dx(g2d), dx);
 
//discretize the function on the grid and transfer the result to the device
const dg::DVec x = dg::construct<dg::DVec>( dg::evaluate( function, g2d));
 
//allocate memory for the result
dg::DVec y=x;
 
//apply the derivative to x and store result in y
dg::blas2::symv(dx, x,y);
 
//or equivalently
dg::blas2::symv(1., dx, x, 0., y);

Parameters

alpha	A Scalar
M	The Matrix.
x	input vector
beta	A Scalar
y	contains the solution on output (may not alias x)

Attention

y may not alias x, the only exception is if MatrixType has the AnyVectorTag

If y on input contains a NaN or Inf, it may contain NaN or Inf on output as well even if beta is zero.

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example scalar or Container: Scalars and containers act as diagonal matrices. std::function<void( const ContainerType0&, ContainerType1&)> dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec> dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec> dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec> dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec> In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

ContainerType

Any 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

symv() [2/2]

template<class MatrixType , class ContainerType1 , class ContainerType2 >

void dg::blas2::symv ( MatrixType &&  M, const ContainerType1 &  x, ContainerType2 &  y  )

inline

\( y = M x\)

This routine computes

\[ y = M x \]

where \( M\) is a matrix (or functor that is called like M(x,y)). There is nothing that prevents you from making the matrix M non-symmetric or even non-linear. In this sense the term "symv" (symmetrix-Matrix-Vector multiplication) is misleading. For better code readability we introduce aliases: dg::blas2::gemv (general Matrix-Vector multiplication) and dg::apply (general, possibly non-linear functor application)

This code snippet demonstrates how to derive a function on a device

//define a function to derive
double function(double x, double y){
    return sin(x)*sin(y);
}
 
//create a grid of the domain [0,2]x[0,2] with 20 cells in x and y and 3 polynomial coefficients
dg::Grid2d g2d( 0, 2, 0, 2, 3, 20, 20, dg::DIR);
 
//declare a device matrix
dg::DMatrix dx;
 
//create the x derivative on the grid and transfer the matrix to the device
dg::blas2::transfer( dg::create::dx(g2d), dx);
 
//discretize the function on the grid and transfer the result to the device
const dg::DVec x = dg::construct<dg::DVec>( dg::evaluate( function, g2d));
 
//allocate memory for the result
dg::DVec y=x;
 
//apply the derivative to x and store result in y
dg::blas2::symv(dx, x,y);
 
//or equivalently
dg::blas2::symv(1., dx, x, 0., y);

Parameters

M	The Matrix.
x	input vector
y	contains the solution on output (may not alias x)

Attention

y may not alias x, the only exception is if MatrixType has the AnyVectorTag and ContainerType1 ==ContainerType2

If y on input contains a NaN or Inf it is not a prioriy clear if it will contain a NaN or Inf on output as well. Our own matrix formats overwrite y and work correctly but for third-party libraries it is worth double-checking.

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example scalar or Container: Scalars and containers act as diagonal matrices. std::function<void( const ContainerType0&, ContainerType1&)> dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec> dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec> dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec> dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec> In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

ContainerType

Any 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

transfer()

template<class MatrixType , class AnotherMatrixType >

void dg::blas2::transfer ( const MatrixType &  x, AnotherMatrixType &  y  )

inline

\( y = x\); Generic way to copy and/or convert a Matrix type to a different Matrix type

e.g. from CPU to GPU, or double to float, etc.

Template Parameters

MatrixType

Any class for which a specialization of TensorTraits exists and defines a tensor_category derived from AnyMatrixTag. Furthermore, any functor/lambda type with signature void operator()( const ContainerType0&, ContainerType1&) . For example scalar or Container: Scalars and containers act as diagonal matrices. std::function<void( const ContainerType0&, ContainerType1&)> dg::HMatrix and dg::IHMatrix with dg::HVec or std::vector<dg::HVec> dg::DMatrix and dg::IDMatrix with dg::DVec or std::vector<dg::DVec> dg::MHMatrix with dg::MHVec or std::vector<dg::MHVec> dg::MDMatrix with dg::MDVec or std::vector<dg::MDVec> In case of SelfMadeMatrixTag only those blas2 functions that have a corresponding member function in the Matrix class (e.g. symv( const ContainerType0&, ContainerType1&); ) can be called. If a Container has the RecursiveVectorTag, then the matrix is applied to each of the elements unless the type has the SelfMadeMatrixTag or is a Functor type.

AnotherMatrix

Another Matrix type

Parameters

x	source
y	sink

Note

y gets resized properly

This code snippet demonstrates how to derive a function on a device

//define a function to derive
double function(double x, double y){
    return sin(x)*sin(y);
}
 
//create a grid of the domain [0,2]x[0,2] with 20 cells in x and y and 3 polynomial coefficients
dg::Grid2d g2d( 0, 2, 0, 2, 3, 20, 20, dg::DIR);
 
//declare a device matrix
dg::DMatrix dx;
 
//create the x derivative on the grid and transfer the matrix to the device
dg::blas2::transfer( dg::create::dx(g2d), dx);
 
//discretize the function on the grid and transfer the result to the device
const dg::DVec x = dg::construct<dg::DVec>( dg::evaluate( function, g2d));
 
//allocate memory for the result
dg::DVec y=x;
 
//apply the derivative to x and store result in y
dg::blas2::symv(dx, x,y);
 
//or equivalently
dg::blas2::symv(1., dx, x, 0., y);