blas2.h

Go to the documentation of this file.

#pragma once
 
#include "backend/tensor_traits.h"
#include "backend/tensor_traits_std.h"
#include "backend/tensor_traits_thrust.h"
#include "backend/blas2_dispatch_scalar.h"
#include "backend/blas2_dispatch_shared.h"
#include "backend/blas2_sparseblockmat.h"
#include "backend/blas2_selfmade.h"
#include "backend/blas2_densematrix.h"
#ifdef MPI_VERSION
#include "backend/blas2_dispatch_mpi.h"
#endif //MPI_VERSION
#include "backend/blas2_dispatch_vector.h"
 
 
namespace dg{
namespace blas2{
 
namespace detail{
 
template< class ContainerType1, class MatrixType, class ContainerType2>
inline std::vector<int64_t> doDot_superacc( int* status, const ContainerType1& x, const MatrixType& m, const ContainerType2& y)
{
    static_assert(
            dg::is_vector_v<ContainerType1> && dg::is_vector_v<MatrixType> &&
            dg::is_vector_v<ContainerType2> ,
        "The container types must have a vector data layout (AnyVector)!");
    //check ContainerTypes: must be either scalar or same base category
    using vector_type = find_if_t<dg::is_not_scalar, ContainerType1, ContainerType1, ContainerType2>;
    using vector_category  = get_tensor_category<vector_type>;
    static_assert(
            dg::is_scalar_or_same_base_category<ContainerType1, vector_category>::value &&
            dg::is_scalar_or_same_base_category<ContainerType2, vector_category>::value,
        "All container types must be either Scalar or have compatible Vector categories (AnyVector or Same base class)!");
    return doDot_superacc( status, x, m, y, get_tensor_category<MatrixType>(), vector_category());
}
 
}//namespace detail
 
template< class ContainerType1, class MatrixType, class ContainerType2>
inline auto dot( const ContainerType1& x, const MatrixType& m, const ContainerType2& y)
{
    if constexpr (std::is_floating_point_v<get_value_type<ContainerType1>> &&
                  std::is_floating_point_v<get_value_type<MatrixType>>   &&
                  std::is_floating_point_v<get_value_type<ContainerType2>>)
    {
        int status = 0;
        std::vector<int64_t> acc = dg::blas2::detail::doDot_superacc( &status,
            x,m,y);
        if( status != 0)
            throw dg::Error(dg::Message(_ping_)<<"dg::blas2::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, m, y);
    }
}
 
template< class MatrixType, class ContainerType>
inline get_value_type<MatrixType> dot( const MatrixType& m, const ContainerType& x)
{
    return dg::blas2::dot( x, m, x);
}
namespace detail{
//resolve tags in two stages: first the matrix and then the container type
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void doSymv( get_value_type<ContainerType1> alpha,
                  MatrixType&& M,
                  const ContainerType1& x,
                  get_value_type<ContainerType1> beta,
                  ContainerType2& y,
                  AnyScalarTag)
{
    dg::blas1::pointwiseDot( alpha, std::forward<MatrixType>(M), x, beta, y);
}
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void doSymv( MatrixType&& M,
                  const ContainerType1& x,
                  ContainerType2& y,
                  AnyScalarTag)
{
    dg::blas1::pointwiseDot( std::forward<MatrixType>(M), x, y);
}
 
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void doSymv( get_value_type<ContainerType1> alpha,
                  MatrixType&& M,
                  const ContainerType1& x,
                  get_value_type<ContainerType1> beta,
                  ContainerType2& y,
                  AnyMatrixTag)
{
    static_assert( std::is_same_v<get_execution_policy<ContainerType1>,
                                  get_execution_policy<ContainerType2>>,
                                "Vector types must have same execution policy");
    // We want to allow double matrix on complex vector ...
    //static_assert( std::is_same_v<get_value_type<ContainerType1>,
    //                              get_value_type<MatrixType>> &&
    //               std::is_same_v<get_value_type<ContainerType2>,
    //                              get_value_type<MatrixType>>,
    //                            "Vector and Matrix types must have same value type");
    static_assert( std::is_same_v<get_tensor_category<ContainerType1>,
                                  get_tensor_category<ContainerType2>>,
                                "Vector types must have same data layout");
    dg::blas2::detail::doSymv( alpha, std::forward<MatrixType>(M), x, beta, y,
            get_tensor_category<MatrixType>(),
            get_tensor_category<ContainerType1>());
}
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void doSymv( MatrixType&& M,
                  const ContainerType1& x,
                  ContainerType2& y,
                  AnyMatrixTag)
{
    static_assert( std::is_same_v<get_execution_policy<ContainerType1>,
                                  get_execution_policy<ContainerType2>>,
                                "Vector types must have same execution policy");
    // We want to allow double matrix on complex vector ...
    //static_assert( std::is_same_v<get_value_type<ContainerType1>,
    //                              get_value_type<MatrixType>> &&
    //               std::is_same_v<get_value_type<ContainerType2>,
    //                              get_value_type<MatrixType>>,
    //                            "Vector and Matrix types must have same value type");
    static_assert( std::is_same_v<get_tensor_category<ContainerType1>,
                                  get_tensor_category<ContainerType2>>,
                                "Vector types must have same data layout");
    dg::blas2::detail::doSymv( std::forward<MatrixType>(M), x, y,
            get_tensor_category<MatrixType>(),
            get_tensor_category<ContainerType1>());
}
template< class FunctorType, class MatrixType, class ContainerType1, class ContainerType2>
inline void doStencil(
                  FunctorType f,
                  MatrixType&& M,
                  const ContainerType1& x,
                  ContainerType2& y,
                  AnyMatrixTag)
{
    static_assert( std::is_same_v<get_execution_policy<ContainerType1>,
                                  get_execution_policy<ContainerType2>>,
                                "Vector types must have same execution policy");
    // We want to allow double matrix on complex vector ...
    //static_assert( std::is_same_v<get_value_type<ContainerType1>,
    //                              get_value_type<MatrixType>> &&
    //               std::is_same_v<get_value_type<ContainerType2>,
    //                              get_value_type<MatrixType>>,
    //                            "Vector and Matrix types must have same value type");
    static_assert( std::is_same_v<get_tensor_category<ContainerType1>,
                                  get_tensor_category<ContainerType2>>,
                                "Vector types must have same data layout");
    dg::blas2::detail::doStencil( f, std::forward<MatrixType>(M), x, y,
            get_tensor_category<MatrixType>(),
            get_tensor_category<ContainerType1>());
}
 
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void doSymv( get_value_type<ContainerType1> alpha,
                  MatrixType&& M,
                  const ContainerType1& x,
                  get_value_type<ContainerType1> beta,
                  ContainerType2& y,
                  NotATensorTag)
{
    // if a compiler error message brings you here then you need to overload
    // the parenthesis operator for your class
    // void operator()( value_type alpha, const ContainerType1& x, value_type beta, ContainerType2& y)
    M(alpha,x,beta,y);
}
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void doSymv( MatrixType&& M,
                  const ContainerType1& x,
                  ContainerType2& y,
                  NotATensorTag)
{
    // if a compiler error message brings you here then you need to overload
    // the parenthesis operator for your class
    // void operator()( const ContainerType1& x, ContainerType2& y)
    M(x,y);
}
 
}//namespace detail
 
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void symv( get_value_type<ContainerType1> alpha,
                  MatrixType&& M,
                  const ContainerType1& x,
                  get_value_type<ContainerType1> beta,
                  ContainerType2& y)
{
    if(alpha == (get_value_type<ContainerType1>)0) {
        dg::blas1::scal( y, beta);
        return;
    }
    dg::blas2::detail::doSymv( alpha, std::forward<MatrixType>(M), x, beta, y, get_tensor_category<MatrixType>());
}
 
 
 
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void symv( MatrixType&& M,
                  const ContainerType1& x,
                  ContainerType2& y)
{
    dg::blas2::detail::doSymv( std::forward<MatrixType>(M), x, y, get_tensor_category<MatrixType>());
}
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void gemv( get_value_type<ContainerType1> alpha,
                  MatrixType&& M,
                  const ContainerType1& x,
                  get_value_type<ContainerType1> beta,
                  ContainerType2& y)
{
    dg::blas2::symv( alpha, std::forward<MatrixType>(M), x, beta, y);
}
 
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void gemv( MatrixType&& M,
                  const ContainerType1& x,
                  ContainerType2& y)
{
    dg::blas2::symv( std::forward<MatrixType>(M), x, y);
}
 
template< class Stencil, class ContainerType, class ...ContainerTypes>
inline void parallel_for( Stencil f, unsigned N, ContainerType&& x, ContainerTypes&&... xs)
{
    // Is the assumption that results are automatically ready on return still true?
    // Do we have to introduce barriers around this function?
    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::blas2::detail::doParallelFor(tensor_category(), f, N, std::forward<ContainerType>(x), std::forward<ContainerTypes>(xs)...);
}
template< class FunctorType, class MatrixType, class ContainerType1, class ContainerType2>
inline void stencil(
                  FunctorType f,
                  MatrixType&& M,
                  const ContainerType1& x,
                  ContainerType2& y)
{
    dg::blas2::detail::doStencil( f, std::forward<MatrixType>(M), x, y, get_tensor_category<MatrixType>());
}
template<class MatrixType, class AnotherMatrixType>
inline void transfer( const MatrixType& x, AnotherMatrixType& y)
{
    dg::blas2::detail::doTransfer( x,y,
            get_tensor_category<MatrixType>(),
            get_tensor_category<AnotherMatrixType>());
}
 
} //namespace blas2
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void apply( get_value_type<ContainerType1> alpha,
                  MatrixType&& M,
                  const ContainerType1& x,
                  get_value_type<ContainerType1> beta,
                  ContainerType2& y)
{
    dg::blas2::symv( alpha, std::forward<MatrixType>(M), x, beta, y);
}
 
template< class MatrixType, class ContainerType1, class ContainerType2>
inline void apply( MatrixType&& M,
                  const ContainerType1& x,
                  ContainerType2& y)
{
    dg::blas2::symv( std::forward<MatrixType>(M), x, y);
}
} //namespace dg