hello.cpp

#include <vexcl/vexcl.hpp>

vex::vector<cl_float2>
hilbert(vex::vector<float> &in, std::size_t rows, std::size_t cols) {
    // No need to pass context in:
    const auto &ctx = in.queue_list();

    // Use function instead of vector. Should save global memory IO.
    VEX_FUNCTION(h, float(size_t/*prm1 = idx*/, size_t/*prm2 = rows*/),
            VEX_STRINGIZE_SOURCE(
                float v = prm1 > prm2 / 2 ? 0 : 2;
                if (prm1 == 0) v = 1;
                if (prm2 % 2 == 0 && prm1 == prm2 / 2) v = 1;
                return v;
                )
            );

#if 0
    // Make FFT work on columns:
    vex::FFT<float, cl_float2> fft(ctx, {rows, cols}, { vex::fft::forward, vex::fft::none });
    vex::FFT<cl_float2, cl_float2> ifft(ctx, {rows, cols}, { vex::fft::inverse, vex::fft::none });

    // This way result will be initialized with proper queues (context):
    vex::vector<cl_float2> result(ctx, rows * cols);

    result = ifft(fft(in) * vex::reshape(h(vex::element_index(), rows), vex::extents[rows][cols], vex::extents[0]));
#else
    // FFT works on rows hence the swapped values for cols and rows
    vex::FFT<float, cl_float2> fft(ctx, {cols, rows}, { vex::fft::none, vex::fft::forward });
    vex::FFT<cl_float2, cl_float2> ifft(ctx, {cols, rows}, { vex::fft::none, vex::fft::inverse });

    // This way result will be initialized with proper queues (context):
    vex::vector<cl_float2> result(ctx, rows * cols);

    // Lots of reshaping needed due to row major matrix when working on columns.
    // Is there a better way?
    result = vex::reshape(
      ifft(fft(vex::reshape(in, vex::extents[cols][rows], vex::extents[1][0])) *
        vex::reshape(h(vex::element_index(), rows), vex::extents[cols][rows], vex::extents[1])
      ), vex::extents[rows][cols], vex::extents[1][0]
    );
#endif

    return result;
}

int main() {
    vex::Context ctx(vex::Filter::Env && vex::Filter::Count(1));
    std::cout << ctx << std::endl;

    const size_t n = 1024;
    const size_t m = 1;

    vex::vector<float> x(ctx, n * m);
    x = 1e-3 * vex::element_index();

    vex::vector<cl_float2> y = hilbert(x, m, n);

    std::cout << y << std::endl;
}