blas1.h

Go to the documentation of this file.

#pragma once
 
#include "backend/predicate.h"
#include "backend/tensor_traits.h"
#include "backend/tensor_traits_scalar.h"
#include "backend/tensor_traits_thrust.h"
#include "backend/tensor_traits_std.h"
#include "backend/blas1_dispatch_scalar.h"
#include "backend/blas1_dispatch_shared.h"
#ifdef MPI_VERSION
#include "backend/mpi_vector.h"
#include "backend/blas1_dispatch_mpi.h"
#endif
#include "backend/blas1_dispatch_vector.h"
#include "backend/blas1_dispatch_map.h"
#include "subroutines.h"
 
namespace dg{
 
namespace blas1
{
template< class ContainerType, class BinarySubroutine, class Functor, class ContainerType0, class ...ContainerTypes>
inline void evaluate( ContainerType& y, BinarySubroutine f, Functor g, const ContainerType0& x0, const ContainerTypes& ...xs);
 
 
template<class Functor, class ContainerType, class ...ContainerTypes>
auto vdot( Functor f, const ContainerType& x, const ContainerTypes& ...xs) ->
    std::invoke_result_t<Functor, dg::get_value_type<ContainerType>, dg::get_value_type<ContainerTypes>...>
{
    // The reason it is called vdot and not dot is because of the amgiuity when vdot is called
    // with two arguments
    using T = std::invoke_result_t<Functor, dg::get_value_type<ContainerType>, dg::get_value_type<ContainerTypes>...>;
 
    int status = 0;
    if constexpr( std::is_integral_v<T>) // e.g. T = int
    {
        std::array<T, 1> fpe;
        dg::blas1::detail::doDot_fpe( &status, fpe, f, x, xs ...);
        if( fpe[0] - fpe[0] != T(0))
            throw dg::Error(dg::Message(_ping_)
                <<"dg::blas1::vdot (integral type) failed "
                <<"since one of the inputs contains NaN or Inf");
        return fpe[0];
    }
    else
    {
        constexpr size_t N = 3;
        std::array<T, N> fpe;
        dg::blas1::detail::doDot_fpe( &status, fpe, f, x, xs ...);
        for( unsigned u=0; u<N; u++)
        {
            if( fpe[u] - fpe[u] != T(0))
                throw dg::Error(dg::Message(_ping_)
                    <<"dg::blas1::vdot (floating type) failed "
                    <<"since one of the inputs contains NaN or Inf");
        }
        return exblas::cpu::Round(fpe);
    }
}
 
template< class ContainerType1, class ContainerType2>
inline auto dot( const ContainerType1& x, const ContainerType2& y)
{
    if constexpr (std::is_floating_point_v<get_value_type<ContainerType1>> &&
                  std::is_floating_point_v<get_value_type<ContainerType2>>)
    {
        int status = 0;
        std::vector<int64_t> acc = dg::blas1::detail::doDot_superacc( &status,
            x,y);
        if( status != 0)
            throw dg::Error(dg::Message(_ping_)<<"dg::blas1::dot failed "
                <<"since one of the inputs contains NaN or Inf");
        return exblas::cpu::Round(acc.data());
    }
    else
    {
        return dg::blas1::vdot( dg::Product(), x, y);
    }
}
 
 
template< class ContainerType, class OutputType, class BinaryOp, class UnaryOp
    = IDENTITY>
inline OutputType reduce( const ContainerType& x, OutputType zero, BinaryOp
        binary_op, UnaryOp unary_op = UnaryOp())
{
    //init must indeed have the same type as the values of Container since op must be associative
    // The generalization would be a transform_reduce combining subroutine and reduce
    return dg::blas1::detail::doReduce(
            dg::get_tensor_category<ContainerType>(), x, zero, binary_op,
            unary_op);
}
 
template<class ContainerTypeIn, class ContainerTypeOut>
inline void copy( const ContainerTypeIn& source, ContainerTypeOut& target){
    if( std::is_same_v<ContainerTypeIn, ContainerTypeOut> && &source==(const ContainerTypeIn*)&target)
        return;
    dg::blas1::subroutine( dg::equals(), source, target);
}
 
template< class ContainerType, class value_type>
inline void scal( ContainerType& x, value_type alpha)
{
    if( alpha == value_type(1))
        return;
    dg::blas1::subroutine( dg::Scal<value_type>(alpha), x );
}
 
template< class ContainerType, class value_type>
inline void plus( ContainerType& x, value_type alpha)
{
    if( alpha == value_type(0))
        return;
    dg::blas1::subroutine( dg::Plus(alpha), x );
}
 
template< class ContainerType, class ContainerType1, class value_type, class value_type1>
inline void axpby( value_type alpha, const ContainerType1& x, value_type1 beta, ContainerType& y)
{
    if( alpha == value_type(0) ) {
        scal( y, beta);
        return;
    }
    if( std::is_same_v<ContainerType, ContainerType1> && &x==(const ContainerType1*)&y){
        dg::blas1::scal( y, (alpha+beta));
        return;
    }
    dg::blas1::subroutine( dg::Axpby(alpha, beta),  x, y);
}
 
template< class ContainerType, class ContainerType1, class ContainerType2, class value_type, class value_type1, class value_type2>
inline void axpbypgz( value_type alpha, const ContainerType1& x, value_type1 beta, const ContainerType2& y, value_type2 gamma, ContainerType& z)
{
    if( alpha == value_type(0) )
    {
        axpby( beta, y, gamma, z);
        return;
    }
    else if( beta == value_type1(0) )
    {
        axpby( alpha, x, gamma, z);
        return;
    }
    if( std::is_same_v<ContainerType1, ContainerType2> && &x==(const ContainerType1*)&y){
        dg::blas1::axpby( alpha+beta, x, gamma, z);
        return;
    }
    else if( std::is_same_v<ContainerType1, ContainerType> && &x==(const ContainerType1*)&z){
        dg::blas1::axpby( beta, y, alpha+gamma, z);
        return;
    }
    else if( std::is_same_v<ContainerType2, ContainerType> && &y==(const ContainerType2*)&z){
        dg::blas1::axpby( alpha, x, beta+gamma, z);
        return;
    }
    dg::blas1::subroutine( dg::Axpbypgz(alpha, beta, gamma),  x, y, z);
}
 
template< class ContainerType, class ContainerType1, class ContainerType2, class value_type, class value_type1>
inline void axpby( value_type alpha, const ContainerType1& x, value_type1 beta, const ContainerType2& y, ContainerType& z)
{
    dg::blas1::evaluate( z , dg::equals(), dg::PairSum(), alpha, x, beta, y);
}
 
template< class ContainerType, class ContainerType1, class ContainerType2, class value_type, class value_type1>
inline void pointwiseDot( value_type alpha, const ContainerType1& x1, const ContainerType2& x2, value_type1 beta, ContainerType& y)
{
    if( alpha == value_type(0) ) {
        dg::blas1::scal(y, beta);
        return;
    }
    //not sure this is necessary performance-wise, subroutine does allow aliases
    if( std::is_same_v<ContainerType, ContainerType1> && &x1==(const ContainerType1*)&y){
        dg::blas1::subroutine( dg::AxyPby(alpha,beta), x2, y );
 
        return;
    }
    if( std::is_same_v<ContainerType, ContainerType2> && &x2==(const ContainerType2*)&y){
        dg::blas1::subroutine( dg::AxyPby(alpha,beta), x1, y );
 
        return;
    }
    dg::blas1::subroutine( dg::PointwiseDot(alpha,beta), x1, x2, y );
}
 
template< class ContainerType, class ContainerType1, class ContainerType2>
inline void pointwiseDot( const ContainerType1& x1, const ContainerType2& x2, ContainerType& y)
{
    dg::blas1::evaluate( y, dg::equals(), dg::PairSum(), x1,x2);
}
 
template< class ContainerType, class ContainerType1, class ContainerType2, class ContainerType3, class value_type, class value_type1>
inline void pointwiseDot( value_type alpha, const ContainerType1& x1, const ContainerType2& x2, const ContainerType3& x3, value_type1 beta, ContainerType& y)
{
    if( alpha == value_type(0) ) {
        dg::blas1::scal(y, beta);
        return;
    }
    dg::blas1::subroutine( dg::PointwiseDot(alpha,beta), x1, x2, x3, y );
}
 
template< class ContainerType, class ContainerType1, class ContainerType2, class value_type, class value_type1>
inline void pointwiseDivide( value_type alpha, const ContainerType1& x1, const ContainerType2& x2, value_type1 beta, ContainerType& y)
{
    if( alpha == value_type(0) ) {
        dg::blas1::scal(y, beta);
        return;
    }
    if( std::is_same_v<ContainerType, ContainerType1> && &x1==(const ContainerType1*)&y){
        dg::blas1::subroutine( dg::PointwiseDivide(alpha,beta), x2, y );
 
        return;
    }
    dg::blas1::subroutine( dg::PointwiseDivide(alpha, beta), x1, x2, y );
}
 
template< class ContainerType, class ContainerType1, class ContainerType2>
inline void pointwiseDivide( const ContainerType1& x1, const ContainerType2& x2, ContainerType& y)
{
    dg::blas1::evaluate( y, dg::equals(), dg::divides(), x1, x2);
}
 
template<class ContainerType, class ContainerType1, class ContainerType2, class ContainerType3, class ContainerType4, class value_type, class value_type1, class value_type2>
void pointwiseDot(  value_type alpha, const ContainerType1& x1, const ContainerType2& y1,
                    value_type1 beta,  const ContainerType3& x2, const ContainerType4& y2,
                    value_type2 gamma, ContainerType & z)
{
    if( alpha==value_type(0)){
        pointwiseDot( beta, x2,y2, gamma, z);
        return;
    }
    else if( beta==value_type1(0)){
        pointwiseDot( alpha, x1,y1, gamma, z);
        return;
    }
    dg::blas1::subroutine( dg::PointwiseDot2(alpha, beta, gamma), x1, y1, x2, y2, z );
}
 
template< class ContainerType, class ContainerType1, class UnaryOp>
inline void transform( const ContainerType1& x, ContainerType& y, UnaryOp op )
{
    dg::blas1::subroutine( dg::Evaluate<dg::equals, UnaryOp>(dg::equals(),op), y, x);
}
 
template< class ContainerType, class BinarySubroutine, class Functor, class ContainerType0, class ...ContainerTypes>
inline void evaluate( ContainerType& y, BinarySubroutine f, Functor g, const ContainerType0& x0, const ContainerTypes& ...xs)
{
    dg::blas1::subroutine( dg::Evaluate<BinarySubroutine, Functor>(f,g), y, x0, xs...);
}
 
 
namespace detail{
 
template< class T, size_t N, class Functor, class ContainerType, class ...ContainerTypes>
inline void doDot_fpe( int* status, std::array<T,N>& fpe, Functor f,
    const ContainerType& x, const ContainerTypes& ...xs)
{
    static_assert( ( dg::is_vector_v<ContainerType> && ...
                  && dg::is_vector_v<ContainerTypes>),
        "All container types must have a vector data layout (AnyVector)!");
    using vector_type = find_if_t<dg::is_not_scalar, ContainerType, ContainerType, ContainerTypes...>;
    using tensor_category  = get_tensor_category<vector_type>;
    static_assert( ( dg::is_scalar_or_same_base_category<ContainerType, tensor_category>::value &&
              ... && dg::is_scalar_or_same_base_category<ContainerTypes, tensor_category>::value),
        "All container types must be either Scalar or have compatible Vector categories (AnyVector or Same base class)!");
    return doDot_fpe( tensor_category(), status, fpe, f, x, xs ...);
 
}
 
template< class ContainerType1, class ContainerType2>
inline std::vector<int64_t> doDot_superacc( int * status, const ContainerType1& x, const ContainerType2& y)
{
    static_assert( ( dg::is_vector_v<ContainerType1> && dg::is_vector_v<ContainerType2>),
        "All container types must have a vector data layout (AnyVector)!");
    using vector_type = find_if_t<dg::is_not_scalar, ContainerType1, ContainerType1, ContainerType2>;
    using tensor_category  = get_tensor_category<vector_type>;
    static_assert( ( dg::is_scalar_or_same_base_category<ContainerType1, tensor_category>::value
                  && dg::is_scalar_or_same_base_category<ContainerType2, tensor_category>::value),
        "All container types must be either Scalar or have compatible Vector categories (AnyVector or Same base class)!");
    return doDot_superacc( status, x, y, tensor_category());
}
 
}//namespace detail
 
template< class Subroutine, class ContainerType, class ...ContainerTypes>
inline void subroutine( Subroutine f, ContainerType&& x, ContainerTypes&&... xs)
{
    static_assert( ( dg::is_vector_v<ContainerType> &&
               ...&& dg::is_vector_v<ContainerTypes>),
        "All container types must have a vector data layout (AnyVector)!");
    using vector_type = find_if_t<dg::is_not_scalar, ContainerType, ContainerType, ContainerTypes...>;
    using tensor_category  = get_tensor_category<vector_type>;
    static_assert( ( dg::is_scalar_or_same_base_category<ContainerType, tensor_category>::value &&
              ... && dg::is_scalar_or_same_base_category<ContainerTypes, tensor_category>::value),
        "All container types must be either Scalar or have compatible Vector categories (AnyVector or Same base class)!");
    dg::blas1::detail::doSubroutine(tensor_category(), f, std::forward<ContainerType>(x), std::forward<ContainerTypes>(xs)...);
}
 
template< class ContainerType0, class BinarySubroutine, class Functor, class ContainerType1, class ...ContainerTypes>
inline void kronecker( ContainerType0& y, BinarySubroutine f, Functor g, const ContainerType1& x0, const ContainerTypes& ...xs)
{
    static_assert( ( (dg::is_vector_v<ContainerType0> &&
                      dg::is_vector_v<ContainerType1>) &&
                      ... && dg::is_vector_v<ContainerTypes>),
        "All container types must have a vector data layout (AnyVector)!");
    using vector_type = find_if_t<dg::is_not_scalar, ContainerType1, ContainerType1, ContainerTypes...>;
    using tensor_category  = get_tensor_category<vector_type>;
    static_assert( (
            (dg::is_scalar_or_same_base_category<ContainerType0, tensor_category>::value &&
             dg::is_scalar_or_same_base_category<ContainerType1, tensor_category>::value) &&
            ... && dg::is_scalar_or_same_base_category<ContainerTypes, tensor_category>::value),
        "All container types must be either Scalar or have compatible Vector categories (AnyVector or Same base class)!");
    dg::blas1::detail::doKronecker(tensor_category(), y, f, g, x0, xs...);
}
 
}//namespace blas1
 
template<class from_ContainerType, class ContainerType, class ...Params>
inline void assign( const from_ContainerType& from, ContainerType& to, Params&& ... ps)
{
    dg::detail::doAssign<from_ContainerType, ContainerType, Params...>( from, to, get_tensor_category<from_ContainerType>(), get_tensor_category<ContainerType>(), std::forward<Params>(ps)...);
}
 
template<class ContainerType, class from_ContainerType, class ...Params>
inline ContainerType construct( const from_ContainerType& from, Params&& ... ps)
{
    return dg::detail::doConstruct<ContainerType, from_ContainerType, Params...>( from, get_tensor_category<ContainerType>(), get_tensor_category<from_ContainerType>(), std::forward<Params>(ps)...);
}
 
template<class ContainerType, class Functor, class ...ContainerTypes>
auto kronecker( Functor&& f, const ContainerType& x0, const ContainerTypes& ... xs)
{
    using tensor_category  = get_tensor_category<ContainerType>;
    return dg::detail::doKronecker( tensor_category(), std::forward<Functor>(f), x0, xs...);
}
 
 
} //namespace dg