All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
fft_operations.hpp
Go to the documentation of this file.
#ifndef VIENNACL_LINALG_CUDA_FFT_OPERATIONS_HPP_
#define VIENNACL_LINALG_CUDA_FFT_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/matrix.hpp>
#include <viennacl/vector.hpp>
#include "viennacl/forwards.h"
#include "viennacl/scalar.hpp"
#include "viennacl/tools/tools.hpp"
#include "viennacl/linalg/cuda/common.hpp"
#include "viennacl/linalg/host_based/vector_operations.hpp"
#include "viennacl/linalg/host_based/fft_operations.hpp"
namespace viennacl
{
namespace linalg
{
namespace cuda
{
namespace detail
{
namespace fft
{
const vcl_size_t MAX_LOCAL_POINTS_NUM = 512;
inline vcl_size_t num_bits(vcl_size_t size)
{
vcl_size_t bits_datasize = 0;
vcl_size_t ds = 1;
while (ds < size)
{
ds = ds << 1;
bits_datasize++;
}
return bits_datasize;
}
inline vcl_size_t next_power_2(vcl_size_t n)
{
n = n - 1;
vcl_size_t power = 1;
while (power < sizeof(vcl_size_t) * 8)
{
n = n | (n >> power);
power *= 2;
}
return n + 1;
}
} //namespace fft
} //namespace detail
// addition
inline __host__ __device__ float2 operator+(float2 a, float2 b)
{
return make_float2(a.x + b.x, a.y + b.y);
}
// subtract
inline __host__ __device__ float2 operator-(float2 a, float2 b)
{
return make_float2(a.x - b.x, a.y - b.y);
}
// division
template<typename SCALARTYPE>
inline __device__ float2 operator/(float2 a,SCALARTYPE b)
{
return make_float2(a.x/b, a.y/b);
}
//multiplication
inline __device__ float2 operator*(float2 in1, float2 in2)
{
return make_float2(in1.x * in2.x - in1.y * in2.y, in1.x * in2.y + in1.y * in2.x);
}
// addition
inline __host__ __device__ double2 operator+(double2 a, double2 b)
{
return make_double2(a.x + b.x, a.y + b.y);
}
// subtraction
inline __host__ __device__ double2 operator-(double2 a, double2 b)
{
return make_double2(a.x - b.x, a.y - b.y);
}
// division
template<typename SCALARTYPE>
inline __host__ __device__ double2 operator/(double2 a,SCALARTYPE b)
{
return make_double2(a.x/b, a.y/b);
}
//multiplication
inline __host__ __device__ double2 operator*(double2 in1, double2 in2)
{
return make_double2(in1.x * in2.x - in1.y * in2.y, in1.x * in2.y + in1.y * in2.x);
}
inline __device__ unsigned int get_reorder_num(unsigned int v, unsigned int bit_size)
{
v = ((v >> 1) & 0x55555555) | ((v & 0x55555555) << 1);
v = ((v >> 2) & 0x33333333) | ((v & 0x33333333) << 2);
v = ((v >> 4) & 0x0F0F0F0F) | ((v & 0x0F0F0F0F) << 4);
v = ((v >> 8) & 0x00FF00FF) | ((v & 0x00FF00FF) << 8);
v = (v >> 16) | (v << 16);
v = v >> (32 - bit_size);
return v;
}
template<typename Numeric2T, typename NumericT>
__global__ void fft_direct(
const Numeric2T * input,
Numeric2T * output,
unsigned int size,
unsigned int stride,
unsigned int batch_num,
NumericT sign,
bool is_row_major)
{
const NumericT NUM_PI(3.14159265358979323846);
for (unsigned int batch_id = 0; batch_id < batch_num; batch_id++)
{
for (unsigned int k = blockIdx.x * blockDim.x + threadIdx.x; k < size; k += gridDim.x * blockDim.x)
{
Numeric2T f;
f.x = 0;
f.y = 0;
for (unsigned int n = 0; n < size; n++)
{
Numeric2T in;
if (!is_row_major)
in = input[batch_id * stride + n]; //input index here
else
in = input[n * stride + batch_id];//input index here
NumericT sn,cs;
NumericT arg = sign * 2 * NUM_PI * k / size * n;
sn = sin(arg);
cs = cos(arg);
Numeric2T ex;
ex.x = cs;
ex.y = sn;
Numeric2T tmp;
tmp.x = in.x * ex.x - in.y * ex.y;
tmp.y = in.x * ex.y + in.y * ex.x;
f = f + tmp;
}
if (!is_row_major)
output[batch_id * stride + k] = f; // output index here
else
output[k * stride + batch_id] = f;// output index here
}
}
}
template<typename NumericT, unsigned int AlignmentV>
void direct(viennacl::vector<NumericT, AlignmentV> const & in,
viennacl::vector<NumericT, AlignmentV> & out,
vcl_size_t size, vcl_size_t stride, vcl_size_t batch_num,
NumericT sign = NumericT(-1),
viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::DATA_ORDER data_order = viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
fft_direct<<<128,128>>>(reinterpret_cast<const numeric2_type *>(viennacl::cuda_arg(in)),
reinterpret_cast< numeric2_type *>(viennacl::cuda_arg(out)),
static_cast<unsigned int>(size),
static_cast<unsigned int>(stride),
static_cast<unsigned int>(batch_num),
sign,
bool(data_order != viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR));
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_direct");
}
template<typename NumericT, unsigned int AlignmentV>
void direct(viennacl::matrix<NumericT, viennacl::row_major, AlignmentV> const & in,
viennacl::matrix<NumericT, viennacl::row_major, AlignmentV> & out,
vcl_size_t size, vcl_size_t stride, vcl_size_t batch_num,
NumericT sign = NumericT(-1),
viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::DATA_ORDER data_order = viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
fft_direct<<<128,128>>>(reinterpret_cast<const numeric2_type *>(viennacl::cuda_arg(in)),
reinterpret_cast< numeric2_type *>(viennacl::cuda_arg(out)),
static_cast<unsigned int>(size),
static_cast<unsigned int>(stride),
static_cast<unsigned int>(batch_num),
sign,
bool(data_order != viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR));
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_direct");
}
template<typename NumericT>
__global__ void fft_reorder(NumericT * input,
unsigned int bit_size,
unsigned int size,
unsigned int stride,
unsigned int batch_num,
bool is_row_major)
{
unsigned int glb_id = blockIdx.x * blockDim.x + threadIdx.x;
unsigned int glb_sz = gridDim.x * blockDim.x;
for (unsigned int batch_id = 0; batch_id < batch_num; batch_id++)
{
for (unsigned int i = glb_id; i < size; i += glb_sz)
{
unsigned int v = get_reorder_num(i, bit_size);
if (i < v)
{
if (!is_row_major)
{
NumericT tmp = input[batch_id * stride + i]; // index
input[batch_id * stride + i] = input[batch_id * stride + v];//index
input[batch_id * stride + v] = tmp;//index
}
else
{
NumericT tmp = input[i * stride + batch_id];
input[i * stride + batch_id] = input[v * stride + batch_id];
input[v * stride + batch_id] = tmp;
}
}
}
}
}
/***
* This function performs reorder of input data. Indexes are sorted in bit-reversal order.
* Such reordering should be done before in-place FFT.
*/
template<typename NumericT, unsigned int AlignmentV>
void reorder(viennacl::vector<NumericT, AlignmentV> & in,
vcl_size_t size, vcl_size_t stride, vcl_size_t bits_datasize, vcl_size_t batch_num,
viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::DATA_ORDER data_order = viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
fft_reorder<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(in)),
static_cast<unsigned int>(bits_datasize),
static_cast<unsigned int>(size),
static_cast<unsigned int>(stride),
static_cast<unsigned int>(batch_num),
static_cast<bool>(data_order));
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_reorder");
}
template<typename Numeric2T, typename NumericT>
__global__ void fft_radix2_local(Numeric2T * input,
unsigned int bit_size,
unsigned int size,
unsigned int stride,
unsigned int batch_num,
NumericT sign,
bool is_row_major)
{
__shared__ Numeric2T lcl_input[1024];
unsigned int grp_id = blockIdx.x;
unsigned int grp_num = gridDim.x;
unsigned int lcl_sz = blockDim.x;
unsigned int lcl_id = threadIdx.x;
const NumericT NUM_PI(3.14159265358979323846);
for (unsigned int batch_id = grp_id; batch_id < batch_num; batch_id += grp_num)
{
for (unsigned int p = lcl_id; p < size; p += lcl_sz)
{
unsigned int v = get_reorder_num(p, bit_size);
if (!is_row_major)
lcl_input[v] = input[batch_id * stride + p];
else
lcl_input[v] = input[p * stride + batch_id];
}
__syncthreads();
//performs Cooley-Tukey FFT on local arrayfft
for (unsigned int s = 0; s < bit_size; s++)
{
unsigned int ss = 1 << s;
NumericT cs, sn;
for (unsigned int tid = lcl_id; tid < size; tid += lcl_sz)
{
unsigned int group = (tid & (ss - 1));
unsigned int pos = ((tid >> s) << (s + 1)) + group;
Numeric2T in1 = lcl_input[pos];
Numeric2T in2 = lcl_input[pos + ss];
NumericT arg = group * sign * NUM_PI / ss;
sn = sin(arg);
cs = cos(arg);
Numeric2T ex;
ex.x = cs;
ex.y = sn;
Numeric2T tmp;
tmp.x = in2.x * ex.x - in2.y * ex.y;
tmp.y = in2.x * ex.y + in2.y * ex.x;
lcl_input[pos + ss] = in1 - tmp;
lcl_input[pos] = in1 + tmp;
}
__syncthreads();
}
//copy local array back to global memory
for (unsigned int p = lcl_id; p < size; p += lcl_sz)
{
if (!is_row_major)
input[batch_id * stride + p] = lcl_input[p]; //index
else
input[p * stride + batch_id] = lcl_input[p];
}
}
}
template<typename Numeric2T, typename NumericT>
__global__ void fft_radix2(Numeric2T * input,
unsigned int s,
unsigned int bit_size,
unsigned int size,
unsigned int stride,
unsigned int batch_num,
NumericT sign,
bool is_row_major)
{
unsigned int ss = 1 << s;
unsigned int half_size = size >> 1;
NumericT cs, sn;
const NumericT NUM_PI(3.14159265358979323846);
unsigned int glb_id = blockIdx.x * blockDim.x + threadIdx.x;
unsigned int glb_sz = gridDim.x * blockDim.x;
for (unsigned int batch_id = 0; batch_id < batch_num; batch_id++)
{
for (unsigned int tid = glb_id; tid < half_size; tid += glb_sz)
{
unsigned int group = (tid & (ss - 1));
unsigned int pos = ((tid >> s) << (s + 1)) + group;
Numeric2T in1;
Numeric2T in2;
unsigned int offset;
if (!is_row_major)
{
offset = batch_id * stride + pos;
in1 = input[offset]; //index
in2 = input[offset + ss];//index
}
else
{
offset = pos * stride + batch_id;
in1 = input[offset]; //index
in2 = input[offset + ss * stride];//index
}
NumericT arg = group * sign * NUM_PI / ss;
sn = sin(arg);
cs = cos(arg);
Numeric2T ex;
ex.x = cs;
ex.y = sn;
Numeric2T tmp;
tmp.x = in2.x * ex.x - in2.y * ex.y;
tmp.y = in2.x * ex.y + in2.y * ex.x;
if (!is_row_major)
input[offset + ss] = in1 - tmp; //index
else
input[offset + ss * stride] = in1 - tmp; //index
input[offset] = in1 + tmp; //index
}
}
}
template<typename NumericT, unsigned int AlignmentV>
void radix2(viennacl::vector<NumericT, AlignmentV> & in,
vcl_size_t size, vcl_size_t stride, vcl_size_t batch_num, NumericT sign = NumericT(-1),
viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::DATA_ORDER data_order = viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
unsigned int bit_size = viennacl::linalg::cuda::detail::fft::num_bits(size);
if (size <= viennacl::linalg::cuda::detail::fft::MAX_LOCAL_POINTS_NUM)
{
fft_radix2_local<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(in)),
static_cast<unsigned int>(bit_size),
static_cast<unsigned int>(size),
static_cast<unsigned int>(stride),
static_cast<unsigned int>(batch_num),
static_cast<NumericT>(sign),
static_cast<bool>(data_order));
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_radix2_local");
}
else
{
fft_reorder<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(in)),
static_cast<unsigned int>(bit_size),
static_cast<unsigned int>(size),
static_cast<unsigned int>(stride),
static_cast<unsigned int>(batch_num),
static_cast<bool>(data_order));
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_reorder");
for (vcl_size_t step = 0; step < bit_size; step++)
{
fft_radix2<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(in)),
static_cast<unsigned int>(step),
static_cast<unsigned int>(bit_size),
static_cast<unsigned int>(size),
static_cast<unsigned int>(stride),
static_cast<unsigned int>(batch_num),
sign,
static_cast<bool>(data_order));
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_radix2");
}
}
}
template<typename NumericT, unsigned int AlignmentV>
void radix2(viennacl::matrix<NumericT, viennacl::row_major, AlignmentV>& in,
vcl_size_t size, vcl_size_t stride, vcl_size_t batch_num, NumericT sign = NumericT(-1),
viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::DATA_ORDER data_order = viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
unsigned int bit_size = viennacl::linalg::cuda::detail::fft::num_bits(size);
if (size <= viennacl::linalg::cuda::detail::fft::MAX_LOCAL_POINTS_NUM)
{
fft_radix2_local<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(in)),
static_cast<unsigned int>(bit_size),
static_cast<unsigned int>(size),
static_cast<unsigned int>(stride),
static_cast<unsigned int>(batch_num),
sign,
static_cast<bool>(data_order));
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_radix2_local");
}
else
{
fft_reorder<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(in)),
static_cast<unsigned int>(bit_size),
static_cast<unsigned int>(size),
static_cast<unsigned int>(stride),
static_cast<unsigned int>(batch_num),
static_cast<bool>(data_order));
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_reorder");
for (vcl_size_t step = 0; step < bit_size; step++)
{
fft_radix2<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(in)),
static_cast<unsigned int>(step),
static_cast<unsigned int>(bit_size),
static_cast<unsigned int>(size),
static_cast<unsigned int>(stride),
static_cast<unsigned int>(batch_num),
sign,
static_cast<bool>(data_order));
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_radix2");
}
}
}
template<typename Numeric2T, typename NumericT>
__global__ void bluestein_post(Numeric2T * Z, Numeric2T * out, unsigned int size, NumericT sign)
{
unsigned int glb_id = blockIdx.x * blockDim.x + threadIdx.x;
unsigned int glb_sz =gridDim.x * blockDim.x;
unsigned int double_size = size << 1;
NumericT sn_a, cs_a;
const NumericT NUM_PI(3.14159265358979323846);
for (unsigned int i = glb_id; i < size; i += glb_sz)
{
unsigned int rm = i * i % (double_size);
NumericT angle = (NumericT)rm / size * (-NUM_PI);
sn_a = sin(angle);
cs_a= cos(angle);
Numeric2T b_i;
b_i.x = cs_a;
b_i.y = sn_a;
out[i].x = Z[i].x * b_i.x - Z[i].y * b_i.y;
out[i].y = Z[i].x * b_i.y + Z[i].y * b_i.x;
}
}
template<typename Numeric2T, typename NumericT>
__global__ void bluestein_pre(Numeric2T * input, Numeric2T * A, Numeric2T * B,
unsigned int size, unsigned int ext_size, NumericT sign)
{
unsigned int glb_id = blockIdx.x * blockDim.x + threadIdx.x;
unsigned int glb_sz = gridDim.x * blockDim.x;
unsigned int double_size = size << 1;
NumericT sn_a, cs_a;
const NumericT NUM_PI(3.14159265358979323846);
for (unsigned int i = glb_id; i < size; i += glb_sz)
{
unsigned int rm = i * i % (double_size);
NumericT angle = (NumericT)rm / size * NUM_PI;
sn_a = sin(-angle);
cs_a= cos(-angle);
Numeric2T a_i;
a_i.x =cs_a;
a_i.y =sn_a;
Numeric2T b_i;
b_i.x =cs_a;
b_i.y =-sn_a;
A[i].x = input[i].x * a_i.x - input[i].y * a_i.y;
A[i].y = input[i].x * a_i.y + input[i].y * a_i.x;
B[i] = b_i;
// very bad instruction, to be fixed
if (i)
B[ext_size - i] = b_i;
}
}
template<typename NumericT>
__global__ void zero2(NumericT * input1, NumericT * input2, unsigned int size)
{
for (unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; i < size; i += gridDim.x * blockDim.x)
{
input1[i].x = 0;
input1[i].y = 0;
input2[i].x = 0;
input2[i].y = 0;
}
}
template<typename NumericT, unsigned int AlignmentV>
void bluestein(viennacl::vector<NumericT, AlignmentV> & in,
viennacl::vector<NumericT, AlignmentV> & out, vcl_size_t /*batch_num*/)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
vcl_size_t size = in.size() >> 1;
vcl_size_t ext_size = viennacl::linalg::cuda::detail::fft::next_power_2(2 * size - 1);
viennacl::vector<NumericT, AlignmentV> A(ext_size << 1);
viennacl::vector<NumericT, AlignmentV> B(ext_size << 1);
viennacl::vector<NumericT, AlignmentV> Z(ext_size << 1);
zero2<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(A)),
reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(B)),
static_cast<unsigned int>(ext_size));
VIENNACL_CUDA_LAST_ERROR_CHECK("zero2");
bluestein_pre<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(in)),
reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(A)),
reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(B)),
static_cast<unsigned int>(size),
static_cast<unsigned int>(ext_size),
NumericT(1));
VIENNACL_CUDA_LAST_ERROR_CHECK("bluestein_pre");
viennacl::linalg::convolve_i(A, B, Z);
bluestein_post<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(Z)),
reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(out)),
static_cast<unsigned int>(size),
NumericT(1));
VIENNACL_CUDA_LAST_ERROR_CHECK("bluestein_post");
}
template<typename NumericT>
__global__ void fft_mult_vec(const NumericT * input1,
const NumericT * input2,
NumericT * output,
unsigned int size)
{
for (unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; i < size; i += gridDim.x * blockDim.x)
{
NumericT in1 = input1[i];
NumericT in2 = input2[i];
output[i] = in1 * in2;
}
}
template<typename NumericT, unsigned int AlignmentV>
void multiply_complex(viennacl::vector<NumericT, AlignmentV> const & input1,
viennacl::vector<NumericT, AlignmentV> const & input2,
viennacl::vector<NumericT, AlignmentV> & output)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
vcl_size_t size = input1.size() / 2;
fft_mult_vec<<<128,128>>>(reinterpret_cast<const numeric2_type *>(viennacl::cuda_arg(input1)),
reinterpret_cast<const numeric2_type *>(viennacl::cuda_arg(input2)),
reinterpret_cast< numeric2_type *>(viennacl::cuda_arg(output)),
static_cast<unsigned int>(size));
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_mult_vec");
}
template<typename Numeric2T, typename NumericT>
__global__ void fft_div_vec_scalar(Numeric2T * input1, unsigned int size, NumericT factor)
{
for (unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; i < size; i += gridDim.x*blockDim.x)
input1[i] = input1[i]/factor;
}
template<typename NumericT, unsigned int AlignmentV>
void normalize(viennacl::vector<NumericT, AlignmentV> & input)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
vcl_size_t size = input.size() >> 1;
NumericT norm_factor = static_cast<NumericT>(size);
fft_div_vec_scalar<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(input)),
static_cast<unsigned int>(size),
norm_factor);
VIENNACL_CUDA_LAST_ERROR_CHECK("fft_div_vec_scalar");
}
template<typename NumericT>
__global__ void transpose(const NumericT * input,
NumericT * output,
unsigned int row_num,
unsigned int col_num)
{
unsigned int size = row_num * col_num;
for (unsigned int i =blockIdx.x * blockDim.x + threadIdx.x; i < size; i+= gridDim.x * blockDim.x)
{
unsigned int row = i / col_num;
unsigned int col = i - row*col_num;
unsigned int new_pos = col * row_num + row;
output[new_pos] = input[i];
}
}
template<typename NumericT, unsigned int AlignmentV>
void transpose(viennacl::matrix<NumericT, viennacl::row_major, AlignmentV> const & input,
viennacl::matrix<NumericT, viennacl::row_major, AlignmentV> & output)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
transpose<<<128,128>>>(reinterpret_cast<const numeric2_type *>(viennacl::cuda_arg(input)),
reinterpret_cast< numeric2_type *>(viennacl::cuda_arg(output)),
static_cast<unsigned int>(input.internal_size1()>>1),
static_cast<unsigned int>(input.internal_size2()>>1));
VIENNACL_CUDA_LAST_ERROR_CHECK("transpose");
}
template<typename NumericT>
__global__ void transpose_inplace(
NumericT * input,
unsigned int row_num,
unsigned int col_num)
{
unsigned int size = row_num * col_num;
for (unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; i < size; i+= gridDim.x * blockDim.x)
{
unsigned int row = i / col_num;
unsigned int col = i - row*col_num;
unsigned int new_pos = col * row_num + row;
if (i < new_pos)
{
NumericT val = input[i];
input[i] = input[new_pos];
input[new_pos] = val;
}
}
}
template<typename NumericT, unsigned int AlignmentV>
void transpose(viennacl::matrix<NumericT, viennacl::row_major, AlignmentV> & input)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
transpose_inplace<<<128,128>>>(reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(input)),
static_cast<unsigned int>(input.internal_size1()>>1),
static_cast<unsigned int>(input.internal_size2() >> 1));
VIENNACL_CUDA_LAST_ERROR_CHECK("transpose_inplace");
}
template<typename RealT,typename ComplexT>
__global__ void real_to_complex(const RealT * in, ComplexT * out, unsigned int size)
{
for (unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; i < size; i += gridDim.x * blockDim.x)
{
ComplexT val;
val.x = in[i];
val.y = 0;
out[i] = val;
}
}
template<typename NumericT>
void real_to_complex(viennacl::vector_base<NumericT> const & in,
viennacl::vector_base<NumericT> & out, vcl_size_t size)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
real_to_complex<<<128,128>>>(viennacl::cuda_arg(in),
reinterpret_cast<numeric2_type *>(viennacl::cuda_arg(out)),
static_cast<unsigned int>(size));
VIENNACL_CUDA_LAST_ERROR_CHECK("real_to_complex");
}
template<typename ComplexT,typename RealT>
__global__ void complex_to_real(const ComplexT * in, RealT * out, unsigned int size)
{
for (unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; i < size; i += gridDim.x * blockDim.x)
out[i] = in[i].x;
}
template<typename NumericT>
void complex_to_real(viennacl::vector_base<NumericT> const & in,
viennacl::vector_base<NumericT>& out, vcl_size_t size)
{
typedef typename viennacl::linalg::cuda::detail::type_to_type2<NumericT>::type numeric2_type;
complex_to_real<<<128,128>>>(reinterpret_cast<const numeric2_type *>(viennacl::cuda_arg(in)),
viennacl::cuda_arg(out),
static_cast<unsigned int>(size));
VIENNACL_CUDA_LAST_ERROR_CHECK("complex_to_real");
}
template<typename NumericT>
__global__ void reverse_inplace(NumericT * vec, unsigned int size)
{
for (unsigned int i = blockIdx.x * blockDim.x + threadIdx.x; i < (size >> 1); i+=gridDim.x * blockDim.x)
{
NumericT val1 = vec[i];
NumericT val2 = vec[size - i - 1];
vec[i] = val2;
vec[size - i - 1] = val1;
}
}
template<typename NumericT>
void reverse(viennacl::vector_base<NumericT>& in)
{
vcl_size_t size = in.size();
reverse_inplace<<<128,128>>>(viennacl::cuda_arg(in), static_cast<unsigned int>(size));
VIENNACL_CUDA_LAST_ERROR_CHECK("reverse_inplace");
}
} //namespace cuda
} //namespace linalg
} //namespace viennacl
#endif /* FFT_OPERATIONS_HPP_ */
viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::DATA_ORDER
DATA_ORDER
Definition: fft_operations.hpp:49
viennacl::is_row_major
Helper class for checking whether a matrix has a row-major layout.
Definition: forwards.h:484
matrix.hpp
Implementation of the dense matrix class.
viennacl::linalg::cuda::bluestein_post
__global__ void bluestein_post(Numeric2T *Z, Numeric2T *out, unsigned int size, NumericT sign)
Definition: fft_operations.hpp:538
tools.hpp
Various little tools used here and there in ViennaCL.
viennacl::linalg::cuda::real_to_complex
__global__ void real_to_complex(const RealT *in, ComplexT *out, unsigned int size)
Definition: fft_operations.hpp:782
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::fft_direct
__global__ void fft_direct(const Numeric2T *input, Numeric2T *output, unsigned int size, unsigned int stride, unsigned int batch_num, NumericT sign, bool is_row_major)
Definition: fft_operations.hpp:140
forwards.h
This file provides the forward declarations for the main types used within ViennaCL.
step
endcode *Final step
Definition: least-squares.cpp:145
viennacl::matrix
A dense matrix class.
Definition: forwards.h:375
viennacl::linalg::cuda::fft_mult_vec
__global__ void fft_mult_vec(const NumericT *input1, const NumericT *input2, NumericT *output, unsigned int size)
Definition: fft_operations.hpp:657
NumericT
float NumericT
Definition: bisect.cpp:40
viennacl::linalg::cuda::operator+
__host__ __device__ float2 operator+(float2 a, float2 b)
Definition: fft_operations.hpp:80
viennacl::linalg::cuda::detail::fft::MAX_LOCAL_POINTS_NUM
const vcl_size_t MAX_LOCAL_POINTS_NUM
Definition: fft_operations.hpp:45
viennacl::linalg::cuda::detail::fft::num_bits
vcl_size_t num_bits(vcl_size_t size)
Definition: fft_operations.hpp:47
viennacl::linalg::cuda::reverse_inplace
__global__ void reverse_inplace(NumericT *vec, unsigned int size)
Definition: fft_operations.hpp:832
viennacl::linalg::cuda::transpose
__global__ void transpose(const NumericT *input, NumericT *output, unsigned int row_num, unsigned int col_num)
Definition: fft_operations.hpp:713
viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR
Definition: fft_operations.hpp:51
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::linalg::cuda::radix2
void radix2(viennacl::vector< NumericT, AlignmentV > &in, vcl_size_t size, vcl_size_t stride, vcl_size_t batch_num, NumericT sign=NumericT(-1), viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::DATA_ORDER data_order=viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR)
Radix-2 1D algorithm for computing Fourier transformation.
Definition: fft_operations.hpp:442
viennacl::linalg::cuda::bluestein
void bluestein(viennacl::vector< NumericT, AlignmentV > &in, viennacl::vector< NumericT, AlignmentV > &out, vcl_size_t)
Bluestein's algorithm for computing Fourier transformation.
Definition: fft_operations.hpp:622
viennacl::linalg::cuda::fft_div_vec_scalar
__global__ void fft_div_vec_scalar(Numeric2T *input1, unsigned int size, NumericT factor)
Definition: fft_operations.hpp:690
viennacl::linalg::cuda::fft_reorder
__global__ void fft_reorder(NumericT *input, unsigned int bit_size, unsigned int size, unsigned int stride, unsigned int batch_num, bool is_row_major)
Definition: fft_operations.hpp:241
viennacl::linalg::cuda::operator*
__device__ float2 operator*(float2 in1, float2 in2)
Definition: fft_operations.hpp:98
viennacl::vector_base< NumericT >
viennacl::vcl_size_t
std::size_t vcl_size_t
Definition: forwards.h:75
viennacl::vector< NumericT, AlignmentV >
viennacl::linalg::cuda::operator/
__device__ float2 operator/(float2 a, SCALARTYPE b)
Definition: fft_operations.hpp:92
viennacl::linalg::cuda::normalize
void normalize(viennacl::vector< NumericT, AlignmentV > &input)
Normalize vector on with his own size.
Definition: fft_operations.hpp:700
viennacl::row
vector_expression< const matrix_base< NumericT, F >, const unsigned int, op_row > row(const matrix_base< NumericT, F > &A, unsigned int i)
Definition: matrix.hpp:910
viennacl::linalg::convolve_i
void convolve_i(viennacl::vector< SCALARTYPE, ALIGNMENT > &input1, viennacl::vector< SCALARTYPE, ALIGNMENT > &input2, viennacl::vector< SCALARTYPE, ALIGNMENT > &output)
viennacl::linalg::cuda::zero2
__global__ void zero2(NumericT *input1, NumericT *input2, unsigned int size)
Definition: fft_operations.hpp:602
viennacl::linalg::cuda::fft_radix2
__global__ void fft_radix2(Numeric2T *input, unsigned int s, unsigned int bit_size, unsigned int size, unsigned int stride, unsigned int batch_num, NumericT sign, bool is_row_major)
Definition: fft_operations.hpp:371
common.hpp
Common routines for CUDA execution.
vector.hpp
The vector type with operator-overloads and proxy classes is defined here. Linear algebra operations ...
viennacl::linalg::cuda::get_reorder_num
__device__ unsigned int get_reorder_num(unsigned int v, unsigned int bit_size)
Definition: fft_operations.hpp:128
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::bluestein_pre
__global__ void bluestein_pre(Numeric2T *input, Numeric2T *A, Numeric2T *B, unsigned int size, unsigned int ext_size, NumericT sign)
Definition: fft_operations.hpp:564
fft_operations.hpp
Implementations of Fast Furier Transformation using a plain single-threaded or OpenMP-enabled executi...
VIENNACL_CUDA_LAST_ERROR_CHECK
#define VIENNACL_CUDA_LAST_ERROR_CHECK(message)
Definition: common.hpp:30
viennacl::matrix_base< NumericT >::internal_size2
size_type internal_size2() const
Returns the internal number of columns. Usually required for launching OpenCL kernels only...
Definition: matrix_def.hpp:240
viennacl::linalg::cuda::transpose_inplace
__global__ void transpose_inplace(NumericT *input, unsigned int row_num, unsigned int col_num)
Definition: fft_operations.hpp:746
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::cuda::reverse
void reverse(viennacl::vector_base< NumericT > &in)
Reverse vector to oposite order and save it in input vector.
Definition: fft_operations.hpp:847
viennacl::matrix_base< NumericT >::internal_size1
size_type internal_size1() const
Returns the internal number of rows. Usually required for launching OpenCL kernels only...
Definition: matrix_def.hpp:238
viennacl::linalg::cuda::multiply_complex
void multiply_complex(viennacl::vector< NumericT, AlignmentV > const &input1, viennacl::vector< NumericT, AlignmentV > const &input2, viennacl::vector< NumericT, AlignmentV > &output)
Mutiply two complex vectors and store result in output.
Definition: fft_operations.hpp:674
viennacl::linalg::cuda::complex_to_real
__global__ void complex_to_real(const ComplexT *in, RealT *out, unsigned int size)
Definition: fft_operations.hpp:809
viennacl::linalg::cuda::operator-
__host__ __device__ float2 operator-(float2 a, float2 b)
Definition: fft_operations.hpp:86
viennacl::linalg::cuda::reorder
void reorder(viennacl::vector< NumericT, AlignmentV > &in, vcl_size_t size, vcl_size_t stride, vcl_size_t bits_datasize, vcl_size_t batch_num, viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::DATA_ORDER data_order=viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR)
Definition: fft_operations.hpp:282
viennacl::linalg::cuda::direct
void direct(viennacl::vector< NumericT, AlignmentV > const &in, viennacl::vector< NumericT, AlignmentV > &out, vcl_size_t size, vcl_size_t stride, vcl_size_t batch_num, NumericT sign=NumericT(-1), viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::DATA_ORDER data_order=viennacl::linalg::host_based::detail::fft::FFT_DATA_ORDER::ROW_MAJOR)
Direct 1D algorithm for computing Fourier transformation.
Definition: fft_operations.hpp:197
viennacl::linalg::cuda::detail::fft::next_power_2
vcl_size_t next_power_2(vcl_size_t n)
Definition: fft_operations.hpp:61
viennacl::linalg::cuda::fft_radix2_local
__global__ void fft_radix2_local(Numeric2T *input, unsigned int bit_size, unsigned int size, unsigned int stride, unsigned int batch_num, NumericT sign, bool is_row_major)
Definition: fft_operations.hpp:298
viennacl::linalg::cuda::detail::type_to_type2
Definition: common.hpp:178
scalar.hpp
Implementation of the ViennaCL scalar class.
fft
ScalarType fft(std::vector< ScalarType > &in, std::vector< ScalarType > &out, unsigned int, unsigned int, unsigned int batch_size)
Definition: fft_1d.cpp:719
viennacl::linalg::detail::sign
SCALARTYPE sign(SCALARTYPE val)
Definition: qr-method-common.hpp:71
vector_operations.hpp
Implementations of NMF operations using a plain single-threaded or OpenMP-enabled execution on CPU...