Return to Question

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?)

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?)

__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;
    }
}

__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;
    }
}

1. 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?)
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).
3. 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?

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?)