Optimize merge_sort algorithm for largest data sizes

[SergeyKopienko]

In this PR we extends the approach from #1933 to merge_sort algorithm.

[danhoeflinger]

I imagine you can replace these branches for the vast majority of work-groups / work-items with a single branch to check if it is a "full" subproblem.

[SergeyKopienko]

1. There is no work-group / work-item separation at all, only planar sycl::range
2. This code exists before too, I simple repacked it into this class too reduce the amount of code duplications.

[danhoeflinger]

I guess we might end up with base diagonals at the {0,0} point for individual subproblems, but if the subproblems are big enough, it is not too much wasted effort perhaps...
Nothing really to do here I don't think but I want to point it out.

[SergeyKopienko]

I am tried to explain this design in commenthttps://github.com//pull/1977#discussion_r1926601240

[danhoeflinger]

Actually... should this be oneapi::dpl::__internal::__dpl_ceiling_div(__n, __chunk * __max_wg_size)-1; ?

As you mention in a comment below, we don't need to calculate {0,0}, nor do we need the final diagonal.

Think about the situation with a single subproblem, but
__max_wg_size = 1024
__chunk = 4
__n = 8191

I would think we only want a single base diagonal, a bound for the first workgroup input data. With this function we get 2.

[SergeyKopienko]

It simple doesn't make sense to think about this because all split-points which is out of merge matrix calculated as constant {n1, n2}, where n1 is the size of the first merging data and n2 is the size of the second merging data. So in any case no extra calculations here.

About your example - we strongly get sizes starting from '8192' and next sizes multiplied by 2 :
8192, 16384, 32768 and so on.

For base diagonal for the float data type we will start from 2097152:
2097152, 4194304, 8388608 and so on.

So I believe no real problem here.

[danhoeflinger]

This makes evenly divided base diagonals, but it seems like maybe we would prefer to explicitly place them at the edges of workgroups instead?

[SergeyKopienko]

Decreasing of amount of base diagonals gives us poor performance, as far as I have checked.

[SergeyKopienko]

Additional comments.

I am tried to comment this design in

oneDPL/include/oneapi/dpl/pstl/hetero/dpcpp/parallel_backend_sycl_merge_sort.h

Line 429 in d356335

__lookup_sp(const std::size_t __linear_id_in_steps_range, const nd_range_params& __nd_range_params,

There is two reasons:

1. We save series of base diagonals data in common storage for a lot of logical merges.
2. When we calculates base diagonals, we works with their indexes in storage through all the data without additional index recalculation.
3. It follows from (1) and (2) that we cannot save both dates for the first and the last base diagonals of every merge matrix, we can save on one of them.
4. As far as we know that split-point of first base diagonal always is {0, 0} I made choice to save last base diagonal data.

[danhoeflinger]

What I would suggest for now, is to take a look at the values within __get_starting_size_limit_for_large_submitter, and see if we get reasonable results for different types when using something like this:

template <typename _Tp>
constexpr std::size_t
__get_starting_size_limit_for_large_submitter()
{
    return MagicConstantCacheSize / sizeof(_Tp) ; 
}

for some fixed value of MagicConstantCacheSize . We can have this be empirically driven to start, but at least it would be consistent with itself and have a basis in something that we can understand from the hardware.

[danhoeflinger]

It seems like this branch is asking if there is a right (upper) bound base diagonal which has successfully defined.

If this is correct, can we instead replace it with a similar check when initializing __sp_right to bound it by the end of the sequence if it is "missing"?

It seems like we are throwing away the left bound information when we don't have the right bound when we should just keep the left bound and use the right bound as the end of sequence.

If I'm misunderstanding the point of this check, please explain it.

[danhoeflinger]

After further consideration I think that this check, combined with an assumption that each subproblem (other than the first) starts with a base diagonal is what allows the final steps_between_two_base_diags steps of each subproblem to work, but also does so in an inefficient way, by falling back to the search bounded by the full subproblem, rather than by base diagonals.

I see a couple issues here:
(1) I'm not sure we can safely assume that subproblems always start with base diagonals. If this assumption is part of our algorithm, we need to make that clear in the setup, and also make it guaranteed that this is always the case.
(2) It seems we could still bound the problem by the left side base diagonal, and the end of the subproblem. As we have it now, this is very inefficient, and at regular intervals throughout the sequence based on subproblems. It probably hurts our performance a lot.

If this isn't how it works then I don't understand how we handle the elements at the end of each subproblem.

[danhoeflinger]

I think __steps is miscalculated here, and should be using __nsorted * 2 instead of __rng_size to represent the subproblem instead.

Lets look at what is done in this function, and how its adjusted later:

It calculates the __base_diag_count for the subproblem (merging two __n_sorted sized sequences).
It calculates the __steps for the full problem, all subproblems combined, by using the __rng_size.
It calculates the __steps_between_two_base_diags combining data from the subproblem (__base_diag_count) and from the full problem (__steps).

Then later, outside the function there is a correction which calculates the number of subproblems, and scales two of the fields by the number of subproblems:

__nd_range_params_this.base_diag_count *= __portions;
__nd_range_params_this.steps *= __portions;

From context, it seems this function is supposed to return values for the subproblem only.

[danhoeflinger]

I'll note I think it would be better if this function returned values for the full problem rather than the subproblem, which needed no correction later, but that is just my preference.

[adamfidel]

The type of __data_area.n1 is known in this context to be _IndexT. Can we replace decltype(__data_area.n1){0} with _IndexT{0}?

[adamfidel]

What does the _w suffix in this function's name denote? Could you leave a comment here and for __serial_merge_w?

[adamfidel]

Suggested change

	// +1 doesn't required here, because we need to calculate split points for each base diagonal
	// +1 is not required here, because we need to calculate split points for each base diagonal