OpenCL approximation of Pi using vectorized implementation

Question

Yesterday I started following the hands-on OpenCL course. 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 size 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?)

Answer 1

Addressing question #2, you can move the parallel sum out of the for loop, and only do it once.

Use another double16 variable (initialized to all zeros), accumulate the sum into it during the loop body. Then, once the loop is done, do the parallel sum to initialize work_sum before the last loop that handles that last few elements.