subroutines.h

Go to the documentation of this file.

#pragma once
#include "dg/topology/functions.h"
#include "dg/backend/config.h"
 
namespace dg{
struct IDENTITY
{
    template<class T>
    DG_DEVICE T operator()(T x)const{return x;}
};
 
 
 
struct equals
{
    template< class T1, class T2>
DG_DEVICE void operator()( T1 x, T2& y) const
    {
        y = x;
    }
};
struct plus_equals
{
    template< class T1, class T2>
DG_DEVICE void operator()( T1 x, T2& y) const
    {
        y += x;
    }
};
struct minus_equals
{
    template< class T1, class T2>
DG_DEVICE void operator()( T1 x, T2& y) const
    {
        y -= x;
    }
};
struct times_equals
{
    template< class T1, class T2>
DG_DEVICE void operator()( T1 x, T2& y) const
    {
        y *= x;
    }
};
struct divides_equals
{
    template< class T1, class T2>
DG_DEVICE void operator()( T1 x, T2& y) const
    {
        y /= x;
    }
};
 
 
 
struct Sum
{
    template< class T1, class ...Ts>
DG_DEVICE auto operator()( T1 x, Ts... rest) const
    {
        // unfortunately the fold expression ( x + ... + rest)
        // does currently not guarantee the order of execution
        // so we need to wait for DR 2611 to be implemented in g++ to use it
        return sum( x, rest ...);
    }
    private:
    template< class T1, class ...Ts>
DG_DEVICE auto sum( T1 x, Ts... rest) const
    {
        return x + sum( rest...);
    }
 
    template<class T1>
DG_DEVICE auto sum( T1 x1) const
    {
        return x1;
    }
};
 
 
struct Product
{
    template< class T1, class ...Ts>
DG_DEVICE auto operator()( T1 x, Ts... rest) const
    {
        // manual implement ( x * ... * rest) until DR 2611 is resolved
        return prod(x, rest...);
    }
    private:
    template< class T1, class ...Ts>
DG_DEVICE auto prod( T1 x, Ts... rest) const
    {
        return x * prod( rest...);
    }
 
    template<class T1>
DG_DEVICE auto prod( T1 x1) const
    {
        return x1;
    }
};
 
struct PairSum
{
    template< class T1, class T2, class ...Ts>
DG_DEVICE auto operator()( T1 a, T2 x, Ts... rest) const
    {
        return sum( a, x, rest...);
    }
    private:
    template<class T1, class T2, class ...Ts>
DG_DEVICE auto sum( T1 alpha, T2 x, Ts... rest) const
    {
        return DG_FMA( alpha, x, sum(rest...));
    }
 
    template< class T1, class T2>
DG_DEVICE auto sum( T1 alpha, T2 x) const
    {
        return alpha*x;
    }
};
struct TripletSum
{
    template< class T0, class T1, class T2, class ...Ts>
DG_DEVICE auto operator()( T0 a, T1 x1, T2 y1, Ts... rest) const
    {
        return sum( a, x1, y1, rest...);
    }
    private:
    template<class T0, class T1, class T2, class ...Ts>
DG_DEVICE auto sum( T0 alpha, T1 x, T2 y, Ts... rest) const
    {
        return DG_FMA( alpha*x, y, sum(rest...));
    }
 
    template<class T0, class T1, class T2>
DG_DEVICE auto sum( T0 alpha, T1 x, T2 y) const
    {
        return (alpha*x)*y;
    }
};
 
 
 
struct EmbeddedPairSum
{
    template< class T1, class ...Ts>
DG_DEVICE void operator()( T1& y, T1& yt, T1 b, T1 bt, Ts... rest) const
    {
        y = b*y;
        yt = bt*yt;
        sum( y, yt, rest...);
    }
    private:
    template< class T1,  class ...Ts>
DG_DEVICE void sum( T1& y_1, T1& yt_1, T1 b, T1 bt, T1 k, Ts... rest) const
    {
        y_1 = DG_FMA( b, k, y_1);
        yt_1 = DG_FMA( bt, k, yt_1);
        sum( y_1, yt_1, rest...);
    }
 
    template< class T1>
DG_DEVICE void sum( T1& y_1, T1& yt_1, T1 b, T1 bt, T1 k) const
    {
        y_1 = DG_FMA( b, k, y_1);
        yt_1 = DG_FMA( bt, k, yt_1);
    }
};
 
 
//The only reason the following classes exist is that nvcc does not allow
//device lambdas or local classes inside host functions
 
template<class BinarySub, class Functor>
struct Evaluate
{
    Evaluate( BinarySub sub, Functor g): m_f( sub), m_g( g) {}
#ifdef __CUDACC__
// cuda compiler spits out a lot of warnings if
// e.g. dg::transform is used on host vectors with host function
// hd_warning_disable is unfortunately undocumented, but let's try
// If it ever causes trouble we can remove it again
// it just suppresses compiler warnings:
// https://stackoverflow.com/questions/55481202/how-to-disable-cuda-host-device-warning-for-just-one-function
#pragma hd_warning_disable
#endif
    template< class T, class... Ts>
DG_DEVICE void operator() ( T& y, Ts... xs){
        m_f(m_g(xs...), y);
    }
    private:
    BinarySub m_f;
    Functor m_g;
};
 
 
template<class T>
struct Scal
{
    Scal( T a): m_a(a){}
    template<class T1>
DG_DEVICE
    void operator()( T1& y)const{
        y *= m_a;
    }
    private:
    T m_a;
};
 
template<class T>
struct Plus
{
    Plus( T a): m_a(a){}
    template<class T1>
DG_DEVICE
    void operator()( T1& y) const{
        y += m_a;
    }
    private:
    T m_a;
};
 
template<class T0, class T1>
struct Axpby
{
    Axpby( T0 a, T1 b): m_a(a), m_b(b){}
    template<class T2, class T3>
DG_DEVICE
    void operator()( T2 x, T3& y)const {
        y *= m_b;
        y = DG_FMA( m_a, x, y);
    }
    private:
    T0 m_a;
    T1 m_b;
};
template<class T0, class T1>
struct AxyPby
{
    AxyPby( T0 a, T1 b): m_a(a), m_b(b){}
    template<class T2, class T3>
DG_DEVICE
    void operator()( T2 x, T3& y)const {
        T3 tmp = y;
        y *= m_b;
        y = DG_FMA( m_a*x, tmp, y);
    }
    private:
    T0 m_a;
    T1 m_b;
};
 
template<class T0, class T1, class T2>
struct Axpbypgz
{
    Axpbypgz( T0 a, T1 b, T2 g): m_a(a), m_b(b), m_g(g){}
    template<class T3, class T4, class T5>
DG_DEVICE
    void operator()( T3 x, T4 y, T5& z)const{
        z *= m_g;
        z = DG_FMA( m_a, x, z);
        z = DG_FMA( m_b, y, z);
    }
    private:
    T0 m_a;
    T1 m_b;
    T2 m_g;
};
 
template<class T0, class T1>
struct PointwiseDot
{
    PointwiseDot( T0 a, T1 b): m_a(a), m_b(b) {}
    template<class T3, class T4, class T5>
DG_DEVICE void operator()( T3 x, T4 y, T5& z)const{
        z *= m_b;
        z = DG_FMA( m_a*x, y, z);
    }
    template<class T3, class T4, class T5, class T6>
DG_DEVICE
    void operator()( T3 x1, T4 x2, T5 x3, T6& y)const{
        y *= m_b;
        y = DG_FMA( m_a*x1, x2*x3, y);
    }
    private:
    T0 m_a;
    T1 m_b;
};
template<class T0, class T1, class T2>
struct PointwiseDot2
{
    PointwiseDot2( T0 a, T1 b, T2 g): m_a(a), m_b(b), m_g(g) {}
    template<class T3, class T4, class T5, class T6, class T7>
DG_DEVICE
    void operator()( T3 x1, T4 y1, T5 x2, T6 y2, T7& z)const{
        z *= m_g;
        z = DG_FMA( m_a*x1, y1, z);
        z = DG_FMA( m_b*x2, y2, z);
    }
    private:
    T0 m_a;
    T1 m_b;
    T2 m_g;
};
 
struct divides
{
    template< class T1, class T2>
DG_DEVICE auto operator()( T1 x1, T2 x2) const
    {
        return x1/x2;
    }
};
 
template<class T0, class T1>
struct PointwiseDivide
{
    PointwiseDivide( T0 a, T1 b): m_a(a), m_b(b){}
    template<class T3, class T4>
DG_DEVICE
    void operator()( T3 y, T4& z)const{
        T4 tmp = z;
        z *= m_b;
        z = DG_FMA( m_a, tmp/y, z);
    }
    template<class T3, class T4, class T5>
DG_DEVICE
    void operator()( T3 x, T4 y, T5& z)const{
        z *= m_b;
        z = DG_FMA( m_a, x/y, z);
    }
    private:
    T0 m_a;
    T1 m_b;
};
namespace detail
{
template<class F, class G>
struct Compose
{
    Compose( F f, G g):m_f(f), m_g(g){}
    template<class ...Xs>
    auto operator() ( Xs&& ... xs){
        return m_f(m_g(std::forward<Xs>(xs)...));
    }
    template<class ...Xs>
    auto operator() ( Xs&& ... xs) const {
        return m_f(m_g(std::forward<Xs>(xs)...));
    }
    private:
    F m_f;
    G m_g;
};
}//namespace detail
 
 
template <class UnaryOp, class Functor>
auto compose( UnaryOp f, Functor g) {
    return detail::Compose<UnaryOp,Functor>( f, g);
    //a C++-14 way of generating a generic lambda with a parameter pack. Taken from:
    //https://stackoverflow.com/questions/19071268/function-composition-in-c-c11
    //return [f,g](auto&&... xs){ return f(g(std::forward<decltype(xs)>(xs)...));};
}
template <class UnaryOp, typename... Functors>
auto compose(UnaryOp f0, Functors... fs) {
    return compose( f0 , compose(fs...));
}
 
}//namespace dg