All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
vector_operations.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_CUDA_VECTOR_OPERATIONS_HPP_
#define VIENNACL_LINALG_CUDA_VECTOR_OPERATIONS_HPP_
/* =========================================================================
Copyright (c) 2010-2016, Institute for Microelectronics,
Institute for Analysis and Scientific Computing,
TU Wien.
Portions of this software are copyright by UChicago Argonne, LLC.
-----------------
ViennaCL - The Vienna Computing Library
-----------------
Project Head: Karl Rupp rupp@iue.tuwien.ac.at
(A list of authors and contributors can be found in the manual)
License: MIT (X11), see file LICENSE in the base directory
============================================================================= */
#include <cmath>
#include "viennacl/forwards.h"
#include "viennacl/scalar.hpp"
#include "viennacl/tools/tools.hpp"
#include "viennacl/meta/predicate.hpp"
#include "viennacl/meta/enable_if.hpp"
#include "viennacl/traits/size.hpp"
#include "viennacl/traits/start.hpp"
#include "viennacl/traits/stride.hpp"
#include "viennacl/linalg/cuda/common.hpp"
namespace viennacl
{
namespace linalg
{
namespace cuda
{
//
// Introductory note: By convention, all dimensions are already checked in the dispatcher frontend. No need to double-check again in here!
//
template<typename DestNumericT, typename SrcNumericT>
__global__ void convert_kernel(DestNumericT * dest, unsigned int start_dest, unsigned int inc_dest, unsigned int size_dest,
SrcNumericT const * src, unsigned int start_src, unsigned int inc_src)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size_dest;
i += gridDim.x * blockDim.x)
dest[i*inc_dest+start_dest] = src[i*inc_src+start_src];
}
template<typename DestNumericT, typename SrcNumericT>
void convert(vector_base<DestNumericT> & dest, vector_base<SrcNumericT> const & src)
{
convert_kernel<<<128, 128>>>(viennacl::cuda_arg(dest),
static_cast<unsigned int>(viennacl::traits::start(dest)),
static_cast<unsigned int>(viennacl::traits::stride(dest)),
static_cast<unsigned int>(viennacl::traits::size(dest)),
viennacl::cuda_arg(src),
static_cast<unsigned int>(viennacl::traits::start(src)),
static_cast<unsigned int>(viennacl::traits::stride(src)) );
VIENNACL_CUDA_LAST_ERROR_CHECK("convert_kernel");
}
// gpu scalar
template<typename NumericT>
__global__ void av_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
const NumericT * fac2,
unsigned int options2,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2)
{
NumericT alpha = *fac2;
if (options2 & (1 << 0))
alpha = -alpha;
if (options2 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] / alpha;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] * alpha;
}
}
// cpu scalar
template<typename NumericT>
__global__ void av_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
NumericT fac2,
unsigned int options2,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2)
{
NumericT alpha = fac2;
if (options2 & (1 << 0))
alpha = -alpha;
if (options2 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] / alpha;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] * alpha;
}
}
template<typename NumericT, typename ScalarType1>
void av(vector_base<NumericT> & vec1,
vector_base<NumericT> const & vec2, ScalarType1 const & alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha)
{
typedef NumericT value_type;
unsigned int options_alpha = detail::make_options(len_alpha, reciprocal_alpha, flip_sign_alpha);
value_type data_alpha = alpha;
if (flip_sign_alpha)
data_alpha = -data_alpha;
if (reciprocal_alpha)
data_alpha = static_cast<value_type>(1) / data_alpha;
value_type temporary_alpha = 0;
if (viennacl::is_cpu_scalar<ScalarType1>::value)
temporary_alpha = alpha;
av_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg<value_type>(detail::arg_reference(alpha, temporary_alpha)),
options_alpha,
viennacl::cuda_arg(vec2),
static_cast<unsigned int>(viennacl::traits::start(vec2)),
static_cast<unsigned int>(viennacl::traits::stride(vec2)) );
VIENNACL_CUDA_LAST_ERROR_CHECK("av_kernel");
}
// alpha and beta on GPU
template<typename NumericT>
__global__ void avbv_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
const NumericT * fac2,
unsigned int options2,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2,
const NumericT * fac3,
unsigned int options3,
const NumericT * vec3,
unsigned int start3,
unsigned int inc3)
{
NumericT alpha = *fac2;
if (options2 & (1 << 0))
alpha = -alpha;
NumericT beta = *fac3;
if (options3 & (1 << 0))
beta = -beta;
if (options2 & (1 << 1))
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] * beta;
}
}
else
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] * beta;
}
}
}
// alpha on CPU, beta on GPU
template<typename NumericT>
__global__ void avbv_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
NumericT fac2,
unsigned int options2,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2,
const NumericT * fac3,
unsigned int options3,
const NumericT * vec3,
unsigned int start3,
unsigned int inc3)
{
NumericT alpha = fac2;
if (options2 & (1 << 0))
alpha = -alpha;
NumericT beta = *fac3;
if (options3 & (1 << 0))
beta = -beta;
if (options2 & (1 << 1))
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] * beta;
}
}
else
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] * beta;
}
}
}
// alpha on GPU, beta on CPU
template<typename NumericT>
__global__ void avbv_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
const NumericT * fac2,
unsigned int options2,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2,
NumericT fac3,
unsigned int options3,
const NumericT * vec3,
unsigned int start3,
unsigned int inc3)
{
NumericT alpha = *fac2;
if (options2 & (1 << 0))
alpha = -alpha;
NumericT beta = fac3;
if (options3 & (1 << 0))
beta = -beta;
if (options2 & (1 << 1))
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] * beta;
}
}
else
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] * beta;
}
}
}
// alpha and beta on CPU
template<typename NumericT>
__global__ void avbv_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
NumericT fac2,
unsigned int options2,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2,
NumericT fac3,
unsigned int options3,
const NumericT * vec3,
unsigned int start3,
unsigned int inc3)
{
NumericT alpha = fac2;
if (options2 & (1 << 0))
alpha = -alpha;
NumericT beta = fac3;
if (options3 & (1 << 0))
beta = -beta;
if (options2 & (1 << 1))
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] * beta;
}
}
else
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] * beta;
}
}
}
template<typename NumericT, typename ScalarT1, typename ScalarT2>
void avbv(vector_base<NumericT> & vec1,
vector_base<NumericT> const & vec2, ScalarT1 const & alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha,
vector_base<NumericT> const & vec3, ScalarT2 const & beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
{
typedef NumericT value_type;
unsigned int options_alpha = detail::make_options(len_alpha, reciprocal_alpha, flip_sign_alpha);
value_type data_alpha = alpha;
if (flip_sign_alpha)
data_alpha = -data_alpha;
if (reciprocal_alpha)
data_alpha = static_cast<value_type>(1) / data_alpha;
value_type temporary_alpha = 0;
if (viennacl::is_cpu_scalar<ScalarT1>::value)
temporary_alpha = alpha;
unsigned int options_beta = detail::make_options(len_beta, reciprocal_beta, flip_sign_beta);
value_type temporary_beta = 0;
if (viennacl::is_cpu_scalar<ScalarT2>::value)
temporary_beta = beta;
avbv_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg<value_type>(detail::arg_reference(alpha, temporary_alpha)),
options_alpha,
viennacl::cuda_arg(vec2),
static_cast<unsigned int>(viennacl::traits::start(vec2)),
static_cast<unsigned int>(viennacl::traits::stride(vec2)),
viennacl::cuda_arg<value_type>(detail::arg_reference(beta, temporary_beta)),
options_beta,
viennacl::cuda_arg(vec3),
static_cast<unsigned int>(viennacl::traits::start(vec3)),
static_cast<unsigned int>(viennacl::traits::stride(vec3)) );
VIENNACL_CUDA_LAST_ERROR_CHECK("avbv_kernel");
}
// alpha and beta on GPU
template<typename NumericT>
__global__ void avbv_v_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
const NumericT * fac2,
unsigned int options2,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2,
const NumericT * fac3,
unsigned int options3,
const NumericT * vec3,
unsigned int start3,
unsigned int inc3)
{
NumericT alpha = *fac2;
if (options2 & (1 << 0))
alpha = -alpha;
NumericT beta = *fac3;
if (options3 & (1 << 0))
beta = -beta;
if (options2 & (1 << 1))
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] * beta;
}
}
else
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] * beta;
}
}
}
// alpha on CPU, beta on GPU
template<typename NumericT>
__global__ void avbv_v_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
NumericT fac2,
unsigned int options2,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2,
const NumericT * fac3,
unsigned int options3,
const NumericT * vec3,
unsigned int start3,
unsigned int inc3)
{
NumericT alpha = fac2;
if (options2 & (1 << 0))
alpha = -alpha;
NumericT beta = *fac3;
if (options3 & (1 << 0))
beta = -beta;
if (options2 & (1 << 1))
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] * beta;
}
}
else
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] * beta;
}
}
}
// alpha on GPU, beta on CPU
template<typename NumericT>
__global__ void avbv_v_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
const NumericT * fac2,
unsigned int options2,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2,
NumericT fac3,
unsigned int options3,
const NumericT * vec3,
unsigned int start3,
unsigned int inc3)
{
NumericT alpha = *fac2;
if (options2 & (1 << 0))
alpha = -alpha;
NumericT beta = fac3;
if (options3 & (1 << 0))
beta = -beta;
if (options2 & (1 << 1))
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] * beta;
}
}
else
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] * beta;
}
}
}
// alpha and beta on CPU
template<typename NumericT>
__global__ void avbv_v_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
NumericT fac2,
unsigned int options2,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2,
NumericT fac3,
unsigned int options3,
const NumericT * vec3,
unsigned int start3,
unsigned int inc3)
{
NumericT alpha = fac2;
if (options2 & (1 << 0))
alpha = -alpha;
NumericT beta = fac3;
if (options3 & (1 << 0))
beta = -beta;
if (options2 & (1 << 1))
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] / alpha + vec3[i*inc3+start3] * beta;
}
}
else
{
if (options3 & (1 << 1))
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] / beta;
}
else
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] += vec2[i*inc2+start2] * alpha + vec3[i*inc3+start3] * beta;
}
}
}
template<typename NumericT, typename ScalarT1, typename ScalarT2>
void avbv_v(vector_base<NumericT> & vec1,
vector_base<NumericT> const & vec2, ScalarT1 const & alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha,
vector_base<NumericT> const & vec3, ScalarT2 const & beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
{
typedef NumericT value_type;
unsigned int options_alpha = detail::make_options(len_alpha, reciprocal_alpha, flip_sign_alpha);
value_type data_alpha = alpha;
if (flip_sign_alpha)
data_alpha = -data_alpha;
if (reciprocal_alpha)
data_alpha = static_cast<value_type>(1) / data_alpha;
value_type temporary_alpha = 0;
if (viennacl::is_cpu_scalar<ScalarT1>::value)
temporary_alpha = alpha;
unsigned int options_beta = detail::make_options(len_beta, reciprocal_beta, flip_sign_beta);
value_type temporary_beta = 0;
if (viennacl::is_cpu_scalar<ScalarT2>::value)
temporary_beta = beta;
avbv_v_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg<value_type>(detail::arg_reference(alpha, temporary_alpha)),
options_alpha,
viennacl::cuda_arg(vec2),
static_cast<unsigned int>(viennacl::traits::start(vec2)),
static_cast<unsigned int>(viennacl::traits::stride(vec2)),
viennacl::cuda_arg<value_type>(detail::arg_reference(beta, temporary_beta)),
options_beta,
viennacl::cuda_arg(vec3),
static_cast<unsigned int>(viennacl::traits::start(vec3)),
static_cast<unsigned int>(viennacl::traits::stride(vec3)) );
}
template<typename NumericT>
__global__ void vector_assign_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
unsigned int internal_size1,
NumericT alpha)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = (i < size1) ? alpha : 0;
}
template<typename NumericT, typename ScalarT1>
void vector_assign(vector_base<NumericT> & vec1, ScalarT1 const & alpha, bool up_to_internal_size = false)
{
typedef NumericT value_type;
value_type temporary_alpha = 0;
if (viennacl::is_cpu_scalar<ScalarT1>::value)
temporary_alpha = alpha;
unsigned int size = up_to_internal_size ? static_cast<unsigned int>(vec1.internal_size()) : static_cast<unsigned int>(viennacl::traits::size(vec1));
vector_assign_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
size,
static_cast<unsigned int>(vec1.internal_size()), //Note: Do NOT use traits::internal_size() here, because vector proxies don't require padding.
viennacl::cuda_arg<value_type>(detail::arg_reference(alpha, temporary_alpha)) );
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_assign_kernel");
}
template<typename NumericT>
__global__ void vector_swap_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
NumericT * vec2,
unsigned int start2,
unsigned int inc2)
{
NumericT tmp;
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
{
tmp = vec2[i*inc2+start2];
vec2[i*inc2+start2] = vec1[i*inc1+start1];
vec1[i*inc1+start1] = tmp;
}
}
template<typename NumericT>
void vector_swap(vector_base<NumericT> & vec1, vector_base<NumericT> & vec2)
{
vector_swap_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(vec2),
static_cast<unsigned int>(viennacl::traits::start(vec2)),
static_cast<unsigned int>(viennacl::traits::stride(vec2)) );
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_swap_kernel");
}
template<typename NumericT>
__global__ void element_op_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
NumericT const * vec2,
unsigned int start2,
unsigned int inc2,
NumericT const * vec3,
unsigned int start3,
unsigned int inc3,
unsigned int op_type
)
{
if (op_type == 2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
{
vec1[i*inc1+start1] = pow(vec2[i*inc2+start2], vec3[i*inc3+start3]);
}
}
else if (op_type == 1)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
{
vec1[i*inc1+start1] = vec2[i*inc2+start2] / vec3[i*inc3+start3];
}
}
else if (op_type == 0)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
{
vec1[i*inc1+start1] = vec2[i*inc2+start2] * vec3[i*inc3+start3];
}
}
}
template<typename NumericT>
__global__ void element_op_int_kernel(NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
NumericT const * vec2,
unsigned int start2,
unsigned int inc2,
NumericT const * vec3,
unsigned int start3,
unsigned int inc3,
unsigned int op_type
)
{
if (op_type == 1)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
{
vec1[i*inc1+start1] = vec2[i*inc2+start2] / vec3[i*inc3+start3];
}
}
else if (op_type == 0)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x;
i < size1;
i += gridDim.x * blockDim.x)
{
vec1[i*inc1+start1] = vec2[i*inc2+start2] * vec3[i*inc3+start3];
}
}
}
template<typename NumericT, typename OpT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_binary<OpT> > const & proxy)
{
unsigned int op_type = 2; //0: product, 1: division, 2: power
if (viennacl::is_division<OpT>::value)
op_type = 1;
else if (viennacl::is_product<OpT>::value)
op_type = 0;
element_op_int_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs())),
viennacl::cuda_arg(proxy.rhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.rhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.rhs())),
op_type
);
VIENNACL_CUDA_LAST_ERROR_CHECK("element_op_kernel");
}
template<typename OpT>
void element_op(vector_base<float> & vec1,
vector_expression<const vector_base<float>, const vector_base<float>, op_element_binary<OpT> > const & proxy)
{
unsigned int op_type = 2; //0: product, 1: division, 2: power
if (viennacl::is_division<OpT>::value)
op_type = 1;
else if (viennacl::is_product<OpT>::value)
op_type = 0;
element_op_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs())),
viennacl::cuda_arg(proxy.rhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.rhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.rhs())),
op_type
);
VIENNACL_CUDA_LAST_ERROR_CHECK("element_op_kernel");
}
template<typename OpT>
void element_op(vector_base<double> & vec1,
vector_expression<const vector_base<double>, const vector_base<double>, op_element_binary<OpT> > const & proxy)
{
unsigned int op_type = 2; //0: product, 1: division, 2: power
if (viennacl::is_division<OpT>::value)
op_type = 1;
else if (viennacl::is_product<OpT>::value)
op_type = 0;
element_op_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs())),
viennacl::cuda_arg(proxy.rhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.rhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.rhs())),
op_type
);
VIENNACL_CUDA_LAST_ERROR_CHECK("element_op_kernel");
}
// Note: Trying to automate things with macros or template metaprogramming failed (preprocessor with nvcc did not work as expected), so this is terribly hand-rolled code
// Question (Karl Rupp): Why is CUDA code always such a hassle when trying to use it in a library context?
// acos
template<typename NumericT>
__global__ void vec_element_acos_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = acos(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_acos> > const & proxy)
{
typedef NumericT value_type;
vec_element_acos_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_acos_kernel");
}
// asin
template<typename NumericT>
__global__ void vec_element_asin_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = asin(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_asin> > const & proxy)
{
vec_element_asin_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_asin_kernel");
}
// atan
template<typename NumericT>
__global__ void vec_element_atan_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = atan(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_atan> > const & proxy)
{
vec_element_atan_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_atan_kernel");
}
// ceil
template<typename NumericT>
__global__ void vec_element_ceil_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = ceil(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_ceil> > const & proxy)
{
vec_element_ceil_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_ceil_kernel");
}
// cos
template<typename NumericT>
__global__ void vec_element_cos_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = cos(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_cos> > const & proxy)
{
vec_element_cos_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_cos_kernel");
}
// cosh
template<typename NumericT>
__global__ void vec_element_cosh_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = cosh(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_cosh> > const & proxy)
{
vec_element_cosh_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_cosh_kernel");
}
// exp
template<typename NumericT>
__global__ void vec_element_exp_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = exp(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_exp> > const & proxy)
{
vec_element_exp_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_exp_kernel");
}
// fabs
template<typename NumericT>
__global__ void vec_element_fabs_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = fabs(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_fabs> > const & proxy)
{
vec_element_fabs_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_fabs_kernel");
}
// abs
template<typename NumericT>
__global__ void vec_element_abs_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = abs(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_abs> > const & proxy)
{
vec_element_abs_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_abs_kernel");
}
// floor
template<typename NumericT>
__global__ void vec_element_floor_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = floor(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_floor> > const & proxy)
{
vec_element_floor_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_floor_kernel");
}
// log
template<typename NumericT>
__global__ void vec_element_log_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = log(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_log> > const & proxy)
{
vec_element_log_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_log_kernel");
}
// log10
template<typename NumericT>
__global__ void vec_element_log10_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = log10(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_log10> > const & proxy)
{
vec_element_log10_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_log10_kernel");
}
// sin
template<typename NumericT>
__global__ void vec_element_sin_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = sin(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_sin> > const & proxy)
{
vec_element_sin_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_sin_kernel");
}
// sinh
template<typename NumericT>
__global__ void vec_element_sinh_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = sinh(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_sinh> > const & proxy)
{
vec_element_sinh_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_sinh_kernel");
}
// sqrt
template<typename NumericT>
__global__ void vec_element_sqrt_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = sqrt(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_sqrt> > const & proxy)
{
vec_element_sqrt_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_sqrt_kernel");
}
// tan
template<typename NumericT>
__global__ void vec_element_tan_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = tan(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_tan> > const & proxy)
{
vec_element_tan_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_tan_kernel");
}
// tanh
template<typename NumericT>
__global__ void vec_element_tanh_kernel(
NumericT * vec1, unsigned int start1, unsigned int inc1, unsigned int size1,
NumericT const * vec2, unsigned int start2, unsigned int inc2)
{
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += gridDim.x * blockDim.x)
vec1[i*inc1+start1] = tanh(vec2[i*inc2+start2]);
}
template<typename NumericT>
void element_op(vector_base<NumericT> & vec1,
vector_expression<const vector_base<NumericT>, const vector_base<NumericT>, op_element_unary<op_tanh> > const & proxy)
{
vec_element_tanh_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(proxy.lhs()),
static_cast<unsigned int>(viennacl::traits::start(proxy.lhs())),
static_cast<unsigned int>(viennacl::traits::stride(proxy.lhs()))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vec_element_tanh_kernel");
}
template<typename NumericT>
__global__ void inner_prod_kernel(const NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
const NumericT * vec2,
unsigned int start2,
unsigned int inc2,
unsigned int size2,
NumericT * group_buffer)
{
__shared__ NumericT tmp_buffer[128];
unsigned int group_start1 = (blockIdx.x * size1) / (gridDim.x) * inc1 + start1;
unsigned int group_start2 = (blockIdx.x * size2) / (gridDim.x) * inc2 + start2;
unsigned int group_size1 = ((blockIdx.x + 1) * size1) / (gridDim.x)
- ( blockIdx.x * size1) / (gridDim.x);
NumericT tmp = 0;
for (unsigned int i = threadIdx.x; i < group_size1; i += blockDim.x)
tmp += vec1[i*inc1+group_start1] * vec2[i*inc2+group_start2];
tmp_buffer[threadIdx.x] = tmp;
// parallel reduction
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] += tmp_buffer[threadIdx.x+stride];
}
if (threadIdx.x == 0)
group_buffer[blockIdx.x] = tmp_buffer[0];
}
// sums the array 'vec1' and writes to result. Makes use of a single work-group only.
template<typename NumericT>
__global__ void vector_sum_kernel_floats(
const NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
unsigned int option, //0: use fmax, 1: just sum, 2: sum and return sqrt of sum
NumericT * result)
{
__shared__ NumericT tmp_buffer[128];
NumericT thread_sum = 0;
for (unsigned int i = threadIdx.x; i<size1; i += blockDim.x)
{
if (option > 0)
thread_sum += vec1[i*inc1+start1];
else
thread_sum = fmax(thread_sum, fabs(vec1[i*inc1+start1]));
}
tmp_buffer[threadIdx.x] = thread_sum;
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
{
if (option > 0)
tmp_buffer[threadIdx.x] += tmp_buffer[threadIdx.x + stride];
else
tmp_buffer[threadIdx.x] = fmax(tmp_buffer[threadIdx.x], tmp_buffer[threadIdx.x + stride]);
}
}
if (threadIdx.x == 0)
{
if (option == 2)
*result = sqrt(tmp_buffer[0]);
else
*result = tmp_buffer[0];
}
}
template<typename NumericT>
__global__ void vector_sum_kernel_integers(
const NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
unsigned int option, //0: use max, 1: just sum
NumericT * result)
{
__shared__ NumericT tmp_buffer[128];
NumericT thread_sum = 0;
for (unsigned int i = threadIdx.x; i<size1; i += blockDim.x)
{
if (option > 0)
thread_sum += vec1[i*inc1+start1];
else
thread_sum = thread_sum > abs(vec1[i*inc1+start1]) ? thread_sum : abs(vec1[i*inc1+start1]);
}
tmp_buffer[threadIdx.x] = thread_sum;
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
{
if (option > 0)
tmp_buffer[threadIdx.x] += tmp_buffer[threadIdx.x + stride];
else
tmp_buffer[threadIdx.x] = tmp_buffer[threadIdx.x] > tmp_buffer[threadIdx.x + stride] ? tmp_buffer[threadIdx.x] : tmp_buffer[threadIdx.x + stride];
}
}
if (threadIdx.x == 0)
*result = tmp_buffer[0];
}
template<typename NumericT>
__global__ void vector_sum_kernel_unsigned_integers(
const NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
unsigned int option, //0: use max, 1: just sum
NumericT * result)
{
__shared__ NumericT tmp_buffer[128];
NumericT thread_sum = 0;
for (unsigned int i = threadIdx.x; i<size1; i += blockDim.x)
{
if (option > 0)
thread_sum += vec1[i*inc1+start1];
else
thread_sum = (thread_sum > vec1[i*inc1+start1]) ? thread_sum : vec1[i*inc1+start1];
}
tmp_buffer[threadIdx.x] = thread_sum;
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
{
if (option > 0)
tmp_buffer[threadIdx.x] += tmp_buffer[threadIdx.x + stride];
else
tmp_buffer[threadIdx.x] = tmp_buffer[threadIdx.x] > tmp_buffer[threadIdx.x + stride] ? tmp_buffer[threadIdx.x] : tmp_buffer[threadIdx.x + stride];
}
}
if (threadIdx.x == 0)
*result = tmp_buffer[0];
}
namespace detail
{
struct vector_sum_kernel_launcher_integers
{
template<typename NumericT, typename ScalarT>
static void apply(vector_base<NumericT> const & temp,
unsigned int option,
ScalarT & result)
{
typedef NumericT value_type;
vector_sum_kernel_integers<<<1, 128>>>(viennacl::cuda_arg(temp),
static_cast<unsigned int>(viennacl::traits::start(temp)),
static_cast<unsigned int>(viennacl::traits::stride(temp)),
static_cast<unsigned int>(viennacl::traits::size(temp)),
static_cast<unsigned int>(option),
viennacl::cuda_arg(result) );
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_sum_kernel");
}
};
struct vector_sum_kernel_launcher_unsigned_integers
{
template<typename NumericT, typename ScalarT>
static void apply(vector_base<NumericT> const & temp,
unsigned int option,
ScalarT & result)
{
typedef NumericT value_type;
vector_sum_kernel_unsigned_integers<<<1, 128>>>(viennacl::cuda_arg(temp),
static_cast<unsigned int>(viennacl::traits::start(temp)),
static_cast<unsigned int>(viennacl::traits::stride(temp)),
static_cast<unsigned int>(viennacl::traits::size(temp)),
static_cast<unsigned int>(option),
viennacl::cuda_arg(result) );
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_sum_kernel");
}
};
struct vector_sum_kernel_launcher_floats
{
template<typename NumericT, typename ScalarT>
static void apply(vector_base<NumericT> const & temp,
unsigned int option,
ScalarT & result)
{
typedef NumericT value_type;
vector_sum_kernel_floats<<<1, 128>>>(viennacl::cuda_arg(temp),
static_cast<unsigned int>(viennacl::traits::start(temp)),
static_cast<unsigned int>(viennacl::traits::stride(temp)),
static_cast<unsigned int>(viennacl::traits::size(temp)),
static_cast<unsigned int>(option),
viennacl::cuda_arg(result) );
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_sum_kernel");
}
};
template<typename NumericT>
struct vector_sum_kernel_launcher : public vector_sum_kernel_launcher_integers {};
template<>
struct vector_sum_kernel_launcher<unsigned char> : public vector_sum_kernel_launcher_unsigned_integers {};
template<>
struct vector_sum_kernel_launcher<unsigned short> : public vector_sum_kernel_launcher_unsigned_integers {};
template<>
struct vector_sum_kernel_launcher<unsigned int> : public vector_sum_kernel_launcher_unsigned_integers {};
template<>
struct vector_sum_kernel_launcher<unsigned long> : public vector_sum_kernel_launcher_unsigned_integers {};
template<>
struct vector_sum_kernel_launcher<float> : public vector_sum_kernel_launcher_floats {};
template<>
struct vector_sum_kernel_launcher<double> : public vector_sum_kernel_launcher_floats {};
}
//implementation of inner product:
//namespace {
template<typename NumericT, typename ScalarT>
void inner_prod_impl(vector_base<NumericT> const & vec1,
vector_base<NumericT> const & vec2,
ScalarT & result)
{
typedef NumericT value_type;
static const unsigned int work_groups = 128;
static viennacl::vector<value_type> temp(work_groups);
inner_prod_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(vec2),
static_cast<unsigned int>(viennacl::traits::start(vec2)),
static_cast<unsigned int>(viennacl::traits::stride(vec2)),
static_cast<unsigned int>(viennacl::traits::size(vec2)),
viennacl::cuda_arg(temp)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("inner_prod_kernel");
detail::vector_sum_kernel_launcher<NumericT>::apply(temp, 1, result);
}
template<typename NumericT>
void inner_prod_cpu(vector_base<NumericT> const & vec1,
vector_base<NumericT> const & vec2,
NumericT & result)
{
typedef NumericT value_type;
const unsigned int work_groups = 128;
viennacl::vector<value_type> temp(work_groups);
inner_prod_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(vec2),
static_cast<unsigned int>(viennacl::traits::start(vec2)),
static_cast<unsigned int>(viennacl::traits::stride(vec2)),
static_cast<unsigned int>(viennacl::traits::size(vec2)),
viennacl::cuda_arg(temp)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("inner_prod_kernel");
// Now copy partial results from GPU back to CPU and run reduction there:
std::vector<value_type> temp_cpu(work_groups);
viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
result = 0;
for (typename std::vector<value_type>::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
result += *it;
}
#define VIENNACL_MDOT_WORKGROUP_SIZE 128
#define VIENNACL_MDOT_WORKGROUP_NUM 128
// M = 2:
template<typename NumericT>
__global__ void inner_prod_2_kernel(const NumericT *x, unsigned int startx, unsigned int stridex, unsigned int sizex,
const NumericT *y0, unsigned int start0, unsigned int stride0,
const NumericT *y1, unsigned int start1, unsigned int stride1,
NumericT *group_results)
{
__shared__ NumericT tmp_buffer[2*VIENNACL_MDOT_WORKGROUP_SIZE];
unsigned int entries_per_thread = (sizex - 1) / (blockDim.x * gridDim.x) + 1;
unsigned int vec_start_index = blockIdx.x * blockDim.x * entries_per_thread;
unsigned int vec_stop_index = min((blockIdx.x + 1) * blockDim.x * entries_per_thread, sizex); // don't go beyond size of x
NumericT entry_x = 0;
NumericT group_sum0 = 0;
NumericT group_sum1 = 0;
for (unsigned int i = vec_start_index + threadIdx.x; i < vec_stop_index; i += blockDim.x) {
entry_x = x[i * stridex + startx]; // load only once from global memory!
group_sum0 += entry_x * y0[i * stride0 + start0];
group_sum1 += entry_x * y1[i * stride1 + start1];
}
tmp_buffer[threadIdx.x] = group_sum0;
tmp_buffer[threadIdx.x + blockDim.x] = group_sum1;
// parallel reduction
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2) {
__syncthreads();
if (threadIdx.x < stride) {
tmp_buffer[threadIdx.x ] += tmp_buffer[threadIdx.x+stride ];
tmp_buffer[threadIdx.x + blockDim.x] += tmp_buffer[threadIdx.x+stride + blockDim.x];
}
}
// write result of group to group_results
if (threadIdx.x == 0) {
group_results[blockIdx.x] = tmp_buffer[0];
group_results[blockIdx.x + gridDim.x] = tmp_buffer[blockDim.x];
}
}
// M = 3:
template<typename NumericT>
__global__ void inner_prod_3_kernel(const NumericT *x, unsigned int startx, unsigned int stridex, unsigned int sizex,
const NumericT *y0, unsigned int start0, unsigned int stride0,
const NumericT *y1, unsigned int start1, unsigned int stride1,
const NumericT *y2, unsigned int start2, unsigned int stride2,
NumericT *group_results)
{
__shared__ NumericT tmp_buffer[3*VIENNACL_MDOT_WORKGROUP_SIZE];
unsigned int entries_per_thread = (sizex - 1) / (blockDim.x * gridDim.x) + 1;
unsigned int vec_start_index = blockIdx.x * blockDim.x * entries_per_thread;
unsigned int vec_stop_index = min((blockIdx.x + 1) * blockDim.x * entries_per_thread, sizex); // don't go beyond vec size
NumericT entry_x = 0;
NumericT group_sum0 = 0;
NumericT group_sum1 = 0;
NumericT group_sum2 = 0;
for (unsigned int i = vec_start_index + threadIdx.x; i < vec_stop_index; i += blockDim.x) {
entry_x = x[i * stridex + startx]; // load only once from global memory!
group_sum0 += entry_x * y0[i * stride0 + start0];
group_sum1 += entry_x * y1[i * stride1 + start1];
group_sum2 += entry_x * y2[i * stride2 + start2];
}
tmp_buffer[threadIdx.x] = group_sum0;
tmp_buffer[threadIdx.x + blockDim.x] = group_sum1;
tmp_buffer[threadIdx.x + 2 * blockDim.x] = group_sum2;
// parallel reduction
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2) {
__syncthreads();
if (threadIdx.x < stride) {
tmp_buffer[threadIdx.x ] += tmp_buffer[threadIdx.x+stride ];
tmp_buffer[threadIdx.x + blockDim.x] += tmp_buffer[threadIdx.x+stride + blockDim.x];
tmp_buffer[threadIdx.x + 2 * blockDim.x] += tmp_buffer[threadIdx.x+stride + 2 * blockDim.x];
}
}
// write result of group to group_results
if (threadIdx.x == 0) {
group_results[blockIdx.x ] = tmp_buffer[0];
group_results[blockIdx.x + gridDim.x] = tmp_buffer[ blockDim.x];
group_results[blockIdx.x + 2 * gridDim.x] = tmp_buffer[2 * blockDim.x];
}
}
// M = 4:
template<typename NumericT>
__global__ void inner_prod_4_kernel(const NumericT *x, unsigned int startx, unsigned int stridex, unsigned int sizex,
const NumericT *y0, unsigned int start0, unsigned int stride0,
const NumericT *y1, unsigned int start1, unsigned int stride1,
const NumericT *y2, unsigned int start2, unsigned int stride2,
const NumericT *y3, unsigned int start3, unsigned int stride3,
NumericT *group_results)
{
__shared__ NumericT tmp_buffer[4*VIENNACL_MDOT_WORKGROUP_SIZE];
unsigned int entries_per_thread = (sizex - 1) / (blockDim.x * gridDim.x) + 1;
unsigned int vec_start_index = blockIdx.x * blockDim.x * entries_per_thread;
unsigned int vec_stop_index = min((blockIdx.x + 1) * blockDim.x * entries_per_thread, sizex); // don't go beyond vec size
NumericT entry_x = 0;
NumericT group_sum0 = 0;
NumericT group_sum1 = 0;
NumericT group_sum2 = 0;
NumericT group_sum3 = 0;
for (unsigned int i = vec_start_index + threadIdx.x; i < vec_stop_index; i += blockDim.x) {
entry_x = x[i * stridex + startx]; // load only once from global memory!
group_sum0 += entry_x * y0[i * stride0 + start0];
group_sum1 += entry_x * y1[i * stride1 + start1];
group_sum2 += entry_x * y2[i * stride2 + start2];
group_sum3 += entry_x * y3[i * stride3 + start3];
}
tmp_buffer[threadIdx.x] = group_sum0;
tmp_buffer[threadIdx.x + blockDim.x] = group_sum1;
tmp_buffer[threadIdx.x + 2 * blockDim.x] = group_sum2;
tmp_buffer[threadIdx.x + 3 * blockDim.x] = group_sum3;
// parallel reduction
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2) {
__syncthreads();
if (threadIdx.x < stride) {
tmp_buffer[threadIdx.x ] += tmp_buffer[threadIdx.x+stride ];
tmp_buffer[threadIdx.x + blockDim.x] += tmp_buffer[threadIdx.x+stride + blockDim.x];
tmp_buffer[threadIdx.x + 2 * blockDim.x] += tmp_buffer[threadIdx.x+stride + 2 * blockDim.x];
tmp_buffer[threadIdx.x + 3 * blockDim.x] += tmp_buffer[threadIdx.x+stride + 3 * blockDim.x];
}
}
// write result of group to group_results
if (threadIdx.x == 0) {
group_results[blockIdx.x ] = tmp_buffer[0];
group_results[blockIdx.x + gridDim.x] = tmp_buffer[ blockDim.x];
group_results[blockIdx.x + 2 * gridDim.x] = tmp_buffer[2 * blockDim.x];
group_results[blockIdx.x + 3 * gridDim.x] = tmp_buffer[3 * blockDim.x];
}
}
// M = 8:
template<typename NumericT>
__global__ void inner_prod_8_kernel(const NumericT *x, unsigned int startx, unsigned int stridex, unsigned int sizex,
const NumericT *y0, unsigned int start0, unsigned int stride0,
const NumericT *y1, unsigned int start1, unsigned int stride1,
const NumericT *y2, unsigned int start2, unsigned int stride2,
const NumericT *y3, unsigned int start3, unsigned int stride3,
const NumericT *y4, unsigned int start4, unsigned int stride4,
const NumericT *y5, unsigned int start5, unsigned int stride5,
const NumericT *y6, unsigned int start6, unsigned int stride6,
const NumericT *y7, unsigned int start7, unsigned int stride7,
NumericT *group_results)
{
__shared__ NumericT tmp_buffer[8*VIENNACL_MDOT_WORKGROUP_SIZE];
unsigned int entries_per_thread = (sizex - 1) / (blockDim.x * gridDim.x) + 1;
unsigned int vec_start_index = blockIdx.x * blockDim.x * entries_per_thread;
unsigned int vec_stop_index = min((blockIdx.x + 1) * blockDim.x * entries_per_thread, sizex); // don't go beyond vec size
NumericT entry_x = 0;
NumericT group_sum0 = 0;
NumericT group_sum1 = 0;
NumericT group_sum2 = 0;
NumericT group_sum3 = 0;
NumericT group_sum4 = 0;
NumericT group_sum5 = 0;
NumericT group_sum6 = 0;
NumericT group_sum7 = 0;
for (unsigned int i = vec_start_index + threadIdx.x; i < vec_stop_index; i += blockDim.x) {
entry_x = x[i * stridex + startx]; // load only once from global memory!
group_sum0 += entry_x * y0[i * stride0 + start0];
group_sum1 += entry_x * y1[i * stride1 + start1];
group_sum2 += entry_x * y2[i * stride2 + start2];
group_sum3 += entry_x * y3[i * stride3 + start3];
group_sum4 += entry_x * y4[i * stride4 + start4];
group_sum5 += entry_x * y5[i * stride5 + start5];
group_sum6 += entry_x * y6[i * stride6 + start6];
group_sum7 += entry_x * y7[i * stride7 + start7];
}
tmp_buffer[threadIdx.x] = group_sum0;
tmp_buffer[threadIdx.x + blockDim.x] = group_sum1;
tmp_buffer[threadIdx.x + 2 * blockDim.x] = group_sum2;
tmp_buffer[threadIdx.x + 3 * blockDim.x] = group_sum3;
tmp_buffer[threadIdx.x + 4 * blockDim.x] = group_sum4;
tmp_buffer[threadIdx.x + 5 * blockDim.x] = group_sum5;
tmp_buffer[threadIdx.x + 6 * blockDim.x] = group_sum6;
tmp_buffer[threadIdx.x + 7 * blockDim.x] = group_sum7;
// parallel reduction
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2) {
__syncthreads();
if (threadIdx.x < stride) {
tmp_buffer[threadIdx.x ] += tmp_buffer[threadIdx.x+stride ];
tmp_buffer[threadIdx.x + blockDim.x] += tmp_buffer[threadIdx.x+stride + blockDim.x];
tmp_buffer[threadIdx.x + 2 * blockDim.x] += tmp_buffer[threadIdx.x+stride + 2 * blockDim.x];
tmp_buffer[threadIdx.x + 3 * blockDim.x] += tmp_buffer[threadIdx.x+stride + 3 * blockDim.x];
tmp_buffer[threadIdx.x + 4 * blockDim.x] += tmp_buffer[threadIdx.x+stride + 4 * blockDim.x];
tmp_buffer[threadIdx.x + 5 * blockDim.x] += tmp_buffer[threadIdx.x+stride + 5 * blockDim.x];
tmp_buffer[threadIdx.x + 6 * blockDim.x] += tmp_buffer[threadIdx.x+stride + 6 * blockDim.x];
tmp_buffer[threadIdx.x + 7 * blockDim.x] += tmp_buffer[threadIdx.x+stride + 7 * blockDim.x];
}
}
// write result of group to group_results
if (threadIdx.x == 0) {
group_results[blockIdx.x ] = tmp_buffer[0];
group_results[blockIdx.x + gridDim.x] = tmp_buffer[ blockDim.x];
group_results[blockIdx.x + 2 * gridDim.x] = tmp_buffer[2 * blockDim.x];
group_results[blockIdx.x + 3 * gridDim.x] = tmp_buffer[3 * blockDim.x];
group_results[blockIdx.x + 4 * gridDim.x] = tmp_buffer[4 * blockDim.x];
group_results[blockIdx.x + 5 * gridDim.x] = tmp_buffer[5 * blockDim.x];
group_results[blockIdx.x + 6 * gridDim.x] = tmp_buffer[6 * blockDim.x];
group_results[blockIdx.x + 7 * gridDim.x] = tmp_buffer[7 * blockDim.x];
}
}
// sums the array 'vec1' and writes to result. Makes use of a single work-group only.
template<typename NumericT>
__global__ void vector_multi_sum_kernel(
NumericT const * vec1,
NumericT * result,
unsigned int start_result,
unsigned int inc_result)
{
__shared__ NumericT tmp_buffer[VIENNACL_MDOT_WORKGROUP_SIZE];
tmp_buffer[threadIdx.x] = vec1[threadIdx.x + blockIdx.x * VIENNACL_MDOT_WORKGROUP_SIZE];
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] += tmp_buffer[threadIdx.x + stride];
}
if (threadIdx.x == 0)
result[start_result + inc_result * blockIdx.x] = tmp_buffer[0];
}
template<typename NumericT>
void inner_prod_impl(vector_base<NumericT> const & x,
vector_tuple<NumericT> const & vec_tuple,
vector_base<NumericT> & result)
{
typedef NumericT value_type;
static viennacl::vector<value_type> temp(8 * VIENNACL_MDOT_WORKGROUP_NUM);
vcl_size_t current_index = 0;
while (vec_tuple.const_size() > current_index)
{
switch (vec_tuple.const_size() - current_index)
{
case 7:
case 6:
case 5:
case 4:
{
vector_base<NumericT> const & y0 = vec_tuple.const_at(current_index);
vector_base<NumericT> const & y1 = vec_tuple.const_at(current_index + 1);
vector_base<NumericT> const & y2 = vec_tuple.const_at(current_index + 2);
vector_base<NumericT> const & y3 = vec_tuple.const_at(current_index + 3);
inner_prod_4_kernel<<<VIENNACL_MDOT_WORKGROUP_NUM,
VIENNACL_MDOT_WORKGROUP_SIZE>>>( viennacl::cuda_arg(x),
static_cast<unsigned int>(viennacl::traits::start(x)),
static_cast<unsigned int>(viennacl::traits::stride(x)),
static_cast<unsigned int>(viennacl::traits::size(x)),
viennacl::cuda_arg(y0),
static_cast<unsigned int>(viennacl::traits::start(y0)),
static_cast<unsigned int>(viennacl::traits::stride(y0)),
viennacl::cuda_arg(y1),
static_cast<unsigned int>(viennacl::traits::start(y1)),
static_cast<unsigned int>(viennacl::traits::stride(y1)),
viennacl::cuda_arg(y2),
static_cast<unsigned int>(viennacl::traits::start(y2)),
static_cast<unsigned int>(viennacl::traits::stride(y2)),
viennacl::cuda_arg(y3),
static_cast<unsigned int>(viennacl::traits::start(y3)),
static_cast<unsigned int>(viennacl::traits::stride(y3)),
viennacl::cuda_arg(temp)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("inner_prod_4_kernel");
vector_multi_sum_kernel<<<4, VIENNACL_MDOT_WORKGROUP_NUM>>>(viennacl::cuda_arg(temp),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start(result) + viennacl::traits::stride(result) * current_index),
static_cast<unsigned int>(viennacl::traits::stride(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_multi_sum_kernel");
}
current_index += 4;
break;
case 3:
{
vector_base<NumericT> const & y0 = vec_tuple.const_at(current_index);
vector_base<NumericT> const & y1 = vec_tuple.const_at(current_index + 1);
vector_base<NumericT> const & y2 = vec_tuple.const_at(current_index + 2);
inner_prod_3_kernel<<<VIENNACL_MDOT_WORKGROUP_NUM,
VIENNACL_MDOT_WORKGROUP_SIZE>>>( viennacl::cuda_arg(x),
static_cast<unsigned int>(viennacl::traits::start(x)),
static_cast<unsigned int>(viennacl::traits::stride(x)),
static_cast<unsigned int>(viennacl::traits::size(x)),
viennacl::cuda_arg(y0),
static_cast<unsigned int>(viennacl::traits::start(y0)),
static_cast<unsigned int>(viennacl::traits::stride(y0)),
viennacl::cuda_arg(y1),
static_cast<unsigned int>(viennacl::traits::start(y1)),
static_cast<unsigned int>(viennacl::traits::stride(y1)),
viennacl::cuda_arg(y2),
static_cast<unsigned int>(viennacl::traits::start(y2)),
static_cast<unsigned int>(viennacl::traits::stride(y2)),
viennacl::cuda_arg(temp)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("inner_prod_3_kernel");
vector_multi_sum_kernel<<<3, VIENNACL_MDOT_WORKGROUP_NUM>>>(viennacl::cuda_arg(temp),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start(result) + viennacl::traits::stride(result) * current_index),
static_cast<unsigned int>(viennacl::traits::stride(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_multi_sum_kernel");
}
current_index += 3;
break;
case 2:
{
vector_base<NumericT> const & y0 = vec_tuple.const_at(current_index);
vector_base<NumericT> const & y1 = vec_tuple.const_at(current_index + 1);
inner_prod_2_kernel<<<VIENNACL_MDOT_WORKGROUP_NUM,
VIENNACL_MDOT_WORKGROUP_SIZE>>>( viennacl::cuda_arg(x),
static_cast<unsigned int>(viennacl::traits::start(x)),
static_cast<unsigned int>(viennacl::traits::stride(x)),
static_cast<unsigned int>(viennacl::traits::size(x)),
viennacl::cuda_arg(y0),
static_cast<unsigned int>(viennacl::traits::start(y0)),
static_cast<unsigned int>(viennacl::traits::stride(y0)),
viennacl::cuda_arg(y1),
static_cast<unsigned int>(viennacl::traits::start(y1)),
static_cast<unsigned int>(viennacl::traits::stride(y1)),
viennacl::cuda_arg(temp)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("inner_prod_2_kernel");
vector_multi_sum_kernel<<<2, VIENNACL_MDOT_WORKGROUP_NUM>>>(viennacl::cuda_arg(temp),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start(result) + viennacl::traits::stride(result) * current_index),
static_cast<unsigned int>(viennacl::traits::stride(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_multi_sum_kernel");
}
current_index += 2;
break;
case 1:
{
vector_base<NumericT> const & y0 = vec_tuple.const_at(current_index);
inner_prod_kernel<<<128, 128>>>(viennacl::cuda_arg(x),
static_cast<unsigned int>(viennacl::traits::start(x)),
static_cast<unsigned int>(viennacl::traits::stride(x)),
static_cast<unsigned int>(viennacl::traits::size(x)),
viennacl::cuda_arg(y0),
static_cast<unsigned int>(viennacl::traits::start(y0)),
static_cast<unsigned int>(viennacl::traits::stride(y0)),
static_cast<unsigned int>(viennacl::traits::size(y0)),
viennacl::cuda_arg(temp)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("inner_prod_kernel");
vector_multi_sum_kernel<<<1, 128>>>(viennacl::cuda_arg(temp),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start(result) + viennacl::traits::stride(result) * current_index),
static_cast<unsigned int>(viennacl::traits::stride(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_multi_sum_kernel");
}
current_index += 1;
break;
default:
{
vector_base<NumericT> const & y0 = vec_tuple.const_at(current_index);
vector_base<NumericT> const & y1 = vec_tuple.const_at(current_index + 1);
vector_base<NumericT> const & y2 = vec_tuple.const_at(current_index + 2);
vector_base<NumericT> const & y3 = vec_tuple.const_at(current_index + 3);
vector_base<NumericT> const & y4 = vec_tuple.const_at(current_index + 4);
vector_base<NumericT> const & y5 = vec_tuple.const_at(current_index + 5);
vector_base<NumericT> const & y6 = vec_tuple.const_at(current_index + 6);
vector_base<NumericT> const & y7 = vec_tuple.const_at(current_index + 7);
inner_prod_8_kernel<<<VIENNACL_MDOT_WORKGROUP_NUM,
VIENNACL_MDOT_WORKGROUP_SIZE>>>( viennacl::cuda_arg(x),
static_cast<unsigned int>(viennacl::traits::start(x)),
static_cast<unsigned int>(viennacl::traits::stride(x)),
static_cast<unsigned int>(viennacl::traits::size(x)),
viennacl::cuda_arg(y0),
static_cast<unsigned int>(viennacl::traits::start(y0)),
static_cast<unsigned int>(viennacl::traits::stride(y0)),
viennacl::cuda_arg(y1),
static_cast<unsigned int>(viennacl::traits::start(y1)),
static_cast<unsigned int>(viennacl::traits::stride(y1)),
viennacl::cuda_arg(y2),
static_cast<unsigned int>(viennacl::traits::start(y2)),
static_cast<unsigned int>(viennacl::traits::stride(y2)),
viennacl::cuda_arg(y3),
static_cast<unsigned int>(viennacl::traits::start(y3)),
static_cast<unsigned int>(viennacl::traits::stride(y3)),
viennacl::cuda_arg(y4),
static_cast<unsigned int>(viennacl::traits::start(y4)),
static_cast<unsigned int>(viennacl::traits::stride(y4)),
viennacl::cuda_arg(y5),
static_cast<unsigned int>(viennacl::traits::start(y5)),
static_cast<unsigned int>(viennacl::traits::stride(y5)),
viennacl::cuda_arg(y6),
static_cast<unsigned int>(viennacl::traits::start(y6)),
static_cast<unsigned int>(viennacl::traits::stride(y6)),
viennacl::cuda_arg(y7),
static_cast<unsigned int>(viennacl::traits::start(y7)),
static_cast<unsigned int>(viennacl::traits::stride(y7)),
viennacl::cuda_arg(temp)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("inner_prod_8_kernel");
vector_multi_sum_kernel<<<8, VIENNACL_MDOT_WORKGROUP_NUM>>>(viennacl::cuda_arg(temp),
viennacl::cuda_arg(result),
static_cast<unsigned int>(viennacl::traits::start(result) + viennacl::traits::stride(result) * current_index),
static_cast<unsigned int>(viennacl::traits::stride(result))
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_multi_sum_kernel");
}
current_index += 8;
break;
}
}
}
#undef VIENNACL_MDOT_WORKGROUP_NUM
#undef VIENNACL_MDOT_WORKGROUP_SIZE
template<typename NumericT>
__global__ void norm_kernel_floats(
const NumericT * vec,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
unsigned int norm_selector,
NumericT * group_buffer)
{
__shared__ NumericT tmp_buffer[128];
NumericT tmp = (norm_selector > 2) ? vec[start1] : 0;
unsigned int work_per_thread = (size1 - 1) / (gridDim.x * blockDim.x) + 1;
unsigned int group_start = blockIdx.x * work_per_thread * blockDim.x;
unsigned int group_stop = (blockIdx.x + 1) * work_per_thread * blockDim.x;
group_stop = (group_stop > size1) ? size1 : group_stop;
if (norm_selector == 1) //norm_1
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp += fabs(vec[i*inc1 + start1]);
}
else if (norm_selector == 2) //norm_2
{
NumericT vec_entry = 0;
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
{
vec_entry = vec[i*inc1 + start1];
tmp += vec_entry * vec_entry;
}
}
else if (norm_selector == 0) //norm_inf
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp = fmax(fabs(vec[i*inc1 + start1]), tmp);
}
else if (norm_selector == 3) //min
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp = (vec[i*inc1 + start1] < tmp) ? vec[i*inc1 + start1] : tmp;
}
else if (norm_selector == 4) //max
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp = (vec[i*inc1 + start1] > tmp) ? vec[i*inc1 + start1] : tmp;
}
tmp_buffer[threadIdx.x] = tmp;
if (norm_selector == 1 || norm_selector == 2) //parallel reduction for norm_1 or norm_2:
{
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] += tmp_buffer[threadIdx.x+stride];
}
}
else if (norm_selector == 3)
{
//min:
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] = (tmp_buffer[threadIdx.x+stride] < tmp_buffer[threadIdx.x]) ? tmp_buffer[threadIdx.x+stride] : tmp_buffer[threadIdx.x];
}
}
else if (norm_selector == 4)
{
//max:
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] = (tmp_buffer[threadIdx.x+stride] > tmp_buffer[threadIdx.x]) ? tmp_buffer[threadIdx.x+stride] : tmp_buffer[threadIdx.x];
}
}
else
{
//norm_inf:
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] = fmax(tmp_buffer[threadIdx.x], tmp_buffer[threadIdx.x+stride]);
}
}
if (threadIdx.x == 0)
group_buffer[blockIdx.x] = tmp_buffer[0];
}
template<typename NumericT>
__global__ void norm_kernel_integers(
const NumericT * vec,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
unsigned int norm_selector,
NumericT * group_buffer)
{
__shared__ NumericT tmp_buffer[128];
NumericT tmp = (norm_selector > 2) ? vec[start1] : 0;
unsigned int work_per_thread = (size1 - 1) / (gridDim.x * blockDim.x) + 1;
unsigned int group_start = blockIdx.x * work_per_thread * blockDim.x;
unsigned int group_stop = (blockIdx.x + 1) * work_per_thread * blockDim.x;
group_stop = (group_stop > size1) ? size1 : group_stop;
if (norm_selector == 1) //norm_1
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp += abs(vec[i*inc1 + start1]);
}
else if (norm_selector == 0) //norm_inf
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp = (tmp > abs(vec[i*inc1 + start1])) ? tmp : abs(vec[i*inc1 + start1]);
}
else if (norm_selector == 3) //min
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp = (vec[i*inc1 + start1] < tmp) ? vec[i*inc1 + start1] : tmp;
}
else if (norm_selector == 4) //max
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp = (vec[i*inc1 + start1] > tmp) ? vec[i*inc1 + start1] : tmp;
}
tmp_buffer[threadIdx.x] = tmp;
if (norm_selector == 1 || norm_selector == 2) //parallel reduction for norm_1 or norm_2:
{
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] += tmp_buffer[threadIdx.x+stride];
}
}
else if (norm_selector == 3)
{
//min:
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] = (tmp_buffer[threadIdx.x+stride] < tmp_buffer[threadIdx.x]) ? tmp_buffer[threadIdx.x+stride] : tmp_buffer[threadIdx.x];
}
}
else if (norm_selector == 4)
{
//max:
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] = (tmp_buffer[threadIdx.x+stride] > tmp_buffer[threadIdx.x]) ? tmp_buffer[threadIdx.x+stride] : tmp_buffer[threadIdx.x];
}
}
else
{
//norm_inf:
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] = (tmp_buffer[threadIdx.x] > tmp_buffer[threadIdx.x+stride]) ? tmp_buffer[threadIdx.x] : tmp_buffer[threadIdx.x+stride];
}
}
if (threadIdx.x == 0)
group_buffer[blockIdx.x] = tmp_buffer[0];
}
template<typename NumericT>
__global__ void norm_kernel_unsigned_integers(
const NumericT * vec,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
unsigned int norm_selector,
NumericT * group_buffer)
{
__shared__ NumericT tmp_buffer[128];
NumericT tmp = (norm_selector > 2) ? vec[start1] : 0;
unsigned int work_per_thread = (size1 - 1) / (gridDim.x * blockDim.x) + 1;
unsigned int group_start = blockIdx.x * work_per_thread * blockDim.x;
unsigned int group_stop = (blockIdx.x + 1) * work_per_thread * blockDim.x;
group_stop = (group_stop > size1) ? size1 : group_stop;
if (norm_selector == 1) //norm_1
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp += vec[i*inc1 + start1];
}
else if (norm_selector == 0) //norm_inf
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp = (tmp > vec[i*inc1 + start1]) ? tmp : vec[i*inc1 + start1];
}
else if (norm_selector == 3) //min
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp = (vec[i*inc1 + start1] < tmp) ? vec[i*inc1 + start1] : tmp;
}
else if (norm_selector == 4) //max
{
for (unsigned int i = group_start + threadIdx.x; i < group_stop; i += blockDim.x)
tmp = (vec[i*inc1 + start1] > tmp) ? vec[i*inc1 + start1] : tmp;
}
tmp_buffer[threadIdx.x] = tmp;
if (norm_selector == 1 || norm_selector == 2) //parallel reduction for norm_1 or norm_2:
{
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] += tmp_buffer[threadIdx.x+stride];
}
}
else if (norm_selector == 3)
{
//min:
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] = (tmp_buffer[threadIdx.x+stride] < tmp_buffer[threadIdx.x]) ? tmp_buffer[threadIdx.x+stride] : tmp_buffer[threadIdx.x];
}
}
else if (norm_selector == 4)
{
//max:
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] = (tmp_buffer[threadIdx.x+stride] > tmp_buffer[threadIdx.x]) ? tmp_buffer[threadIdx.x+stride] : tmp_buffer[threadIdx.x];
}
}
else
{
//norm_inf:
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
tmp_buffer[threadIdx.x] = (tmp_buffer[threadIdx.x] > tmp_buffer[threadIdx.x+stride]) ? tmp_buffer[threadIdx.x] : tmp_buffer[threadIdx.x+stride];
}
}
if (threadIdx.x == 0)
group_buffer[blockIdx.x] = tmp_buffer[0];
}
namespace detail
{
struct norm_kernel_launcher_integers
{
template<typename NumericT>
static void apply(vector_base<NumericT> const & vec1,
vector_base<NumericT> & temp,
unsigned int option)
{
norm_kernel_integers<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
static_cast<unsigned int>(option),
viennacl::cuda_arg(temp)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("norm_kernel");
}
};
struct norm_kernel_launcher_unsigned_integers
{
template<typename NumericT>
static void apply(vector_base<NumericT> const & vec1,
vector_base<NumericT> & temp,
unsigned int option)
{
norm_kernel_unsigned_integers<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
static_cast<unsigned int>(option),
viennacl::cuda_arg(temp)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("norm_kernel");
}
};
struct norm_kernel_launcher_floats
{
template<typename NumericT>
static void apply(vector_base<NumericT> const & vec1,
vector_base<NumericT> & temp,
unsigned int option)
{
norm_kernel_floats<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
static_cast<unsigned int>(option),
viennacl::cuda_arg(temp)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("norm_kernel");
}
};
template<typename NumericT>
struct norm_kernel_launcher : public norm_kernel_launcher_integers {};
template<>
struct norm_kernel_launcher<unsigned char> : public norm_kernel_launcher_unsigned_integers {};
template<>
struct norm_kernel_launcher<unsigned short> : public norm_kernel_launcher_unsigned_integers {};
template<>
struct norm_kernel_launcher<unsigned int> : public norm_kernel_launcher_unsigned_integers {};
template<>
struct norm_kernel_launcher<unsigned long> : public norm_kernel_launcher_unsigned_integers {};
template<>
struct norm_kernel_launcher<float> : public norm_kernel_launcher_floats {};
template<>
struct norm_kernel_launcher<double> : public norm_kernel_launcher_floats {};
}
template<typename NumericT>
void norm_1_impl(vector_base<NumericT> const & vec1,
scalar<NumericT> & result)
{
typedef NumericT value_type;
vcl_size_t work_groups = 128;
viennacl::vector<value_type> temp(work_groups);
detail::norm_kernel_launcher<NumericT>::apply(vec1, temp, 1);
detail::vector_sum_kernel_launcher<NumericT>::apply(temp, 1, result);
}
template<typename NumericT>
void norm_1_cpu(vector_base<NumericT> const & vec1,
NumericT & result)
{
typedef NumericT value_type;
vcl_size_t work_groups = 128;
viennacl::vector<value_type> temp(work_groups);
detail::norm_kernel_launcher<NumericT>::apply(vec1, temp, 1);
// Now copy partial results from GPU back to CPU and run reduction there:
std::vector<value_type> temp_cpu(work_groups);
viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
result = 0;
for (typename std::vector<value_type>::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
result += *it;
}
template<typename NumericT>
void norm_2_impl(vector_base<NumericT> const & vec1,
scalar<NumericT> & result)
{
typedef NumericT value_type;
vcl_size_t work_groups = 128;
viennacl::vector<value_type> temp(work_groups);
detail::norm_kernel_launcher<NumericT>::apply(vec1, temp, 2);
detail::vector_sum_kernel_launcher<NumericT>::apply(temp, 2, result);
}
template<typename NumericT>
void norm_2_cpu(vector_base<NumericT> const & vec1,
NumericT & result)
{
typedef NumericT value_type;
vcl_size_t work_groups = 128;
viennacl::vector<value_type> temp(work_groups);
detail::norm_kernel_launcher<NumericT>::apply(vec1, temp, 2);
std::vector<value_type> temp_cpu(work_groups);
viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
result = 0;
for (typename std::vector<value_type>::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
result += *it;
result = std::sqrt(result);
}
template<typename NumericT>
void norm_inf_impl(vector_base<NumericT> const & vec1,
scalar<NumericT> & result)
{
typedef NumericT value_type;
vcl_size_t work_groups = 128;
viennacl::vector<value_type> temp(work_groups);
detail::norm_kernel_launcher<NumericT>::apply(vec1, temp, 0);
detail::vector_sum_kernel_launcher<NumericT>::apply(temp, 0, result);
}
template<typename NumericT>
void norm_inf_cpu(vector_base<NumericT> const & vec1,
NumericT & result)
{
typedef NumericT value_type;
vcl_size_t work_groups = 128;
viennacl::vector<value_type> temp(work_groups);
detail::norm_kernel_launcher<NumericT>::apply(vec1, temp, 0);
std::vector<value_type> temp_cpu(work_groups);
viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
result = 0;
for (typename std::vector<value_type>::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
result = std::max(result, *it);
}
// second reduction stage for min() and max()
template<typename NumericT>
__global__ void vector_maxmin_kernel(
const NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
unsigned int option, //0: use max, 1: use min
NumericT * result)
{
__shared__ NumericT tmp_buffer[128];
NumericT thread_minmax = vec1[start1];
for (unsigned int i = threadIdx.x; i<size1; i += blockDim.x)
{
if (option > 0) //min
thread_minmax = (vec1[i*inc1+start1] < thread_minmax) ? vec1[i*inc1+start1] : thread_minmax;
else
thread_minmax = (vec1[i*inc1+start1] > thread_minmax) ? vec1[i*inc1+start1] : thread_minmax;
}
tmp_buffer[threadIdx.x] = thread_minmax;
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
{
if (option > 0) //min
tmp_buffer[threadIdx.x] = (tmp_buffer[threadIdx.x + stride] < tmp_buffer[threadIdx.x]) ? tmp_buffer[threadIdx.x + stride] : tmp_buffer[threadIdx.x];
else
tmp_buffer[threadIdx.x] = (tmp_buffer[threadIdx.x + stride] > tmp_buffer[threadIdx.x]) ? tmp_buffer[threadIdx.x + stride] : tmp_buffer[threadIdx.x];
}
}
if (threadIdx.x == 0)
*result = tmp_buffer[0];
}
template<typename NumericT>
void max_impl(vector_base<NumericT> const & vec1,
scalar<NumericT> & result)
{
typedef NumericT value_type;
vcl_size_t work_groups = 128;
viennacl::vector<value_type> temp(work_groups, viennacl::traits::context(vec1));
detail::norm_kernel_launcher<NumericT>::apply(vec1, temp, 4);
vector_maxmin_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
static_cast<unsigned int>(0),
viennacl::cuda_arg(result)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_maxmin_kernel");
}
template<typename NumericT>
void max_cpu(vector_base<NumericT> const & vec1,
NumericT & result)
{
typedef NumericT value_type;
vcl_size_t work_groups = 128;
viennacl::vector<value_type> temp(work_groups, viennacl::traits::context(vec1));
detail::norm_kernel_launcher<NumericT>::apply(vec1, temp, 4);
std::vector<value_type> temp_cpu(work_groups);
viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
result = temp[0];
for (typename std::vector<value_type>::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
result = std::max(result, *it);
}
template<typename NumericT>
void min_impl(vector_base<NumericT> const & vec1,
scalar<NumericT> & result)
{
typedef NumericT value_type;
vcl_size_t work_groups = 128;
viennacl::vector<value_type> temp(work_groups, viennacl::traits::context(vec1));
detail::norm_kernel_launcher<NumericT>::apply(vec1, temp, 3);
vector_maxmin_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
static_cast<unsigned int>(1),
viennacl::cuda_arg(result)
);
VIENNACL_CUDA_LAST_ERROR_CHECK("vector_maxmin_kernel");
}
template<typename NumericT>
void min_cpu(vector_base<NumericT> const & vec1,
NumericT & result)
{
typedef NumericT value_type;
vcl_size_t work_groups = 128;
viennacl::vector<value_type> temp(work_groups, viennacl::traits::context(vec1));
detail::norm_kernel_launcher<NumericT>::apply(vec1, temp, 3);
std::vector<value_type> temp_cpu(work_groups);
viennacl::fast_copy(temp.begin(), temp.end(), temp_cpu.begin());
result = temp[0];
for (typename std::vector<value_type>::const_iterator it = temp_cpu.begin(); it != temp_cpu.end(); ++it)
result = std::min(result, *it);
}
template<typename NumericT>
void sum_impl(vector_base<NumericT> const & vec1,
scalar<NumericT> & result)
{
typedef NumericT value_type;
viennacl::vector<NumericT> all_ones = viennacl::scalar_vector<NumericT>(vec1.size(), NumericT(1), viennacl::traits::context(vec1));
viennacl::linalg::cuda::inner_prod_impl(vec1, all_ones, result);
}
template<typename NumericT>
void sum_cpu(vector_base<NumericT> const & vec1,
NumericT & result)
{
typedef NumericT value_type;
viennacl::vector<NumericT> all_ones = viennacl::scalar_vector<NumericT>(vec1.size(), NumericT(1), viennacl::traits::context(vec1));
viennacl::linalg::cuda::inner_prod_cpu(vec1, all_ones, result);
}
//index_norm_inf:
// fixes the problem of not having (f)abs available in a consistent manner
template<typename NumericT>
__device__ NumericT cuda_abs(NumericT val) { return (val < 0) ? -val : val; }
__device__ inline unsigned long cuda_abs(unsigned long val) { return val; }
__device__ inline unsigned int cuda_abs(unsigned int val) { return val; }
__device__ inline unsigned short cuda_abs(unsigned short val) { return val; }
__device__ inline unsigned char cuda_abs(unsigned char val) { return val; }
template<typename NumericT>
__global__ void index_norm_inf_kernel(const NumericT * vec,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
unsigned int * result)
{
__shared__ NumericT float_buffer[128];
__shared__ unsigned int index_buffer[128];
float_buffer[threadIdx.x] = 0;
index_buffer[threadIdx.x] = 0;
//step 1: fill buffer:
NumericT cur_max = NumericT(0);
NumericT tmp;
for (unsigned int i = threadIdx.x; i < size1; i += blockDim.x)
{
tmp = vec[i*inc1+start1];
tmp = cuda_abs(tmp);
if (cur_max < tmp)
{
float_buffer[threadIdx.x] = tmp;
index_buffer[threadIdx.x] = i;
cur_max = tmp;
}
}
//step 2: parallel reduction:
for (unsigned int stride = blockDim.x/2; stride > 0; stride /= 2)
{
__syncthreads();
if (threadIdx.x < stride)
{
//find the first occurring index
if (float_buffer[threadIdx.x] < float_buffer[threadIdx.x+stride])
{
index_buffer[threadIdx.x] = index_buffer[threadIdx.x+stride];
float_buffer[threadIdx.x] = float_buffer[threadIdx.x+stride];
}
}
}
if (threadIdx.x == 0)
*result = index_buffer[0];
}
//This function should return a CPU scalar, otherwise statements like
// vcl_rhs[index_norm_inf(vcl_rhs)]
// are ambiguous
template<typename NumericT>
vcl_size_t index_norm_inf(vector_base<NumericT> const & vec1)
{
typedef NumericT value_type;
viennacl::backend::mem_handle h;
viennacl::backend::memory_create(h, sizeof(unsigned int), viennacl::traits::context(vec1));
index_norm_inf_kernel<<<1, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg<unsigned int>(h)
//reinterpret_cast<unsigned int *>(h.cuda_handle().get())
);
VIENNACL_CUDA_LAST_ERROR_CHECK("index_norm_inf_kernel");
unsigned int ret = 0;
viennacl::backend::memory_read(h, 0, sizeof(unsigned int), &ret);
return static_cast<vcl_size_t>(ret);
}
template<typename NumericT>
__global__ void plane_rotation_kernel(
NumericT * vec1,
unsigned int start1,
unsigned int inc1,
unsigned int size1,
NumericT * vec2,
unsigned int start2,
unsigned int inc2,
unsigned int size2,
NumericT alpha,
NumericT beta)
{
NumericT tmp1 = 0;
NumericT tmp2 = 0;
for (unsigned int i = blockDim.x * blockIdx.x + threadIdx.x; i < size1; i += blockDim.x * gridDim.x)
{
tmp1 = vec1[i*inc1+start1];
tmp2 = vec2[i*inc2+start2];
vec1[i*inc1+start1] = alpha * tmp1 + beta * tmp2;
vec2[i*inc2+start2] = alpha * tmp2 - beta * tmp1;
}
}
template<typename NumericT>
void plane_rotation(vector_base<NumericT> & vec1,
vector_base<NumericT> & vec2,
NumericT alpha, NumericT beta)
{
typedef NumericT value_type;
value_type temporary_alpha = 0;
if (viennacl::is_cpu_scalar<value_type>::value)
temporary_alpha = alpha;
value_type temporary_beta = 0;
if (viennacl::is_cpu_scalar<value_type>::value)
temporary_beta = beta;
plane_rotation_kernel<<<128, 128>>>(viennacl::cuda_arg(vec1),
static_cast<unsigned int>(viennacl::traits::start(vec1)),
static_cast<unsigned int>(viennacl::traits::stride(vec1)),
static_cast<unsigned int>(viennacl::traits::size(vec1)),
viennacl::cuda_arg(vec2),
static_cast<unsigned int>(viennacl::traits::start(vec2)),
static_cast<unsigned int>(viennacl::traits::stride(vec2)),
static_cast<unsigned int>(viennacl::traits::size(vec2)),
viennacl::cuda_arg<value_type>(detail::arg_reference(alpha, temporary_alpha)),
viennacl::cuda_arg<value_type>(detail::arg_reference(beta, temporary_beta)) );
VIENNACL_CUDA_LAST_ERROR_CHECK("plane_rotation_kernel");
}
template<typename NumericT>
__global__ void scan_kernel_1(NumericT const *X,
unsigned int startX,
unsigned int incX,
unsigned int sizeX,
NumericT *Y,
unsigned int startY,
unsigned int incY,
unsigned int scan_offset,
NumericT *carries) // 0 for inclusive scan, 1 for exclusive
{
__shared__ NumericT shared_buffer[256];
NumericT my_value;
unsigned int work_per_thread = (sizeX - 1) / (gridDim.x * blockDim.x) + 1;
unsigned int block_start = work_per_thread * blockDim.x * blockIdx.x;
unsigned int block_stop = work_per_thread * blockDim.x * (blockIdx.x + 1);
unsigned int block_offset = 0;
// run scan on each section
for (unsigned int i = block_start + threadIdx.x; i < block_stop; i += blockDim.x)
{
// load data:
my_value = (i < sizeX) ? X[i * incX + startX] : 0;
// inclusive scan in shared buffer:
for(unsigned int stride = 1; stride < blockDim.x; stride *= 2)
{
__syncthreads();
shared_buffer[threadIdx.x] = my_value;
__syncthreads();
if (threadIdx.x >= stride)
my_value += shared_buffer[threadIdx.x - stride];
}
__syncthreads();
shared_buffer[threadIdx.x] = my_value;
__syncthreads();
// exclusive scan requires us to write a zero value at the beginning of each block
if (scan_offset > 0)
my_value = (threadIdx.x > 0) ? shared_buffer[threadIdx.x - 1] : 0;
// write to output array
if (i < sizeX)
Y[i * incY + startY] = block_offset + my_value;
block_offset += shared_buffer[blockDim.x-1];
}
// write carry:
if (threadIdx.x == 0)
carries[blockIdx.x] = block_offset;
}
// exclusive-scan of carries
template<typename NumericT>
__global__ void scan_kernel_2(NumericT *carries)
{
__shared__ NumericT shared_buffer[256];
// load data:
NumericT my_carry = carries[threadIdx.x];
// exclusive scan in shared buffer:
for(unsigned int stride = 1; stride < blockDim.x; stride *= 2)
{
__syncthreads();
shared_buffer[threadIdx.x] = my_carry;
__syncthreads();
if (threadIdx.x >= stride)
my_carry += shared_buffer[threadIdx.x - stride];
}
__syncthreads();
shared_buffer[threadIdx.x] = my_carry;
__syncthreads();
// write to output array
carries[threadIdx.x] = (threadIdx.x > 0) ? shared_buffer[threadIdx.x - 1] : 0;
}
template<typename NumericT>
__global__ void scan_kernel_3(NumericT *Y,
unsigned int startY,
unsigned int incY,
unsigned int sizeY,
NumericT const *carries)
{
unsigned int work_per_thread = (sizeY - 1) / (gridDim.x * blockDim.x) + 1;
unsigned int block_start = work_per_thread * blockDim.x * blockIdx.x;
unsigned int block_stop = work_per_thread * blockDim.x * (blockIdx.x + 1);
__shared__ NumericT shared_offset;
if (threadIdx.x == 0)
shared_offset = carries[blockIdx.x];
__syncthreads();
// add offset to each element in the block:
for (unsigned int i = block_start + threadIdx.x; i < block_stop; i += blockDim.x)
if (i < sizeY)
Y[i * incY + startY] += shared_offset;
}
namespace detail
{
template<typename NumericT>
void scan_impl(vector_base<NumericT> const & input,
vector_base<NumericT> & output,
bool is_inclusive)
{
vcl_size_t block_num = 128;
vcl_size_t threads_per_block = 128;
viennacl::backend::mem_handle cuda_carries;
viennacl::backend::memory_create(cuda_carries, sizeof(NumericT)*block_num, viennacl::traits::context(input));
// First step: Scan within each thread group and write carries
scan_kernel_1<<<block_num, threads_per_block>>>(viennacl::cuda_arg(input),
static_cast<unsigned int>(viennacl::traits::start(input)),
static_cast<unsigned int>(viennacl::traits::stride(input)),
static_cast<unsigned int>(viennacl::traits::size(input)),
viennacl::cuda_arg(output),
static_cast<unsigned int>(viennacl::traits::start(output)),
static_cast<unsigned int>(viennacl::traits::stride(output)),
static_cast<unsigned int>(is_inclusive ? 0 : 1),
viennacl::cuda_arg<NumericT>(cuda_carries)
);
// Second step: Compute offset for each thread group (exclusive scan for each thread group)
scan_kernel_2<<<1, block_num>>>(viennacl::cuda_arg<NumericT>(cuda_carries));
// Third step: Offset each thread group accordingly
scan_kernel_3<<<block_num, threads_per_block>>>(viennacl::cuda_arg(output),
static_cast<unsigned int>(viennacl::traits::start(output)),
static_cast<unsigned int>(viennacl::traits::stride(output)),
static_cast<unsigned int>(viennacl::traits::size(output)),
viennacl::cuda_arg<NumericT>(cuda_carries)
);
}
}
template<typename NumericT>
void inclusive_scan(vector_base<NumericT> const & input,
vector_base<NumericT> & output)
{
detail::scan_impl(input, output, true);
}
template<typename NumericT>
void exclusive_scan(vector_base<NumericT> const & input,
vector_base<NumericT> & output)
{
detail::scan_impl(input, output, false);
}
} //namespace cuda
} //namespace linalg
} //namespace viennacl
#endif
viennacl::vector_tuple::const_size
vcl_size_t const_size() const
Definition: vector.hpp:1143
viennacl::linalg::cuda::convert_kernel
__global__ void convert_kernel(DestNumericT *dest, unsigned int start_dest, unsigned int inc_dest, unsigned int size_dest, SrcNumericT const *src, unsigned int start_src, unsigned int inc_src)
Definition: vector_operations.hpp:48
viennacl::linalg::cuda::vec_element_abs_kernel
__global__ void vec_element_abs_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1255
viennacl::linalg::cuda::detail::make_options
unsigned int make_options(vcl_size_t length, bool reciprocal, bool flip_sign)
Definition: common.hpp:160
viennacl::linalg::cuda::vector_assign
void vector_assign(vector_base< NumericT > &vec1, ScalarT1 const &alpha, bool up_to_internal_size=false)
Assign a constant value to a vector (-range/-slice)
Definition: vector_operations.hpp:803
viennacl::scalar
This class represents a single scalar value on the GPU and behaves mostly like a built-in scalar type...
Definition: forwards.h:227
viennacl::linalg::cuda::vector_sum_kernel_unsigned_integers
__global__ void vector_sum_kernel_unsigned_integers(const NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, unsigned int option, NumericT *result)
Definition: vector_operations.hpp:1612
viennacl::linalg::cuda::norm_2_cpu
void norm_2_cpu(vector_base< NumericT > const &vec1, NumericT &result)
Computes the l^2-norm of a vector - implementation.
Definition: vector_operations.hpp:2651
viennacl::linalg::cuda::convert
void convert(matrix_base< DestNumericT > &mat1, matrix_base< SrcNumericT > const &mat2)
Definition: matrix_operations.hpp:57
viennacl::linalg::cuda::norm_kernel_floats
__global__ void norm_kernel_floats(const NumericT *vec, unsigned int start1, unsigned int inc1, unsigned int size1, unsigned int norm_selector, NumericT *group_buffer)
Definition: vector_operations.hpp:2238
viennacl::traits::stride1
result_of::size_type< matrix_base< NumericT > >::type stride1(matrix_base< NumericT > const &s)
Definition: stride.hpp:55
viennacl::linalg::cuda::inner_prod_4_kernel
__global__ void inner_prod_4_kernel(const NumericT *x, unsigned int startx, unsigned int stridex, unsigned int sizex, const NumericT *y0, unsigned int start0, unsigned int stride0, const NumericT *y1, unsigned int start1, unsigned int stride1, const NumericT *y2, unsigned int start2, unsigned int stride2, const NumericT *y3, unsigned int start3, unsigned int stride3, NumericT *group_results)
Definition: vector_operations.hpp:1891
viennacl::linalg::cuda::vec_element_asin_kernel
__global__ void vec_element_asin_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1074
size.hpp
Generic size and resize functionality for different vector and matrix types.
viennacl::linalg::cuda::plane_rotation_kernel
__global__ void plane_rotation_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT *vec2, unsigned int start2, unsigned int inc2, unsigned int size2, NumericT alpha, NumericT beta)
Definition: vector_operations.hpp:2996
start.hpp
Extracts the underlying OpenCL start index handle from a vector, a matrix, an expression etc...
viennacl::linalg::cuda::inner_prod_3_kernel
__global__ void inner_prod_3_kernel(const NumericT *x, unsigned int startx, unsigned int stridex, unsigned int sizex, const NumericT *y0, unsigned int start0, unsigned int stride0, const NumericT *y1, unsigned int start1, unsigned int stride1, const NumericT *y2, unsigned int start2, unsigned int stride2, NumericT *group_results)
Definition: vector_operations.hpp:1846
tools.hpp
Various little tools used here and there in ViennaCL.
viennacl::traits::internal_size1
vcl_size_t internal_size1(matrix_base< NumericT > const &mat)
Helper routine for obtaining the internal number of entries per row of a ViennaCL matrix...
Definition: size.hpp:386
viennacl::traits::size1
vcl_size_t size1(MatrixType const &mat)
Generic routine for obtaining the number of rows of a matrix (ViennaCL, uBLAS, etc.)
Definition: size.hpp:163
viennacl::linalg::cuda::av
void av(vector_base< NumericT > &vec1, vector_base< NumericT > const &vec2, ScalarType1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha)
Definition: vector_operations.hpp:144
viennacl::linalg::cuda::vector_multi_sum_kernel
__global__ void vector_multi_sum_kernel(NumericT const *vec1, NumericT *result, unsigned int start_result, unsigned int inc_result)
Definition: vector_operations.hpp:2017
viennacl::linalg::cuda::vector_maxmin_kernel
__global__ void vector_maxmin_kernel(const NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, unsigned int option, NumericT *result)
Definition: vector_operations.hpp:2722
viennacl::linalg::cuda::vec_element_fabs_kernel
__global__ void vec_element_fabs_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1230
viennacl::traits::stride
result_of::size_type< viennacl::vector_base< T > >::type stride(viennacl::vector_base< T > const &s)
Definition: stride.hpp:45
viennacl::linalg::cuda::max_impl
void max_impl(vector_base< NumericT > const &vec1, scalar< NumericT > &result)
Computes the maximum of a vector, both reduction stages run on the GPU.
Definition: vector_operations.hpp:2765
forwards.h
This file provides the forward declarations for the main types used within ViennaCL.
viennacl::traits::start1
result_of::size_type< T >::type start1(T const &obj)
Definition: start.hpp:65
viennacl::linalg::cuda::norm_kernel_unsigned_integers
__global__ void norm_kernel_unsigned_integers(const NumericT *vec, unsigned int start1, unsigned int inc1, unsigned int size1, unsigned int norm_selector, NumericT *group_buffer)
Definition: vector_operations.hpp:2415
stride.hpp
Determines row and column increments for matrices and matrix proxies.
viennacl::backend::memory_read
void memory_read(mem_handle const &src_buffer, vcl_size_t src_offset, vcl_size_t bytes_to_read, void *ptr, bool async=false)
Reads data from a buffer back to main RAM.
Definition: memory.hpp:261
viennacl::linalg::cuda::exclusive_scan
void exclusive_scan(vector_base< NumericT > const &input, vector_base< NumericT > &output)
This function implements an exclusive scan using CUDA.
Definition: vector_operations.hpp:3239
viennacl::linalg::cuda::norm_1_cpu
void norm_1_cpu(vector_base< NumericT > const &vec1, NumericT &result)
Computes the l^1-norm of a vector.
Definition: vector_operations.hpp:2605
viennacl::linalg::detail::max
T max(const T &lhs, const T &rhs)
Maximum.
Definition: util.hpp:59
VIENNACL_MDOT_WORKGROUP_SIZE
#define VIENNACL_MDOT_WORKGROUP_SIZE
Definition: vector_operations.hpp:1803
viennacl::vector_expression
An expression template class that represents a binary operation that yields a vector.
Definition: forwards.h:239
viennacl::linalg::cuda::scan_kernel_2
__global__ void scan_kernel_2(NumericT *carries)
Definition: vector_operations.hpp:3121
viennacl::traits::size2
result_of::size_type< MatrixType >::type size2(MatrixType const &mat)
Generic routine for obtaining the number of columns of a matrix (ViennaCL, uBLAS, etc...
Definition: size.hpp:201
viennacl::linalg::cuda::vec_element_sqrt_kernel
__global__ void vec_element_sqrt_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1412
viennacl::linalg::cuda::avbv_v_kernel
__global__ void avbv_v_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, const NumericT *fac2, unsigned int options2, const NumericT *vec2, unsigned int start2, unsigned int inc2, const NumericT *fac3, unsigned int options3, const NumericT *vec3, unsigned int start3, unsigned int inc3)
Definition: vector_operations.hpp:483
viennacl::linalg::cuda::max_cpu
void max_cpu(vector_base< NumericT > const &vec1, NumericT &result)
Computes the maximum of a vector, first reduction stage on the GPU, second stage on the CPU...
Definition: vector_operations.hpp:2793
viennacl::linalg::cuda::vec_element_tanh_kernel
__global__ void vec_element_tanh_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1464
NumericT
float NumericT
Definition: bisect.cpp:40
VIENNACL_MDOT_WORKGROUP_NUM
#define VIENNACL_MDOT_WORKGROUP_NUM
Definition: vector_operations.hpp:1804
viennacl::linalg::cuda::vector_swap
void vector_swap(vector_base< NumericT > &vec1, vector_base< NumericT > &vec2)
Swaps the contents of two vectors, data is copied.
Definition: vector_operations.hpp:853
viennacl::linalg::cuda::index_norm_inf_kernel
__global__ void index_norm_inf_kernel(const NumericT *vec, unsigned int start1, unsigned int inc1, unsigned int size1, unsigned int *result)
Definition: vector_operations.hpp:2917
viennacl::linalg::cuda::detail::scan_impl
void scan_impl(vector_base< NumericT > const &input, vector_base< NumericT > &output, bool is_inclusive)
Worker routine for scan routines.
Definition: vector_operations.hpp:3181
viennacl::linalg::cuda::inner_prod_kernel
__global__ void inner_prod_kernel(const NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, const NumericT *vec2, unsigned int start2, unsigned int inc2, unsigned int size2, NumericT *group_buffer)
Definition: vector_operations.hpp:1493
viennacl::linalg::cuda::vec_element_tan_kernel
__global__ void vec_element_tan_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1438
viennacl::traits::size
vcl_size_t size(VectorType const &vec)
Generic routine for obtaining the size of a vector (ViennaCL, uBLAS, etc.)
Definition: size.hpp:239
viennacl::traits::start2
result_of::size_type< T >::type start2(T const &obj)
Definition: start.hpp:84
viennacl::linalg::cuda::inclusive_scan
void inclusive_scan(vector_base< NumericT > const &input, vector_base< NumericT > &output)
This function implements an inclusive scan using CUDA.
Definition: vector_operations.hpp:3226
viennacl::is_cpu_scalar
Helper struct for checking whether a type is a host scalar type (e.g. float, double) ...
Definition: forwards.h:448
viennacl::vector_base< NumericT >::begin
iterator begin()
Returns an iterator pointing to the beginning of the vector (STL like)
viennacl::vector_tuple
Tuple class holding pointers to multiple vectors. Mainly used as a temporary object returned from vie...
Definition: forwards.h:269
viennacl::linalg::cuda::norm_inf_impl
void norm_inf_impl(vector_base< NumericT > const &vec1, scalar< NumericT > &result)
Computes the supremum-norm of a vector.
Definition: vector_operations.hpp:2679
viennacl::linalg::cuda::inner_prod_cpu
void inner_prod_cpu(vector_base< NumericT > const &vec1, vector_base< NumericT > const &vec2, NumericT &result)
Computes the inner product of two vectors - implementation. Library users should call inner_prod(vec1...
Definition: vector_operations.hpp:1771
viennacl::traits::start
result_of::size_type< T >::type start(T const &obj)
Definition: start.hpp:44
viennacl::linalg::cuda::vec_element_atan_kernel
__global__ void vec_element_atan_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1100
viennacl::linalg::cuda::sum_cpu
void sum_cpu(vector_base< NumericT > const &vec1, NumericT &result)
Computes the maximum of a vector, first reduction stage on the GPU, second stage on the CPU...
Definition: vector_operations.hpp:2891
viennacl::linalg::cuda::scan_kernel_3
__global__ void scan_kernel_3(NumericT *Y, unsigned int startY, unsigned int incY, unsigned int sizeY, NumericT const *carries)
Definition: vector_operations.hpp:3147
viennacl::vector_base
Common base class for dense vectors, vector ranges, and vector slices.
Definition: vector_def.hpp:104
viennacl::linalg::cuda::vector_sum_kernel_floats
__global__ void vector_sum_kernel_floats(const NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, unsigned int option, NumericT *result)
Definition: vector_operations.hpp:1533
viennacl::linalg::cuda::scan_kernel_1
__global__ void scan_kernel_1(NumericT const *X, unsigned int startX, unsigned int incX, unsigned int sizeX, NumericT *Y, unsigned int startY, unsigned int incY, unsigned int scan_offset, NumericT *carries)
Definition: vector_operations.hpp:3063
viennacl::linalg::cuda::sum_impl
void sum_impl(vector_base< NumericT > const &vec1, scalar< NumericT > &result)
Computes the maximum of a vector, both reduction stages run on the GPU.
Definition: vector_operations.hpp:2874
viennacl::linalg::cuda::avbv_v
void avbv_v(vector_base< NumericT > &vec1, vector_base< NumericT > const &vec2, ScalarT1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, vector_base< NumericT > const &vec3, ScalarT2 const &beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
Definition: vector_operations.hpp:735
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::linalg::cuda::vec_element_log_kernel
__global__ void vec_element_log_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1308
viennacl::linalg::cuda::vector_sum_kernel_integers
__global__ void vector_sum_kernel_integers(const NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, unsigned int option, NumericT *result)
Definition: vector_operations.hpp:1575
viennacl::vector
Definition: forwards.h:266
viennacl::linalg::cuda::vec_element_cosh_kernel
__global__ void vec_element_cosh_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1178
viennacl::is_division
Helper metafunction for checking whether the provided type is viennacl::op_div (for division) ...
Definition: predicate.hpp:466
viennacl::linalg::cuda::av_kernel
__global__ void av_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, const NumericT *fac2, unsigned int options2, const NumericT *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:77
viennacl::linalg::cuda::vec_element_acos_kernel
__global__ void vec_element_acos_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1047
viennacl::linalg::cuda::vec_element_ceil_kernel
__global__ void vec_element_ceil_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1126
viennacl::linalg::cuda::min_cpu
void min_cpu(vector_base< NumericT > const &vec1, NumericT &result)
Computes the maximum of a vector, first reduction stage on the GPU, second stage on the CPU...
Definition: vector_operations.hpp:2847
viennacl::linalg::cuda::norm_kernel_integers
__global__ void norm_kernel_integers(const NumericT *vec, unsigned int start1, unsigned int inc1, unsigned int size1, unsigned int norm_selector, NumericT *group_buffer)
Definition: vector_operations.hpp:2331
viennacl::linalg::cuda::norm_1_impl
void norm_1_impl(vector_base< NumericT > const &vec1, scalar< NumericT > &result)
Computes the l^1-norm of a vector.
Definition: vector_operations.hpp:2587
viennacl::traits::stride2
result_of::size_type< matrix_base< NumericT > >::type stride2(matrix_base< NumericT > const &s)
Definition: stride.hpp:65
viennacl::linalg::cuda::index_norm_inf
vcl_size_t index_norm_inf(vector_base< NumericT > const &vec1)
Computes the index of the first entry that is equal to the supremum-norm in modulus.
Definition: vector_operations.hpp:2972
viennacl::linalg::cuda::element_op_kernel
__global__ void element_op_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2, NumericT const *vec3, unsigned int start3, unsigned int inc3, unsigned int op_type)
Definition: vector_operations.hpp:869
viennacl::linalg::cuda::vec_element_exp_kernel
__global__ void vec_element_exp_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1204
predicate.hpp
All the predicates used within ViennaCL. Checks for expressions to be vectors, etc.
viennacl::linalg::cuda::plane_rotation
void plane_rotation(vector_base< NumericT > &vec1, vector_base< NumericT > &vec2, NumericT alpha, NumericT beta)
Computes a plane rotation of two vectors.
Definition: vector_operations.hpp:3032
viennacl::linalg::cuda::min_impl
void min_impl(vector_base< NumericT > const &vec1, scalar< NumericT > &result)
Computes the maximum of a vector, both reduction stages run on the GPU.
Definition: vector_operations.hpp:2819
viennacl::traits::context
viennacl::context context(T const &t)
Returns an ID for the currently active memory domain of an object.
Definition: context.hpp:40
viennacl::linalg::cuda::element_op
void element_op(matrix_base< NumericT, SizeT > &A, matrix_expression< const matrix_base< NumericT, SizeT >, const matrix_base< NumericT, SizeT >, op_element_binary< OpT > > const &proxy)
Definition: matrix_operations.hpp:548
common.hpp
Common routines for CUDA execution.
viennacl::linalg::cuda::inner_prod_impl
void inner_prod_impl(vector_base< NumericT > const &vec1, vector_base< NumericT > const &vec2, ScalarT &result)
Computes the inner product of two vectors - implementation. Library users should call inner_prod(vec1...
Definition: vector_operations.hpp:1739
viennacl::scalar_vector
Represents a vector consisting of scalars 's' only, i.e. v[i] = s for all i. To be used as an initial...
Definition: vector_def.hpp:87
viennacl::linalg::cuda::vector_swap_kernel
__global__ void vector_swap_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:826
viennacl::linalg::cuda::norm_inf_cpu
void norm_inf_cpu(vector_base< NumericT > const &vec1, NumericT &result)
Computes the supremum-norm of a vector.
Definition: vector_operations.hpp:2699
viennacl::linalg::cuda::avbv_kernel
__global__ void avbv_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, const NumericT *fac2, unsigned int options2, const NumericT *vec2, unsigned int start2, unsigned int inc2, const NumericT *fac3, unsigned int options3, const NumericT *vec3, unsigned int start3, unsigned int inc3)
Definition: vector_operations.hpp:179
viennacl::vector_base::size
size_type size() const
Returns the length of the vector (cf. std::vector)
Definition: vector_def.hpp:118
viennacl::linalg::cuda::inner_prod_2_kernel
__global__ void inner_prod_2_kernel(const NumericT *x, unsigned int startx, unsigned int stridex, unsigned int sizex, const NumericT *y0, unsigned int start0, unsigned int stride0, const NumericT *y1, unsigned int start1, unsigned int stride1, NumericT *group_results)
Definition: vector_operations.hpp:1807
viennacl::linalg::cuda::vec_element_sin_kernel
__global__ void vec_element_sin_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1360
viennacl::backend::mem_handle
Main abstraction class for multiple memory domains. Represents a buffer in either main RAM...
Definition: mem_handle.hpp:89
viennacl::vector_tuple::const_at
VectorType const & const_at(vcl_size_t i) const
Definition: vector.hpp:1146
VIENNACL_CUDA_LAST_ERROR_CHECK
#define VIENNACL_CUDA_LAST_ERROR_CHECK(message)
Definition: common.hpp:30
viennacl::op_element_binary
A tag class representing element-wise binary operations (like multiplication) on vectors or matrices...
Definition: forwards.h:130
viennacl::linalg::cuda::vec_element_cos_kernel
__global__ void vec_element_cos_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1152
viennacl::linalg::cuda::vector_assign_kernel
__global__ void vector_assign_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, unsigned int internal_size1, NumericT alpha)
Definition: vector_operations.hpp:782
viennacl::backend::memory_create
void memory_create(mem_handle &handle, vcl_size_t size_in_bytes, viennacl::context const &ctx, const void *host_ptr=NULL)
Creates an array of the specified size. If the second argument is provided, the buffer is initialized...
Definition: memory.hpp:87
viennacl::linalg::cuda::vec_element_log10_kernel
__global__ void vec_element_log10_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1334
viennacl::cuda_arg
NumericT * cuda_arg(scalar< NumericT > &obj)
Convenience helper function for extracting the CUDA handle from a ViennaCL scalar. Non-const version.
Definition: common.hpp:39
viennacl::linalg::detail::min
T min(const T &lhs, const T &rhs)
Minimum.
Definition: util.hpp:45
viennacl::vector_base::internal_size
size_type internal_size() const
Returns the internal length of the vector, which is given by size() plus the extra memory due to padd...
Definition: vector_def.hpp:120
viennacl::vector_base< NumericT >::end
iterator end()
Returns an iterator pointing to the end of the vector (STL like)
viennacl::is_product
Helper metafunction for checking whether the provided type is viennacl::op_prod (for products/multipl...
Definition: predicate.hpp:436
viennacl::linalg::cuda::vec_element_floor_kernel
__global__ void vec_element_floor_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1282
viennacl::op_element_unary
A tag class representing element-wise unary operations (like sin()) on vectors or matrices...
Definition: forwards.h:134
scalar.hpp
Implementation of the ViennaCL scalar class.
viennacl::linalg::cuda::norm_2_impl
void norm_2_impl(vector_base< NumericT > const &vec1, scalar< NumericT > &result)
Computes the l^2-norm of a vector - implementation.
Definition: vector_operations.hpp:2632
viennacl::linalg::cuda::avbv
void avbv(vector_base< NumericT > &vec1, vector_base< NumericT > const &vec2, ScalarT1 const &alpha, vcl_size_t len_alpha, bool reciprocal_alpha, bool flip_sign_alpha, vector_base< NumericT > const &vec3, ScalarT2 const &beta, vcl_size_t len_beta, bool reciprocal_beta, bool flip_sign_beta)
Definition: vector_operations.hpp:433
viennacl::linalg::cuda::inner_prod_8_kernel
__global__ void inner_prod_8_kernel(const NumericT *x, unsigned int startx, unsigned int stridex, unsigned int sizex, const NumericT *y0, unsigned int start0, unsigned int stride0, const NumericT *y1, unsigned int start1, unsigned int stride1, const NumericT *y2, unsigned int start2, unsigned int stride2, const NumericT *y3, unsigned int start3, unsigned int stride3, const NumericT *y4, unsigned int start4, unsigned int stride4, const NumericT *y5, unsigned int start5, unsigned int stride5, const NumericT *y6, unsigned int start6, unsigned int stride6, const NumericT *y7, unsigned int start7, unsigned int stride7, NumericT *group_results)
Definition: vector_operations.hpp:1942
viennacl::linalg::cuda::vec_element_sinh_kernel
__global__ void vec_element_sinh_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2)
Definition: vector_operations.hpp:1386
viennacl::linalg::cuda::detail::arg_reference
viennacl::backend::mem_handle::cuda_handle_type & arg_reference(viennacl::scalar< NumericT > &s, OtherT)
Definition: common.hpp:188
viennacl::linalg::cuda::cuda_abs
__device__ NumericT cuda_abs(NumericT val)
Definition: vector_operations.hpp:2910
enable_if.hpp
Simple enable-if variant that uses the SFINAE pattern.
viennacl::linalg::min
NumericT min(std::vector< NumericT > const &v1)
Definition: maxmin.hpp:91
viennacl::linalg::cuda::element_op_int_kernel
__global__ void element_op_int_kernel(NumericT *vec1, unsigned int start1, unsigned int inc1, unsigned int size1, NumericT const *vec2, unsigned int start2, unsigned int inc2, NumericT const *vec3, unsigned int start3, unsigned int inc3, unsigned int op_type)
Definition: vector_operations.hpp:915
viennacl::fast_copy
void fast_copy(const const_vector_iterator< SCALARTYPE, ALIGNMENT > &gpu_begin, const const_vector_iterator< SCALARTYPE, ALIGNMENT > &gpu_end, CPU_ITERATOR cpu_begin)