Yesterday I started following the [hands-on OpenCL course](https://handsonopencl.github.io/). I now got to the point where we are requested to reimplement an approximation algorithm for Pi in OpenCL (in steps, up to a vectorized implementation). I have done so and wanted to get a review for it coupled with some questions I had.

My kernel:
        
    __kernel void pi(const long steps_total, const long steps_per_kernel,
    const double step, __global double *global_sums)
    {
        const int global_id  = get_global_id(0);
        const int local_id   = get_local_id(0);
        const int local_size = get_local_size(0);

        const long vector_size = 16;
        const double16 deltas  = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15);

        const long work_dx_start     = global_id * steps_per_kernel;
        const long work_dx_end       = min(work_dx_start + steps_per_kernel, steps_total);
        const long vectorized_dx_end = work_dx_end / vector_size * vector_size;

        double work_sum = 0.0;
        for (long dx = work_dx_start; dx < vectorized_dx_end; dx += vector_size)
        {
            const double16 mid_points = (dx - 0.5 + deltas) * step;
            const double16 partial_sums = (4.0 / (1.0 + mid_points * mid_points));

            work_sum += partial_sums.s0 + partial_sums.s1 + partial_sums.s2 + partial_sums.s3 +
                        partial_sums.s4 + partial_sums.s5 + partial_sums.s6 + partial_sums.s7 +
                        partial_sums.s8 + partial_sums.s9 + partial_sums.sa + partial_sums.sb +
                        partial_sums.sc + partial_sums.sd + partial_sums.se + partial_sums.sf;
        }
        for (long dx = vectorized_dx_end; dx < work_dx_end; dx++)
        {
            const double mid_point = (dx - 0.5) * step;
            work_sum += 4.0 / (1.0 + mid_point * mid_point);
        }

        const double group_sum = work_group_reduce_add(work_sum);
        if (local_id == 0)
        {
            global_sums[global_id / local_size] = group_sum;
        }
    }

I create `max_compute_units * preferred_work_group_size` work-items, with `preferred_work_group_size` as the number of items in a group (and so each kernel executes around `steps/global_work_size` iterations). On the host side the global sums array is finally added again.

Any comments are greatly appreciated (style, optimizations,...), I also have some questions:

 1. It there a nicer way to remove the extra loop at the end in case the iterations are not a multiple of the vector size?
2. Is there a way to sum all elements of an OpenCL vector faster than this (the elements of `partial_sums`)? The best I could find was IIRC calculating the dot product (which only seems to exist for vectors up to size 4, and even then, it was slower on my GPU). Also stepwise reducing to vectors of half the side by using addition between lower and upper has no apparent effect. I guess it may not even be possible to optimize further?
3. On my GPU, using vectors of size 16 performs the best. I take it this can change depending on the device? Is there some way to try to detect/heuristically calculate this up front (statically) and dynamically load a specific kernel? (or other techniques?)