template<typename Lambda, typename... Args>
__global__
void cu_parallel_execution_kernel(size_t N, Lambda lambda, Args... args)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if( idx < N )
    {
        lambda(idx, args...);
    }
}
template<typename Lambda, typename... Args>
void cu_parallel_execution(size_t N, Lambda lambda, Args... args)
{
    cu_parallel_execution_kernel<<<(N+1023)/1024, 1024>>>(N, lambda, args...);
    CU_CHECK_AND_SYNC();
}

void normalize_vectors(vectorfield & vf)
{
    auto lambda = [=] __device__ (size_t idx, Vector3 * vf) {
        vf[idx].normalize();
    };
    int n = vf.size();
    parallel_execution<<<(n+1023)/1024, 1024>>>(n, lambda, vf.data());
}

void add_c_cross(const scalar & c, const vectorfield & a, const vectorfield & b, vectorfield & out)
{
    auto lambda = [=] __device__ (size_t idx, scalar c, const Vector3 * a, const Vector3 * b, Vector3 * out)
    {
        out[idx] += c*a[idx].cross(b[idx]);
    };
    int n = out.size();
    _parallel_execution<<<(n+1023)/1024, 1024>>>(n, lambda, c, a.data(), b.data(), out.data());
}