[pstl] Implement OpenMP backend

[nadiasvertex]

Summary

Implement the PSTL backend using OpenMP.

[MikeDvorskiy]

pstl/include/pstl/internal/parallel_backend_omp.h

343

ISO C++ forbids variable length array.

553

1. What's a goal of the structure?
2. I think better to don't introduce new stuff (struct _MinKItems), we can use std::vector instead (if we really need).

733

I think a container, for example std::vector, is better causing RAII reason.

[nadiasvertex]

This corrects all known problems in the parallel partial sort implementation. The parallel merge algorithm still needs implementing.

[nadiasvertex]

Resolved all review comments.

[nadiasvertex]

This patch brings in fixes made in the oneDPL project. All unit tests in oneDPL now pass.

[nadiasvertex]

Adds the "omp" entry into the PSTL CMakeLists.txt

[nadiasvertex]

Update the linker script for parallel STL.

[nadiasvertex]

The previous patch was somehow broken. I have regenerated it.

[nadiasvertex]

Another attempt at fixing the broken patch.

[nadiasvertex]

Perform clang-format

[nadiasvertex]

Disable lint checks

[nadiasvertex]

[pstl][omp] Rework patch to pacify clang-tidy

[nadiasvertex]

[pstl][omp] Add 'omp' option to list

[nadiasvertex]

[pstl][omp] Fix clang-tidy header guard warnings

[nadiasvertex]

[pstl][omp] Fix header file movement

[MikeDvorskiy]

pstl/include/pstl/internal/omp/parallel_for_each.h

29

On this line we call increment for an iterator value (last-first) time in one thread, because it is a forward iterator..

Probably, there is an approach to avoid such overhead.. At least, I would add //TODO: like "to think how to avoid the overhead"

34

Probably, here we have overhead, because an each loop iteration increments __first from "zero".
Yes, some it is done in parallel but... In any case make sense to think whether we avoid it, if it possible. At least, I would add //TODO: like "to think how to avoid the overhead"

pstl/include/pstl/internal/omp/parallel_invoke.h

40

Nit: redundant braces.

pstl/include/pstl/internal/omp/parallel_merge.h

43

We have redundant parentheses, according to the C++ operation priority list.

48

We have redundant parentheses, according to the C++ operation priority list.

75

The "parallel merge" pattern accepts just random access iterator types. So, we can write here just difference operator, everywhere in the patch.

Probably, it makes sense to add an assert (in the begin of the function) on random "access iterator type".

91

Nit: redundant braces. (here and everywhere in the same places)

pstl/include/pstl/internal/omp/parallel_reduce.h

29

Talking about fallback to serial call (here and everywhere in the patch) - in the TBB backend (and others) we have not such logic (excepting sort and merge algorithms). A point here - a user functor/lambda may be so complicated and, probably, parallel execution is reasonable even for "short" sequences... PSTL relies that is user's responsibility - to run execution in parallel or not. Probably, in the future, if the users will complain, we will change approach, but in any case the behavior should be consistent in the all backends.

34

We have redundant parentheses (here and everywhere in the similar places), according to the C++ operation priority list.

54

According to PSTL design, the trivial pre-checks are processed on the "middle" PSTL layer - in algorithm_iml.h, numeric_impl.h, etc.

If in some places on the "middle layer" it is not true - it should be added. And remove such trivial pre-checks on a threading backend level, m.b. just add an assert instead of.

pstl/include/pstl/internal/omp/parallel_scan.h

55

Regarding the following functions: downsweep, upsweep, __split.
The implementation if ones can be shared between TBB and OpenMP backends...
For now I think is OK. It may be done with the next patch.

pstl/include/pstl/internal/omp/parallel_stable_sort.h

122

Please, see my comment above (about trivial pre-checks).

pstl/include/pstl/internal/omp/parallel_transform_reduce.h

51

Actually, we cannot use OpenMP UDR because it requires an identical value to initialize each accumulator. "__init" parameter is initial value of reduction, not an "identity".
So, the comment should be corrected.

92

Please, see my comment above (about trivial pre-checks).

pstl/include/pstl/internal/omp/parallel_transform_scan.h

28

I believe, "//TODO" should added here- about that parallel implementation will be added later, while we have fallback to a serial call.

pstl/include/pstl/internal/omp/util.h

60 ↗

(On Diff #371932)

I believe the call class __buffer can be shared between the all host backend.
But, I think, for now is OK. It may be done later, in the one of next patch.

[nadiasvertex]

This update applies changes suggested by MikeDvorskiy's review.

[nadiasvertex]

All of the serious comments were addressed. The nits were style differences I disagree with. :-) If these nits are LLVM required style I will adjust them.

pstl/include/pstl/internal/omp/parallel_invoke.h

40

I prefer the braces to make it clear what is in the block and for maintenance reasons.

pstl/include/pstl/internal/omp/parallel_merge.h

43

I prefer them.

48

I prefer them.

[rarutyun]

pstl/include/pstl/internal/omp/parallel_for.h

65

Please modify to omp single nowait here and in all other places where omp single is used within omp parallel. Currently we have two barriers from parallel and single. The inner barrier may be eliminated.

pstl/include/pstl/internal/omp/parallel_merge.h

75

I would even say we should write __xe - __xs because std::distance works with any iterator type. Since we expect random access iterator only to be there we might hide the error with the code as it is. In case of using operator- explicitly we get the compile-time error immediately if this very operator is not defined. And one more thing. Using additional abstraction layers in places where we don't need them (in this case slightly, but) increases compilation time.

pstl/include/pstl/internal/omp/parallel_stable_sort.h

68

Why do you need such a proxy object? Could you call __parallel_move_range directly?

71

If the Callable is required please mark method as const

102

__move_value is named as the type (not the alias) that decreases readability. Wouldn't it be better to rename the variable (probably __move_value_func)? The same is true for `__move_range`

119

I don't need __exec as a name, do you? I personally like the style with comments (/*__exec */) but it's up to you

pstl/include/pstl/internal/omp/parallel_transform_reduce.h

38

Don't see the reason why we need to add such an alias. In some places you use that in others - you don't

[rarutyun]

pstl/include/pstl/internal/omp/parallel_stable_sort.h

26

Why did you prefer to make the whole class template rather than making operator() with template parameters?

30

const at the end please

38

__first2 supposed to be OutputIterator semantically. Let's rename the second template parameter to OutputIterator and the last function parameter to __d_first or something. Otherwise, readability suffers.

67

This comment is applicable if the callable is required. Why do you need to be templatized? Could those template parameters be tied to operator()

pstl/include/pstl/internal/omp/util.h

60 ↗

(On Diff #371932)

I think it should not be a big deal to create the header with the __buffer class within and include it to both TBB backend and OpenMP backend. Otherwise we will commit relatively big code duplication

83 ↗

(On Diff #371932)

I am not sure if we need this variable basing on Mikhail's comments (probably, it's still useful for some algorithms) but if it remains in the code please mark it inline.

88 ↗

(On Diff #371932)

Is there any reason to use trailing return types? As I remember correctly you use those only in this file.

99 ↗

(On Diff #371932)

Why do you use std::intptr_t? it's a type to represent the pointer as the signed value. If you already check (with enable_if) that your type is integral and you have two same types (by function signature) you may case directly to Index, may not you?

104 ↗

(On Diff #371932)

chuck_policy (a least to me) is a very confusing name in the context of Parallel STL. I had an impression that use have some special execution policy or you wrap ExecutionPolicy object with this function or something like that. Cannot propose better from the top of my head. Need to think a little bit about that.

106 ↗

(On Diff #371932)

I would recommend to split those variable to make each in its line for the sake of readability.

106 ↗

(On Diff #371932)

If I don't miss something you always construct chunk_policy object with 3 arguments. Why do you need field initialization? I would prefer to leave this struct trivially_default_constructible unless we have strong reasons not to do that.

[MikeDvorskiy]

pstl/include/pstl/internal/omp/parallel_for.h

5

I think the file header (here and in the all headers) should be arranged with "pstl LLVM" file header, like this: (w/o copyright - the contributors listed in https://github.com/llvm/llvm-project/blob/main/pstl/CREDITS.txt)

===----------------------------------------------------------------------===

Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
See https://llvm.org/LICENSE.txt for license information.
SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception

===----------------------------------------------------------------------===//

[nadiasvertex]

pstl/include/pstl/internal/omp/parallel_merge.h

75

Fixed.

pstl/include/pstl/internal/omp/parallel_stable_sort.h

26

No reason. Happy to change it.

102

I'm not sure why you think adding a _func suffix makes this clearer. I'm not opposed to it, but most functors in the STL aren't suffixed with "_func" and that's not confusing.

119

I agree. This was an accident.

pstl/include/pstl/internal/omp/parallel_transform_reduce.h

38

Isn't calling omp_get_num_threads() slower than reading a variable?

51

This is fixed.

pstl/include/pstl/internal/omp/util.h

60 ↗

(On Diff #371932)

I made this change for the oneDPL commit, but it was deferred. I can make this change here if that is the consensus.

83 ↗

(On Diff #371932)

Many algorithms chunk up the work. If we don't have a variable we don't know how big to make the chunks. I'm happy to make it inline. Also, my testing shows that it should probably be 2048.

88 ↗

(On Diff #371932)

I prefer them. If the preference in PSTL is to not, I can change it.

99 ↗

(On Diff #371932)

Good point! I will fix that!

104 ↗

(On Diff #371932)

chunk_policy is intended to refer to the information needed to apply the chunking algorithm. We could call it chunk_boundaries or __chunk_partition_details.

106 ↗

(On Diff #371932)

I think that the compiler complains if I don't initialize the values. I can remove it and see. It might have been an artifact from other refactoring.

[nadiasvertex]

Add changes from review by rarutyan.

[nadiasvertex]

Commented on, and marked many reviews as fixed.

pstl/include/pstl/internal/omp/parallel_stable_sort.h

68

A callable is required.

[nadiasvertex]

Fix pstl for removal of std::result_of (use std::invoke_result)

[nadiasvertex]

Minor compilation fixes from typos.

[nadiasvertex]

Remove copyright.

[nadiasvertex]

Fix additional uses of std::distance() and remove serial fallbacks.

[nadiasvertex]

Apply the clang-format changes from CI.

[nadiasvertex]

Add come comments, minor formatting changes.

[MikeDvorskiy]

Christopher, one more comment - we discussed and agreed that a better name for macro _PSTL_PAR_BACKEND_OMP is

_PSTL_PAR_BACKEND_OPENMP

(I've counted 3 matches in the patch.)

pstl/include/pstl/internal/omp/parallel_merge.h

54

Talking about the redundant parentheses on line 45 - yes, perhaps, I can agree that it improves readability... But, here, on the line 53, I guess, it worsens readability, IMHO.

BTW,
On the line 103 the similar code written without redundant parentheses -

__parallel_merge_body(
            __mid - __xs, __xe - __mid,...

pstl/include/pstl/internal/omp/parallel_stable_sort.h

38

According to my and Ruslan's comments about usage of "std::distance" in case of random access iterators - yes, as we agreed, you have changed it by difference.... but why here (and on line 55) it was kept?

pstl/include/pstl/internal/omp/parallel_transform_reduce.h

84

The lines (83-87) look a kind of "merge artifacts"....

[nadiasvertex]

Rename the config key to _PSTL_PAR_BACKEND_OPENMP.
Remove redundant parentheses.
Fix merge artifacts.

[nadiasvertex]

Remove this _Size alias in parallel_transform_reduce

[nadiasvertex]

Addressed all known concerns with the last patch.

pstl/include/pstl/internal/omp/parallel_merge.h

54

I did a search and replace. Those parentheses were a side-effect of a mistake I made in the regular expression. Sorry. This is fixed now.

pstl/include/pstl/internal/omp/parallel_stable_sort.h

38

Yup. It wasn't clear to me that this was also required to be a random access iterator. I have made the change here too.

pstl/include/pstl/internal/omp/parallel_transform_reduce.h

38

I misunderstood what you were talking about. I have fixed this.

[nadiasvertex]

Fix typo

[nadiasvertex]

Fix libc++ compile issues.

[nadiasvertex]

Undo clang-format of pstl_config.h.

[rarutyun]

I've put some small comments but they are definitely not a showstopper for this patch to be merged. Whoever merges it please just make sure to apply those.

Other than that, looks good to me. Remaining improvements can be done later.

pstl/CREDITS.txt

22

Please add new line when merge this review

pstl/include/pstl/internal/omp/parallel_for.h

12

In all (or almost all) files the guard macro is not started from underscore. It's inconsistent with the rest of the code and makes somebody think that this macro is public.

It probably doesn't deserve one more patch but I would like to make sure it's properly merged to main branch.

[ldionne]

Thanks for the patch! Please ping me once you've applied those comments and we should be close to ready to ship this.

pstl/include/pstl/internal/omp/parallel_for.h

12

Please fix. This should be _PSTL_INTERNAL_OMP_PARALLEL_FOR_H (and similarly throughout).

15

Missing #include <cstddef>

pstl/include/pstl/internal/omp/parallel_for_each.h

28

It should be sufficient to use std::distance, not ::std::distance. Any reason why you're doing the latter? Applies in a bunch of places.

pstl/include/pstl/internal/omp/parallel_merge.h

93–95

Apply throughout.

pstl/include/pstl/internal/omp/util.h

28 ↗

(On Diff #378448)

This error message isn't too helpful. Something like this would probably be better:

#if !defined(_OPENMP)
#    error The OpenMP PSTL backend requires OpenMP, but it looks like your platform doesn't implement it.
#endif

Is the above message accurate and what you're trying to convey to the user? Also, is there any way that you'd be able to include <omp.h> yet OpenMP wouldn't be available? If not, then the #if !defined(_OPENMP) block is useless altogether.

[nadiasvertex]

Fix the header guards naming standard. Also remove ::std uses.

[nadiasvertex]

Resolved all outstanding comments.

pstl/include/pstl/internal/omp/parallel_for_each.h

28

Someone mentioned for ADL it is better to root it in the global space. I am fine with taking it out.

pstl/include/pstl/internal/omp/util.h

28 ↗

(On Diff #378448)

Probably true. This was reduced from a larger block that checked various openmp versions. I will remove it.

[ldionne]

Leaving the correctness of the backend to @MikeDvorskiy (thanks!), but apart from that this looks good to me. Let's work on enabling the tests in the CI?

[rarutyun]

In D99836#3064347, @ldionne wrote:

Leaving the correctness of the backend to @MikeDvorskiy (thanks!), but apart from that this looks good to me.

Together with Mikhail we reviewed that patch and it looks good. The tests were launched on the (almost) latest code base in oneDPL CI (just several configurations, to be honest). Anyway, no issues observed.

Do you want us to commit this patch?

Let's work on enabling the tests in the CI?

I believe it can be done separately if you don't have opposite opinion

[MikeDvorskiy]

Commit 6069a6a5
by Mikhail Dvorskiy, 10/15/2021 03:36 PM

[pstl] Initial implementation of OpenMP backend, on behalf of Christopher Nelson [e-mail]

A couple of parallel patterns still remains serial - "Parallel partial sort", and "Parallel transform scan" - there are //TODOs in the code.

libcxx/cmake/Modules/DefineLinkerScript.cmake

function(define_linker_script target)
   set(symlink "$<TARGET_LINKER_FILE:${target}>")
   set(soname "$<TARGET_SONAME_FILE_NAME:${target}>")
 
   get_target_property(interface_libs "${target}" INTERFACE_LINK_LIBRARIES)
 
   set(link_libraries)
   if (interface_libs)
     foreach(lib IN LISTS interface_libs)
-      if ("${lib}" STREQUAL "cxx-headers")
+      if ("${lib}" MATCHES "cxx-headers|ParallelSTL")
         continue()
       endif()
       # If ${lib} is not a target, we use a dummy target which we know will
       # have an OUTPUT_NAME property so that CMake doesn't fail when evaluating
       # the non-selected branch of the `IF`. It doesn't matter what it evaluates
       # to because it's not selected, but it must not cause an error.
       # See https://gitlab.kitware.com/cmake/cmake/-/issues/21045.
       set(output_name_tgt "$<IF:$<TARGET_EXISTS:${lib}>,${lib},${target}>")

pstl/CMakeLists.txt

 if (PSTL_PARALLEL_BACKEND STREQUAL "serial")
     message(STATUS "Parallel STL uses the serial backend")
     set(_PSTL_PAR_BACKEND_SERIAL ON)
 elseif (PSTL_PARALLEL_BACKEND STREQUAL "tbb")
     find_package(TBB 2018 REQUIRED tbb OPTIONAL_COMPONENTS tbbmalloc)
     message(STATUS "Parallel STL uses TBB ${TBB_VERSION} (interface version: ${TBB_INTERFACE_VERSION})")
     target_link_libraries(ParallelSTL INTERFACE TBB::tbb)
     set(_PSTL_PAR_BACKEND_TBB ON)
+elseif (PSTL_PARALLEL_BACKEND STREQUAL "omp")
+    message(STATUS "Parallel STL uses the omp backend")
+    target_compile_options(ParallelSTL INTERFACE "-fopenmp=libomp")
+    set(_PSTL_PAR_BACKEND_OMP ON)
 else()
     message(FATAL_ERROR "Requested unknown Parallel STL backend '${PSTL_PARALLEL_BACKEND}'.")
 endif()
 
 set(PSTL_GENERATED_HEADERS_DIR "${CMAKE_CURRENT_BINARY_DIR}/generated_headers")
 set(PSTL_CONFIG_SITE_PATH "${PSTL_GENERATED_HEADERS_DIR}/__pstl_config_site")
 configure_file("include/__pstl_config_site.in"
                "${PSTL_CONFIG_SITE_PATH}"

pstl/CREDITS.txt

 
 N: Intel Corporation
 W: http://www.intel.com
 D: Created the initial implementation.
 
 N: Thomas Rodgers
 E: [e-mail]
 D: Identifier name transformation for inclusion in a Standard C++ library.
+
+N: Christopher Nelson
+E: [e-mail]
+D: Add support for an OpenMP backend.
 No newline at end of file

pstl/include/__pstl_config_site.in

 //===----------------------------------------------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef __PSTL_CONFIG_SITE
 #define __PSTL_CONFIG_SITE
 
 #cmakedefine _PSTL_PAR_BACKEND_SERIAL
 #cmakedefine _PSTL_PAR_BACKEND_TBB
+#cmakedefine _PSTL_PAR_BACKEND_OMP
 #cmakedefine _PSTL_HIDE_FROM_ABI_PER_TU
 
 #endif // __PSTL_CONFIG_SITE

pstl/include/pstl/internal/omp/parallel_for.h

+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+template <class _Index, class _Fp>
+void
+__parallel_for_body(_Index __first, _Index __last, _Fp __f)
+{
+    // initial partition of the iteration space into chunks
+    auto __policy = __omp_backend::__chunk_partitioner(__first, __last);
+
+    // To avoid over-subscription we use taskloop for the nested parallelism
+    _PSTL_PRAGMA(omp taskloop untied mergeable)
+    for (std::size_t __chunk = 0; __chunk < __policy.__n_chunks; ++__chunk)
+    {
+        __omp_backend::__process_chunk(__policy, __first, __chunk, __f);
+    }
+}
+
+//------------------------------------------------------------------------
+// Notation:
+// Evaluation of brick f[i,j) for each subrange [i,j) of [first, last)
+//------------------------------------------------------------------------
+
+template <class _ExecutionPolicy, class _Index, class _Fp>
+void
+__parallel_for(_ExecutionPolicy&&, _Index __first, _Index __last, _Fp __f)
+{
+    if (__should_run_serial(__first, __last))
+    {
+        __f(__first, __last);
+        return;
+    }
+
+    if (omp_in_parallel())
+    {
+        // we don't create a nested parallel region in an existing parallel
+        // region: just create tasks
+        __pstl::__omp_backend::__parallel_for_body(__first, __last, __f);
+    }
+    else
+    {
+        // in any case (nested or non-nested) one parallel region is created and
+        // only one thread creates a set of tasks
+        _PSTL_PRAGMA(omp parallel)
+        _PSTL_PRAGMA(omp single) { __pstl::__omp_backend::__parallel_for_body(__first, __last, __f); }
+    }
+}
+
+} // namespace __omp_backend
+} // namespace __pstl

pstl/include/pstl/internal/omp/parallel_for_each.h

+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+template <class _ForwardIterator, class _Fp>
+void
+__parallel_for_each_body(_ForwardIterator __first, _ForwardIterator __last, _Fp __f)
+{
+    using DifferenceType = typename ::std::iterator_traits<_ForwardIterator>::difference_type;
+    auto __size = ::std::distance(__first, __last);
+
+    _PSTL_PRAGMA(omp taskloop untied mergeable)
+    for (DifferenceType __index = 0; __index < __size; ++__index)
+    {
+        auto __iter = ::std::next(__first, __index);
+        __f(*__iter);
+    }
+}
+
+template <class _ExecutionPolicy, class _ForwardIterator, class _Fp>
+void
+__parallel_for_each(_ExecutionPolicy&&, _ForwardIterator __first, _ForwardIterator __last, _Fp __f)
+{
+    if (omp_in_parallel())
+    {
+        // we don't create a nested parallel region in an existing parallel
+        // region: just create tasks
+        __pstl::__omp_backend::__parallel_for_each_body(__first, __last, __f);
+    }
+    else
+    {
+        // in any case (nested or non-nested) one parallel region is created and
+        // only one thread creates a set of tasks
+        _PSTL_PRAGMA(omp parallel)
+        _PSTL_PRAGMA(omp single) { __pstl::__omp_backend::__parallel_for_each_body(__first, __last, __f); }
+    }
+}
+
+} // namespace __omp_backend
+} // namespace __pstl

pstl/include/pstl/internal/omp/parallel_invoke.h

+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+template <typename _F1, typename _F2>
+void
+__parallel_invoke_body(_F1&& __f1, _F2&& __f2)
+{
+    _PSTL_PRAGMA(omp taskgroup)
+    {
+        _PSTL_PRAGMA(omp task untied mergeable) { std::forward<_F1>(__f1)(); }
+        _PSTL_PRAGMA(omp task untied mergeable) { std::forward<_F2>(__f2)(); }
+    }
+}
+
+template <class _ExecutionPolicy, typename _F1, typename _F2>
+void
+__parallel_invoke(_ExecutionPolicy&&, _F1&& __f1, _F2&& __f2)
+{
+    if (omp_in_parallel())
+    {
+        __parallel_invoke_body(std::forward<_F1>(__f1), std::forward<_F2>(__f2));
+    }
+    else
+    {
+        _PSTL_PRAGMA(omp parallel)
+        _PSTL_PRAGMA(omp single)
+        __parallel_invoke_body(std::forward<_F1>(__f1), std::forward<_F2>(__f2));
+    }
+}
+
+} // namespace __omp_backend
+} // namespace __pstl

pstl/include/pstl/internal/omp/parallel_merge.h

+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+template <typename _RandomAccessIterator1, typename _RandomAccessIterator2, typename _RandomAccessIterator3,
+          typename _Compare, typename _LeafMerge>
+void
+__parallel_merge_body(std::size_t __size_x, std::size_t __size_y, _RandomAccessIterator1 __xs,
+                      _RandomAccessIterator1 __xe, _RandomAccessIterator2 __ys, _RandomAccessIterator2 __ye,
+                      _RandomAccessIterator3 __zs, _Compare __comp, _LeafMerge __leaf_merge)
+{
+
+    if (__size_x + __size_y <= __omp_backend::__default_chunk_size)
+    {
+        __leaf_merge(__xs, __xe, __ys, __ye, __zs, __comp);
+        return;
+    }
+
+    _RandomAccessIterator1 __xm;
+    _RandomAccessIterator2 __ym;
+
+    if (__size_x < __size_y)
+    {
+        __ym = __ys + (__size_y / 2);
+        __xm = ::std::upper_bound(__xs, __xe, *__ym, __comp);
+    }
+    else
+    {
+        __xm = __xs + (__size_x / 2);
+        __ym = ::std::lower_bound(__ys, __ye, *__xm, __comp);
+    }
+
+    auto __zm = __zs + std::distance(__xs, __xm) + std::distance(__ys, __ym);
+
+    _PSTL_PRAGMA(omp task untied mergeable default(none)
+                     firstprivate(__xs, __xm, __ys, __ym, __zs, __comp, __leaf_merge))
+    __parallel_merge_body(std::distance(__xs, __xm), std::distance(__ys, __ym), __xs, __xm, __ys, __ym, __zs, __comp,
+                          __leaf_merge);
+
+    _PSTL_PRAGMA(omp task untied mergeable default(none)
+                     firstprivate(__xm, __xe, __ym, __ye, __zm, __comp, __leaf_merge))
+    __parallel_merge_body(std::distance(__xm, __xe), std::distance(__ym, __ye), __xm, __xe, __ym, __ye, __zm, __comp,
+                          __leaf_merge);
+
+    _PSTL_PRAGMA(omp taskwait)
+}
+
+template <class _ExecutionPolicy, typename _RandomAccessIterator1, typename _RandomAccessIterator2,
+          typename _RandomAccessIterator3, typename _Compare, typename _LeafMerge>
+void
+__parallel_merge(_ExecutionPolicy&& /*__exec*/, _RandomAccessIterator1 __xs, _RandomAccessIterator1 __xe,
+                 _RandomAccessIterator2 __ys, _RandomAccessIterator2 __ye, _RandomAccessIterator3 __zs, _Compare __comp,
+                 _LeafMerge __leaf_merge)
+
+{
+    std::size_t __size_x = std::distance(__xs, __xe);
+    std::size_t __size_y = std::distance(__ys, __ye);
+
+    // If the size is too small, don't bother parallelizing it.
+    if (__size_x + __size_y <= __default_chunk_size)
+    {
+        __leaf_merge(__xs, __xe, __ys, __ye, __zs, __comp);
+        return;
+    }
+
+    /*
+     * Run the merge in parallel by chunking it up. Use the smaller range (if any) as the iteration range, and the
+     * larger range as the search range.
+     */
+
+    if (omp_in_parallel())
+    {
+        __parallel_merge_body(__size_x, __size_y, __xs, __xe, __ys, __ye, __zs, __comp, __leaf_merge);
+    }
+    else
+    {
+        _PSTL_PRAGMA(omp parallel)
+        {
+            _PSTL_PRAGMA(omp single)
+            __parallel_merge_body(__size_x, __size_y, __xs, __xe, __ys, __ye, __zs, __comp, __leaf_merge);
+        }
+    }
+}
+
+} // namespace __omp_backend
+} // namespace __pstl

pstl/include/pstl/internal/omp/parallel_reduce.h

+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+template <class _RandomAccessIterator, class _Value, typename _RealBody, typename _Reduction>
+_Value
+__parallel_reduce_body(_RandomAccessIterator __first, _RandomAccessIterator __last, _Value __identity,
+                       _RealBody __real_body, _Reduction __reduce)
+{
+    if (__should_run_serial(__first, __last))
+    {
+        return __real_body(__first, __last, __identity);
+    }
+
+    auto __middle = __first + ((__last - __first) / 2);
+    _Value __v1(__identity), __v2(__identity);
+    __parallel_invoke_body(
+        [&]() { __v1 = __parallel_reduce_body(__first, __middle, __identity, __real_body, __reduce); },
+        [&]() { __v2 = __parallel_reduce_body(__middle, __last, __identity, __real_body, __reduce); });
+
+    return __reduce(__v1, __v2);
+}
+
+//------------------------------------------------------------------------
+// Notation:
+//      r(i,j,init) returns reduction of init with reduction over [i,j)
+//      c(x,y) combines values x and y that were the result of r
+//------------------------------------------------------------------------
+
+template <class _ExecutionPolicy, class _RandomAccessIterator, class _Value, typename _RealBody, typename _Reduction>
+_Value
+__parallel_reduce(_ExecutionPolicy&&, _RandomAccessIterator __first, _RandomAccessIterator __last, _Value __identity,
+                  _RealBody __real_body, _Reduction __reduction)
+{
+    if (__first == __last)
+        return __identity;
+
+    // Don't bother parallelizing the work if the size is too small.
+    if (__should_run_serial(__first, __last))
+    {
+        return __real_body(__first, __last, __identity);
+    }
+
+    // We don't create a nested parallel region in an existing parallel region:
+    // just create tasks.
+    if (omp_in_parallel())
+    {
+        return __pstl::__omp_backend::__parallel_reduce_body(__first, __last, __identity, __real_body, __reduction);
+    }
+
+    // In any case (nested or non-nested) one parallel region is created and only
+    // one thread creates a set of tasks.
+    _Value __res = __identity;
+
+    _PSTL_PRAGMA(omp parallel)
+    _PSTL_PRAGMA(omp single)
+    {
+        __res = __pstl::__omp_backend::__parallel_reduce_body(__first, __last, __identity, __real_body, __reduction);
+    }
+
+    return __res;
+}
+
+} // namespace __omp_backend
+} // namespace __pstl

pstl/include/pstl/internal/omp/parallel_scan.h

+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+template <typename _Index>
+_Index
+__split(_Index __m)
+{
+    _Index __k = 1;
+    while (2 * __k < __m)
+        __k *= 2;
+    return __k;
+}
+
+template <typename _Index, typename _Tp, typename _Rp, typename _Cp>
+void
+__upsweep(_Index __i, _Index __m, _Index __tilesize, _Tp* __r, _Index __lastsize, _Rp __reduce, _Cp __combine)
+{
+    if (__m == 1)
+        __r[0] = __reduce(__i * __tilesize, __lastsize);
+    else
+    {
+        _Index __k = __split(__m);
+        __omp_backend::__parallel_invoke_body(
+            [=] { __omp_backend::__upsweep(__i, __k, __tilesize, __r, __tilesize, __reduce, __combine); },
+            [=] {
+                __omp_backend::__upsweep(__i + __k, __m - __k, __tilesize, __r + __k, __lastsize, __reduce, __combine);
+            });
+        if (__m == 2 * __k)
+            __r[__m - 1] = __combine(__r[__k - 1], __r[__m - 1]);
+    }
+}
+
+template <typename _Index, typename _Tp, typename _Cp, typename _Sp>
+void
+__downsweep(_Index __i, _Index __m, _Index __tilesize, _Tp* __r, _Index __lastsize, _Tp __initial, _Cp __combine,
+            _Sp __scan)
+{
+    if (__m == 1)
+        __scan(__i * __tilesize, __lastsize, __initial);
+    else
+    {
+        const _Index __k = __split(__m);
+        __omp_backend::__parallel_invoke_body(
+            [=] { __omp_backend::__downsweep(__i, __k, __tilesize, __r, __tilesize, __initial, __combine, __scan); },
+            // Assumes that __combine never throws.
+            // TODO: Consider adding a requirement for user functors to be constant.
+            [=, &__combine]
+            {
+                __omp_backend::__downsweep(__i + __k, __m - __k, __tilesize, __r + __k, __lastsize,
+                                           __combine(__initial, __r[__k - 1]), __combine, __scan);
+            });
+    }
+}
+
+template <typename _ExecutionPolicy, typename _Index, typename _Tp, typename _Rp, typename _Cp, typename _Sp,
+          typename _Ap>
+void
+__parallel_strict_scan_body(_Index __n, _Tp __initial, _Rp __reduce, _Cp __combine, _Sp __scan, _Ap __apex)
+{
+    _Index __p = omp_get_num_threads();
+    const _Index __slack = 4;
+    _Index __tilesize = (__n - 1) / (__slack * __p) + 1;
+    _Index __m = (__n - 1) / __tilesize;
+    __buffer<_Tp> __buf(__m + 1);
+    _Tp* __r = __buf.get();
+
+    __omp_backend::__upsweep(_Index(0), _Index(__m + 1), __tilesize, __r, __n - __m * __tilesize, __reduce, __combine);
+
+    std::size_t __k = __m + 1;
+    _Tp __t = __r[__k - 1];
+    while ((__k &= __k - 1))
+    {
+        __t = __combine(__r[__k - 1], __t);
+    }
+
+    __apex(__combine(__initial, __t));
+    __omp_backend::__downsweep(_Index(0), _Index(__m + 1), __tilesize, __r, __n - __m * __tilesize, __initial,
+                               __combine, __scan);
+}
+
+template <class _ExecutionPolicy, typename _Index, typename _Tp, typename _Rp, typename _Cp, typename _Sp, typename _Ap>
+void
+__parallel_strict_scan(_ExecutionPolicy&&, _Index __n, _Tp __initial, _Rp __reduce, _Cp __combine, _Sp __scan,
+                       _Ap __apex)
+{
+    if (__n <= __default_chunk_size)
+    {
+        _Tp __sum = __initial;
+        if (__n)
+        {
+            __sum = __combine(__sum, __reduce(_Index(0), __n));
+        }
+        __apex(__sum);
+        if (__n)
+        {
+            __scan(_Index(0), __n, __initial);
+        }
+        return;
+    }
+
+    if (omp_in_parallel())
+    {
+        __pstl::__omp_backend::__parallel_strict_scan_body<_ExecutionPolicy>(__n, __initial, __reduce, __combine,
+                                                                             __scan, __apex);
+    }
+    else
+    {
+        _PSTL_PRAGMA(omp parallel)
+        _PSTL_PRAGMA(omp single)
+        {
+            __pstl::__omp_backend::__parallel_strict_scan_body<_ExecutionPolicy>(__n, __initial, __reduce, __combine,
+                                                                                 __scan, __apex);
+        }
+    }
+}
+
+} // namespace __omp_backend
+} // namespace __pstl

pstl/include/pstl/internal/omp/parallel_stable_partial_sort.h

+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+template <typename _RandomAccessIterator, typename _Compare, typename _LeafSort>
+void
+__parallel_stable_partial_sort(_RandomAccessIterator __xs, _RandomAccessIterator __xe, _Compare __comp,
+                               _LeafSort __leaf_sort, std::size_t __nsort)
+{
+    //TODO: "Parallel partial sort needs to be implemented.");
+    __leaf_sort(__xs, __xe, __comp);
+}
+
+} // namespace __omp_backend
+} // namespace __pstl

pstl/include/pstl/internal/omp/parallel_stable_sort.h

+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+namespace __sort_details
+{
+template <typename Iterator1, typename Iterator2>
+struct __move_value
+{
+    void
+    operator()(Iterator1 __x, Iterator2 __z)
+    {
+        *__z = ::std::move(*__x);
+    }
+};
+
+template <typename Iterator1, typename Iterator2>
+Iterator2
+__parallel_move_range(Iterator1 __first1, Iterator1 __last1, Iterator2 __first2)
+{
+    std::size_t __size = std::distance(__first1, __last1);
+
+    // Perform serial moving of small chunks
+
+    if (__size <= __default_chunk_size)
+    {
+        return ::std::move(__first1, __last1, __first2);
+    }
+
+    // Perform parallel moving of larger chunks
+    auto __policy = __omp_backend::__chunk_partitioner(__first1, __last1);
+
+    _PSTL_PRAGMA(omp taskloop)
+    for (std::size_t __chunk = 0; __chunk < __policy.__n_chunks; ++__chunk)
+    {
+        __omp_backend::__process_chunk(__policy, __first1, __chunk,
+                                       [&](auto __chunk_first, auto __chunk_last)
+                                       {
+                                           auto __chunk_offset = ::std::distance(__first1, __chunk_first);
+                                           auto __output = __first2 + __chunk_offset;
+                                           ::std::move(__chunk_first, __chunk_last, __output);
+                                       });
+    }
+
+    return __first2 + __size;
+}
+
+template <typename Iterator1, typename Iterator2>
+struct __move_range
+{
+    Iterator2
+    operator()(Iterator1 __first1, Iterator1 __last1, Iterator2 __first2)
+    {
+        return __parallel_move_range(__first1, __last1, __first2);
+    }
+};
+} // namespace __sort_details
+
+template <typename _RandomAccessIterator, typename _Compare, typename _LeafSort>
+void
+__parallel_stable_sort_body(_RandomAccessIterator __xs, _RandomAccessIterator __xe, _Compare __comp,
+                            _LeafSort __leaf_sort)
+{
+    using _ValueType = typename std::iterator_traits<_RandomAccessIterator>::value_type;
+    using _VecType = typename std::vector<_ValueType>;
+    using _OutputIterator = typename _VecType::iterator;
+    using _MoveValueType = typename __omp_backend::__sort_details::__move_value<_RandomAccessIterator, _OutputIterator>;
+    using _MoveRangeType = __omp_backend::__sort_details::__move_range<_RandomAccessIterator, _OutputIterator>;
+
+    if (__should_run_serial(__xs, __xe))
+    {
+        __leaf_sort(__xs, __xe, __comp);
+    }
+    else
+    {
+        std::size_t __size = std::distance(__xs, __xe);
+        auto __mid = __xs + (__size / 2);
+        __parallel_invoke_body([&]() { __parallel_stable_sort_body(__xs, __mid, __comp, __leaf_sort); },
+                               [&]() { __parallel_stable_sort_body(__mid, __xe, __comp, __leaf_sort); });
+
+        // Perform a parallel merge of the sorted ranges into __output.
+        _VecType __output(__size);
+        _MoveValueType __move_value;
+        _MoveRangeType __move_range;
+        __utils::__serial_move_merge __merge(__size);
+        __parallel_merge_body(
+            std::distance(__xs, __mid), std::distance(__mid, __xe), __xs, __mid, __mid, __xe, __output.begin(), __comp,
+            [&__merge, &__move_value, &__move_range](_RandomAccessIterator __as, _RandomAccessIterator __ae,
+                                                     _RandomAccessIterator __bs, _RandomAccessIterator __be,
+                                                     _OutputIterator __cs, _Compare __comp)
+            { __merge(__as, __ae, __bs, __be, __cs, __comp, __move_value, __move_value, __move_range, __move_range); });
+
+        // Move the values from __output back in the original source range.
+        __sort_details::__parallel_move_range(__output.begin(), __output.end(), __xs);
+    }
+}
+
+template <class _ExecutionPolicy, typename _RandomAccessIterator, typename _Compare, typename _LeafSort>
+void
+__parallel_stable_sort(_ExecutionPolicy&& __exec, _RandomAccessIterator __xs, _RandomAccessIterator __xe,
+                       _Compare __comp, _LeafSort __leaf_sort, std::size_t __nsort = 0)
+{
+    if (__xs >= __xe)
+    {
+        return;
+    }
+
+    auto __count = static_cast<std::size_t>(std::distance(__xs, __xe));
+    if (__count <= __default_chunk_size || __nsort < __count)
+    {
+        __leaf_sort(__xs, __xe, __comp);
+        return;
+    }
+
+    // TODO: the partial sort implementation should
+    // be shared with the other backends.
+
+    if (omp_in_parallel())
+    {
+        if (__count <= __nsort)
+        {
+            __parallel_stable_sort_body(__xs, __xe, __comp, __leaf_sort);
+        }
+        else
+        {
+            __parallel_stable_partial_sort(__xs, __xe, __comp, __leaf_sort, __nsort);
+        }
+    }
+    else
+    {
+        _PSTL_PRAGMA(omp parallel)
+        _PSTL_PRAGMA(omp single)
+        if (__count <= __nsort)
+        {
+            __parallel_stable_sort_body(__xs, __xe, __comp, __leaf_sort);
+        }
+        else
+        {
+            __parallel_stable_partial_sort(__xs, __xe, __comp, __leaf_sort, __nsort);
+        }
+    }
+}
+
+} // namespace __omp_backend
+} // namespace __pstl

pstl/include/pstl/internal/omp/parallel_transform_reduce.h

+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+//------------------------------------------------------------------------
+// parallel_transform_reduce
+//
+// Notation:
+//      r(i,j,init) returns reduction of init with reduction over [i,j)
+//      u(i) returns f(i,i+1,identity) for a hypothetical left identity element
+//      of r c(x,y) combines values x and y that were the result of r or u
+//------------------------------------------------------------------------
+
+template <class _RandomAccessIterator, class _UnaryOp, class _Value, class _Combiner, class _Reduction>
+auto
+__transform_reduce_body(_RandomAccessIterator __first, _RandomAccessIterator __last, _UnaryOp __unary_op, _Value __init,
+                        _Combiner __combiner, _Reduction __reduction)
+{
+    using _Size = std::size_t;
+    const _Size __num_threads = omp_get_num_threads();
+    const _Size __size = __last - __first;
+
+    // Initial partition of the iteration space into chunks. If the range is too small,
+    // this will result in a nonsense policy, so we check on the size as well below.
+    auto __policy = __omp_backend::__chunk_partitioner(__first + __num_threads, __last);
+
+    if (__size <= __num_threads || __policy.__n_chunks < 2)
+    {
+        return __reduction(__first, __last, __init);
+    }
+
+    // Here, we cannot use OpenMP UDR because we must store the init value in
+    // the combiner and it will be used several times. Although there should be
+    // the only one; we manually generate the identity elements for each thread.
+    ::std::vector<_Value> __accums;
+    __accums.reserve(__num_threads);
+
+    // initialize accumulators for all threads
+    for (_Size __i = 0; __i < __num_threads; ++__i)
+    {
+        __accums.emplace_back(__unary_op(__first + __i));
+    }
+
+    // main loop
+    _PSTL_PRAGMA(omp taskloop shared(__accums))
+    for (std::size_t __chunk = 0; __chunk < __policy.__n_chunks; ++__chunk)
+    {
+        __omp_backend::__process_chunk(__policy, __first + __num_threads, __chunk,
+                                       [&](auto __chunk_first, auto __chunk_last)
+                                       {
+                                           auto __thread_num = omp_get_thread_num();
+                                           __accums[__thread_num] =
+                                               __reduction(__chunk_first, __chunk_last, __accums[__thread_num]);
+                                       });
+    }
+
+    // combine by accumulators
+    for (_Size __i = 0; __i < __num_threads; ++__i)
+    {
+        __init = __combiner(__init, __accums[__i]);
+    }
+
+    return __init;
+}
+
+template <class _ExecutionPolicy, class _RandomAccessIterator, class _UnaryOp, class _Value, class _Combiner,
+          class _Reduction>
+_Value
+__parallel_transform_reduce(_ExecutionPolicy&&, _RandomAccessIterator __first, _RandomAccessIterator __last,
+                            _UnaryOp __unary_op, _Value __init, _Combiner __combiner, _Reduction __reduction)
+{
+
+    if (__first == __last)
+    {
+        return __init;
+    }
+
+    _Value __result = __init;
+    if (omp_in_parallel())
+    {
+        // We don't create a nested parallel region in an existing parallel
+        // region: just create tasks
+        __result = __pstl::__omp_backend::__transform_reduce_body(__first, __last, __unary_op, __init, __combiner,
+                                                                  __reduction);
+    }
+    else
+    {
+        // Create a parallel region, and a single thread will create tasks
+        // for the region.
+        _PSTL_PRAGMA(omp parallel)
+        _PSTL_PRAGMA(omp single)
+        {
+            __result = __pstl::__omp_backend::__transform_reduce_body(__first, __last, __unary_op, __init, __combiner,
+                                                                      __reduction);
+        }
+    }
+
+    return __result;
+}
+
+} // namespace __omp_backend
+} // namespace __pstl

pstl/include/pstl/internal/omp/parallel_transform_scan.h

+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+template <class _ExecutionPolicy, class _Index, class _Up, class _Tp, class _Cp, class _Rp, class _Sp>
+_Tp
+__parallel_transform_scan(_ExecutionPolicy&&, _Index __n, _Up __u, _Tp __init, _Cp __combine, _Rp __brick_reduce,
+                          _Sp __scan)
+{
+    return __scan(_Index(0), __n, __init);
+}
+
+} // namespace __omp_backend
+} // namespace __pstl

pstl/include/pstl/internal/parallel_backend.h

 namespace __par_backend = __serial_backend;
 }
 #elif defined(_PSTL_PAR_BACKEND_TBB)
 #    include "parallel_backend_tbb.h"
 namespace __pstl
 {
 namespace __par_backend = __tbb_backend;
 }
+#elif defined(_PSTL_PAR_BACKEND_OMP)
+#    include "parallel_backend_omp.h"
+namespace __pstl
+{
+namespace __par_backend = __omp_backend;
+}
 #else
 _PSTL_PRAGMA_MESSAGE("Parallel backend was not specified");
 #endif
 
 #endif /* _PSTL_PARALLEL_BACKEND_H */

pstl/include/pstl/internal/parallel_backend_omp.h

+// -*- C++ -*-
+// -*-===----------------------------------------------------------------------===//
+//
+// Copyright (C) Christopher Nelson
+//
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef _PSTL_PARALLEL_BACKEND_OMP_H
+#define _PSTL_PARALLEL_BACKEND_OMP_H
+
+#include <algorithm>
+#include <atomic>
+#include <iterator>
+#include <cstddef>
+#include <cstdio>
+#include <memory>
+#include <vector>
+#include <omp.h>
+
+#include "parallel_backend_utils.h"
+#include "unseq_backend_simd.h"
+#include "utils.h"
+
+#if !defined(_OPENMP)
+#    error _OPENMP not defined.
+#endif
+
+// Portability "#pragma" definition
+#ifdef _MSC_VER
+#    define _PSTL_PRAGMA(x) __pragma(x)
+#else
+#    define _PSTL_PRAGMA(x) _Pragma(#    x)
+#endif
+
+_PSTL_HIDE_FROM_ABI_PUSH
+
+namespace __pstl
+{
+namespace __omp_backend
+{
+
+//------------------------------------------------------------------------
+// use to cancel execution
+//------------------------------------------------------------------------
+inline void
+__cancel_execution()
+{
+}
+
+//------------------------------------------------------------------------
+// raw buffer
+//------------------------------------------------------------------------
+
+template <typename _Tp>
+class __buffer
+{
+    std::allocator<_Tp> __allocator_;
+    _Tp* __ptr_;
+    const std::size_t __buf_size_;
+    __buffer(const __buffer&) = delete;
+    void
+    operator=(const __buffer&) = delete;
+
+  public:
+    __buffer(std::size_t __n) : __allocator_(), __ptr_(__allocator_.allocate(__n)), __buf_size_(__n) {}
+
+    operator bool() const { return __ptr_ != nullptr; }
+
+    _Tp*
+    get() const
+    {
+        return __ptr_;
+    }
+    ~__buffer() { __allocator_.deallocate(__ptr_, __buf_size_); }
+};
+
+// Preliminary size of each chunk: requires further discussion
+constexpr std::size_t __default_chunk_size = 512;
+
+// Convenience function to determine when we should run serial.
+template <typename _Iterator, ::std::enable_if_t<!::std::is_integral<_Iterator>::value, bool> = true>
+constexpr auto
+__should_run_serial(_Iterator __first, _Iterator __last) -> bool
+{
+    using _difference_type = typename ::std::iterator_traits<_Iterator>::difference_type;
+    auto __size = ::std::distance(__first, __last);
+    return __size <= static_cast<_difference_type>(__default_chunk_size);
+}
+
+template <typename _Index, ::std::enable_if_t<::std::is_integral<_Index>::value, bool> = true>
+constexpr auto
+__should_run_serial(_Index __first, _Index __last) -> bool
+{
+    using _difference_type = ::std::intptr_t;
+    auto __size = __last - __first;
+    return __size <= static_cast<_difference_type>(__default_chunk_size);
+}
+
+struct __chunk_policy
+{
+    std::size_t __n_chunks{0}, __chunk_size{0}, __first_chunk_size{0};
+};
+
+// The iteration space partitioner according to __requested_chunk_size
+template <class _RandomAccessIterator, class _Size = std::size_t>
+auto
+__chunk_partitioner(_RandomAccessIterator __first, _RandomAccessIterator __last,
+                    _Size __requested_chunk_size = __default_chunk_size) -> __chunk_policy
+{
+    /*
+     * This algorithm improves distribution of elements in chunks by avoiding
+     * small tail chunks. The leftover elements that do not fit neatly into
+     * the chunk size are redistributed to early chunks. This improves
+     * utilization of the processor's prefetch and reduces the number of
+     * tasks needed by 1.
+     */
+
+    const _Size __n = __last - __first;
+    _Size __n_chunks = 0;
+    _Size __chunk_size = 0;
+    _Size __first_chunk_size = 0;
+    if (__n < __requested_chunk_size)
+    {
+        __chunk_size = __n;
+        __first_chunk_size = __n;
+        __n_chunks = 1;
+        return __chunk_policy{__n_chunks, __chunk_size, __first_chunk_size};
+    }
+
+    __n_chunks = (__n / __requested_chunk_size) + 1;
+    __chunk_size = __n / __n_chunks;
+    __first_chunk_size = __chunk_size;
+    const _Size __n_leftover_items = __n - (__n_chunks * __chunk_size);
+
+    if (__n_leftover_items == __chunk_size)
+    {
+        __n_chunks += 1;
+        return __chunk_policy{__n_chunks, __chunk_size, __first_chunk_size};
+    }
+    else if (__n_leftover_items == 0)
+    {
+        __first_chunk_size = __chunk_size;
+        return __chunk_policy{__n_chunks, __chunk_size, __first_chunk_size};
+    }
+
+    const _Size __n_extra_items_per_chunk = __n_leftover_items / __n_chunks;
+    const _Size __n_final_leftover_items = __n_leftover_items - (__n_extra_items_per_chunk * __n_chunks);
+
+    __chunk_size += __n_extra_items_per_chunk;
+    __first_chunk_size = __chunk_size + __n_final_leftover_items;
+
+    return __chunk_policy{__n_chunks, __chunk_size, __first_chunk_size};
+}
+
+template <typename _Iterator, typename _Index, typename _F>
+void
+__process_chunk(const __chunk_policy& __policy, _Iterator __base, _Index __chunk_index, _F __f)
+{
+    auto __this_chunk_size = __chunk_index == 0 ? __policy.__first_chunk_size : __policy.__chunk_size;
+    auto __index = __chunk_index == 0 ? 0
+                                      : (__chunk_index * __policy.__chunk_size) +
+                                            (__policy.__first_chunk_size - __policy.__chunk_size);
+    auto __first = __base + __index;
+    auto __last = __first + __this_chunk_size;
+    __f(__first, __last);
+}
+
+} // namespace __omp_backend
+} // namespace __pstl
+
+//------------------------------------------------------------------------
+// parallel_invoke
+//------------------------------------------------------------------------
+
+#include "./omp/parallel_invoke.h"
+
+//------------------------------------------------------------------------
+// parallel_for
+//------------------------------------------------------------------------
+
+#include "./omp/parallel_for.h"
+
+//------------------------------------------------------------------------
+// parallel_for_each
+//------------------------------------------------------------------------
+
+#include "./omp/parallel_for_each.h"
+
+//------------------------------------------------------------------------
+// parallel_reduce
+//------------------------------------------------------------------------
+
+#include "./omp/parallel_reduce.h"
+#include "./omp/parallel_transform_reduce.h"
+
+//------------------------------------------------------------------------
+// parallel_scan
+//------------------------------------------------------------------------
+
+#include "./omp/parallel_scan.h"
+#include "./omp/parallel_transform_scan.h"
+
+//------------------------------------------------------------------------
+// parallel_stable_sort
+//------------------------------------------------------------------------
+
+#include "./omp/parallel_stable_partial_sort.h"
+#include "./omp/parallel_stable_sort.h"
+
+//------------------------------------------------------------------------
+// parallel_merge
+//------------------------------------------------------------------------
+#include "./omp/parallel_merge.h"
+
+#endif //_PSTL_PARALLEL_BACKEND_OMP_H

pstl/include/pstl/internal/pstl_config.h

 // -*- C++ -*-
 //===----------------------------------------------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef _PSTL_CONFIG_H
 #define _PSTL_CONFIG_H
 
 #include <__pstl_config_site>
 
-// The version is XYYZ, where X is major, YY is minor, and Z is patch (i.e. X.YY.Z)
+// The version is XYYZ, where X is major, YY is minor, and Z is patch (i.e.
+// X.YY.Z)
 #define _PSTL_VERSION 14000
 #define _PSTL_VERSION_MAJOR (_PSTL_VERSION / 1000)
 #define _PSTL_VERSION_MINOR ((_PSTL_VERSION % 1000) / 10)
 #define _PSTL_VERSION_PATCH (_PSTL_VERSION % 10)
 
-#if !defined(_PSTL_PAR_BACKEND_SERIAL) && !defined(_PSTL_PAR_BACKEND_TBB)
-#    error "A parallel backend must be specified"
+#if !defined(_PSTL_PAR_BACKEND_SERIAL) && !defined(_PSTL_PAR_BACKEND_TBB) &&   \
+    !defined(_PSTL_PAR_BACKEND_OMP)
+#error "A parallel backend must be specified"
 #endif
 
 // Check the user-defined macro for warnings
 #if defined(PSTL_USAGE_WARNINGS)
-#    define _PSTL_USAGE_WARNINGS
+#define _PSTL_USAGE_WARNINGS
 #endif
 
 #if !defined(_PSTL_ASSERT)
-#    include <cassert>
-#    define _PSTL_ASSERT(pred) (assert((pred)))
+#include <cassert>
+#define _PSTL_ASSERT(pred) (assert((pred)))
 #endif
 
 // Portability "#pragma" definition
 #ifdef _MSC_VER
-#    define _PSTL_PRAGMA(x) __pragma(x)
+#define _PSTL_PRAGMA(x) __pragma(x)
 #else
-#    define _PSTL_PRAGMA(x) _Pragma(#    x)
+#define _PSTL_PRAGMA(x) _Pragma(#x)
 #endif
 
 #define _PSTL_STRING_AUX(x) #x
 #define _PSTL_STRING(x) _PSTL_STRING_AUX(x)
 #define _PSTL_STRING_CONCAT(x, y) x #y
 
 #ifdef _PSTL_HIDE_FROM_ABI_PER_TU
-#    define _PSTL_HIDE_FROM_ABI_PUSH                                                                                   \
-        _Pragma("clang attribute push(__attribute__((internal_linkage)), apply_to=any(function,record))")
-#    define _PSTL_HIDE_FROM_ABI_POP _Pragma("clang attribute pop")
+#define _PSTL_HIDE_FROM_ABI_PUSH                                               \
+  _Pragma("clang attribute push(__attribute__((internal_linkage)), "           \
+          "apply_to=any(function,record))")
+#define _PSTL_HIDE_FROM_ABI_POP _Pragma("clang attribute pop")
 #else
-#    define _PSTL_HIDE_FROM_ABI_PUSH /* nothing */
-#    define _PSTL_HIDE_FROM_ABI_POP  /* nothing */
+#define _PSTL_HIDE_FROM_ABI_PUSH /* nothing */
+#define _PSTL_HIDE_FROM_ABI_POP  /* nothing */
 #endif
 
-// note that when ICC or Clang is in use, _PSTL_GCC_VERSION might not fully match
-// the actual GCC version on the system.
-#define _PSTL_GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
+// note that when ICC or Clang is in use, _PSTL_GCC_VERSION might not fully
+// match the actual GCC version on the system.
+#define _PSTL_GCC_VERSION                                                      \
+  (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__)
 
 #if defined(__clang__)
 // according to clang documentation, version can be vendor specific
-#    define _PSTL_CLANG_VERSION (__clang_major__ * 10000 + __clang_minor__ * 100 + __clang_patchlevel__)
+#define _PSTL_CLANG_VERSION                                                    \
+  (__clang_major__ * 10000 + __clang_minor__ * 100 + __clang_patchlevel__)
 #endif
 
 // Enable SIMD for compilers that support OpenMP 4.0
-#if (defined(_OPENMP) && _OPENMP >= 201307) || \
-    (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1600) || \
-    (!defined(__INTEL_COMPILER) && _PSTL_GCC_VERSION >= 40900) || \
+#if (defined(_OPENMP) && _OPENMP >= 201307) ||                                 \
+    (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1600) ||                 \
+    (!defined(__INTEL_COMPILER) && _PSTL_GCC_VERSION >= 40900) ||              \
     defined(__clang__)
-#    define _PSTL_PRAGMA_SIMD _PSTL_PRAGMA(omp simd)
-#    define _PSTL_PRAGMA_DECLARE_SIMD _PSTL_PRAGMA(omp declare simd)
-#    define _PSTL_PRAGMA_SIMD_REDUCTION(PRM) _PSTL_PRAGMA(omp simd reduction(PRM))
+#define _PSTL_PRAGMA_SIMD _PSTL_PRAGMA(omp simd)
+#define _PSTL_PRAGMA_DECLARE_SIMD _PSTL_PRAGMA(omp declare simd)
+#define _PSTL_PRAGMA_SIMD_REDUCTION(PRM) _PSTL_PRAGMA(omp simd reduction(PRM))
 #elif !defined(_MSC_VER) //#pragma simd
-#    define _PSTL_PRAGMA_SIMD _PSTL_PRAGMA(simd)
-#    define _PSTL_PRAGMA_DECLARE_SIMD
-#    define _PSTL_PRAGMA_SIMD_REDUCTION(PRM) _PSTL_PRAGMA(simd reduction(PRM))
-#else //no simd
-#    define _PSTL_PRAGMA_SIMD
-#    define _PSTL_PRAGMA_DECLARE_SIMD
-#    define _PSTL_PRAGMA_SIMD_REDUCTION(PRM)
-#endif //Enable SIMD
+#define _PSTL_PRAGMA_SIMD _PSTL_PRAGMA(simd)
+#define _PSTL_PRAGMA_DECLARE_SIMD
+#define _PSTL_PRAGMA_SIMD_REDUCTION(PRM) _PSTL_PRAGMA(simd reduction(PRM))
+#else // no simd
+#define _PSTL_PRAGMA_SIMD
+#define _PSTL_PRAGMA_DECLARE_SIMD
+#define _PSTL_PRAGMA_SIMD_REDUCTION(PRM)
+#endif // Enable SIMD
 
 #if defined(__INTEL_COMPILER)
-#    define _PSTL_PRAGMA_FORCEINLINE _PSTL_PRAGMA(forceinline)
+#define _PSTL_PRAGMA_FORCEINLINE _PSTL_PRAGMA(forceinline)
 #else
-#    define _PSTL_PRAGMA_FORCEINLINE
+#define _PSTL_PRAGMA_FORCEINLINE
 #endif
 
 #if defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1900
-#    define _PSTL_PRAGMA_SIMD_SCAN(PRM) _PSTL_PRAGMA(omp simd reduction(inscan, PRM))
-#    define _PSTL_PRAGMA_SIMD_INCLUSIVE_SCAN(PRM) _PSTL_PRAGMA(omp scan inclusive(PRM))
-#    define _PSTL_PRAGMA_SIMD_EXCLUSIVE_SCAN(PRM) _PSTL_PRAGMA(omp scan exclusive(PRM))
+#define _PSTL_PRAGMA_SIMD_SCAN(PRM)                                            \
+  _PSTL_PRAGMA(omp simd reduction(inscan, PRM))
+#define _PSTL_PRAGMA_SIMD_INCLUSIVE_SCAN(PRM)                                  \
+  _PSTL_PRAGMA(omp scan inclusive(PRM))
+#define _PSTL_PRAGMA_SIMD_EXCLUSIVE_SCAN(PRM)                                  \
+  _PSTL_PRAGMA(omp scan exclusive(PRM))
 #else
-#    define _PSTL_PRAGMA_SIMD_SCAN(PRM)
-#    define _PSTL_PRAGMA_SIMD_INCLUSIVE_SCAN(PRM)
-#    define _PSTL_PRAGMA_SIMD_EXCLUSIVE_SCAN(PRM)
+#define _PSTL_PRAGMA_SIMD_SCAN(PRM)
+#define _PSTL_PRAGMA_SIMD_INCLUSIVE_SCAN(PRM)
+#define _PSTL_PRAGMA_SIMD_EXCLUSIVE_SCAN(PRM)
 #endif
 
-// Should be defined to 1 for environments with a vendor implementation of C++17 execution policies
-#define _PSTL_CPP17_EXECUTION_POLICIES_PRESENT (_MSC_VER >= 1912)
+// Should be defined to 1 for environments with a vendor implementation of C++17
+// execution policies
+#define _PSTL_CPP17_EXECUTION_POLICIES_PRESENT                                 \
+  (_MSC_VER >= 1912 && _MSVC_LANG >= 201703L) ||                               \
+      (_GLIBCXX_RELEASE >= 9 && __GLIBCXX__ >= 20190503 &&                     \
+       __cplusplus >= 201703L)
 
-#if (defined(_MSC_VER) && _MSC_VER >= 1900) || \
-    __cplusplus >= 201300L || \
+#if (defined(_MSC_VER) && _MSC_VER >= 1900) || __cplusplus >= 201300L ||       \
     __cpp_lib_robust_nonmodifying_seq_ops == 201304
-#   define _PSTL_CPP14_2RANGE_MISMATCH_EQUAL_PRESENT
+#define _PSTL_CPP14_2RANGE_MISMATCH_EQUAL_PRESENT
 #endif
-#if (defined(_MSC_VER) && _MSC_VER >= 1900) || \
-    __cplusplus >= 201402L || \
+#if (defined(_MSC_VER) && _MSC_VER >= 1900) || __cplusplus >= 201402L ||       \
     __cpp_lib_make_reverse_iterator == 201402
-#   define _PSTL_CPP14_MAKE_REVERSE_ITERATOR_PRESENT
+#define _PSTL_CPP14_MAKE_REVERSE_ITERATOR_PRESENT
 #endif
 #if (defined(_MSC_VER) && _MSC_VER >= 1900) || __cplusplus >= 201402L
-#   define _PSTL_CPP14_INTEGER_SEQUENCE_PRESENT
+#define _PSTL_CPP14_INTEGER_SEQUENCE_PRESENT
 #endif
-#if (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1700) || \
-    (defined(_MSC_FULL_VER) && _MSC_FULL_VER >= 190023918) || \
+#if (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1700) ||                 \
+    (defined(_MSC_FULL_VER) && _MSC_FULL_VER >= 190023918) ||                  \
     __cplusplus >= 201402L
-#   define _PSTL_CPP14_VARIABLE_TEMPLATES_PRESENT
+#define _PSTL_CPP14_VARIABLE_TEMPLATES_PRESENT
 #endif
 
 #if defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1800
-#   define _PSTL_EARLYEXIT_PRESENT
-#   define _PSTL_MONOTONIC_PRESENT
+#define _PSTL_EARLYEXIT_PRESENT
+#define _PSTL_MONOTONIC_PRESENT
 #endif
 
-#if (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1900) || \
-    (!defined(__INTEL_COMPILER) && _PSTL_GCC_VERSION >= 40900) || \
+#if (defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1900) ||                 \
+    (!defined(__INTEL_COMPILER) && _PSTL_GCC_VERSION >= 40900) ||              \
     (defined(_OPENMP) && _OPENMP >= 201307)
-#    define _PSTL_UDR_PRESENT
+#define _PSTL_UDR_PRESENT
 #endif
 
-#if defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1900 && __INTEL_COMPILER_BUILD_DATE >= 20180626
-#   define _PSTL_UDS_PRESENT
+#if defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1900 &&                   \
+    __INTEL_COMPILER_BUILD_DATE >= 20180626
+#define _PSTL_UDS_PRESENT
 #endif
 
 #if defined(_PSTL_EARLYEXIT_PRESENT)
-#    define _PSTL_PRAGMA_SIMD_EARLYEXIT _PSTL_PRAGMA(omp simd early_exit)
+#define _PSTL_PRAGMA_SIMD_EARLYEXIT _PSTL_PRAGMA(omp simd early_exit)
 #else
-#    define _PSTL_PRAGMA_SIMD_EARLYEXIT
+#define _PSTL_PRAGMA_SIMD_EARLYEXIT
 #endif
 
 #if defined(_PSTL_MONOTONIC_PRESENT)
-#    define _PSTL_PRAGMA_SIMD_ORDERED_MONOTONIC(PRM) _PSTL_PRAGMA(omp ordered simd monotonic(PRM))
-#    define _PSTL_PRAGMA_SIMD_ORDERED_MONOTONIC_2ARGS(PRM1, PRM2) _PSTL_PRAGMA(omp ordered simd monotonic(PRM1, PRM2))
+#define _PSTL_PRAGMA_SIMD_ORDERED_MONOTONIC(PRM)                               \
+  _PSTL_PRAGMA(omp ordered simd monotonic(PRM))
+#define _PSTL_PRAGMA_SIMD_ORDERED_MONOTONIC_2ARGS(PRM1, PRM2)                  \
+  _PSTL_PRAGMA(omp ordered simd monotonic(PRM1, PRM2))
 #else
-#    define _PSTL_PRAGMA_SIMD_ORDERED_MONOTONIC(PRM)
-#    define _PSTL_PRAGMA_SIMD_ORDERED_MONOTONIC_2ARGS(PRM1, PRM2)
+#define _PSTL_PRAGMA_SIMD_ORDERED_MONOTONIC(PRM)
+#define _PSTL_PRAGMA_SIMD_ORDERED_MONOTONIC_2ARGS(PRM1, PRM2)
 #endif
 
 // Declaration of reduction functor, where
 // NAME - the name of the functor
 // OP - type of the callable object with the reduction operation
 // omp_in - refers to the local partial result
 // omp_out - refers to the final value of the combiner operator
 // omp_priv - refers to the private copy of the initial value
 // omp_orig - refers to the original variable to be reduced
-#define _PSTL_PRAGMA_DECLARE_REDUCTION(NAME, OP)                                                                       \
-    _PSTL_PRAGMA(omp declare reduction(NAME:OP : omp_out(omp_in)) initializer(omp_priv = omp_orig))
+#define _PSTL_PRAGMA_DECLARE_REDUCTION(NAME, OP)                               \
+  _PSTL_PRAGMA(omp declare reduction(NAME:OP                                   \
+                                     : omp_out(omp_in))                        \
+                   initializer(omp_priv = omp_orig))
 
 #if defined(__INTEL_COMPILER) && __INTEL_COMPILER >= 1600
-#    define _PSTL_PRAGMA_VECTOR_UNALIGNED _PSTL_PRAGMA(vector unaligned)
+#define _PSTL_PRAGMA_VECTOR_UNALIGNED _PSTL_PRAGMA(vector unaligned)
 #else
-#    define _PSTL_PRAGMA_VECTOR_UNALIGNED
+#define _PSTL_PRAGMA_VECTOR_UNALIGNED
 #endif
 
 // Check the user-defined macro to use non-temporal stores
 #if defined(PSTL_USE_NONTEMPORAL_STORES) && (__INTEL_COMPILER >= 1600)
-#    define _PSTL_USE_NONTEMPORAL_STORES_IF_ALLOWED _PSTL_PRAGMA(vector nontemporal)
+#define _PSTL_USE_NONTEMPORAL_STORES_IF_ALLOWED _PSTL_PRAGMA(vector nontemporal)
 #else
-#    define _PSTL_USE_NONTEMPORAL_STORES_IF_ALLOWED
+#define _PSTL_USE_NONTEMPORAL_STORES_IF_ALLOWED
 #endif
 
-#if defined(_MSC_VER) || defined(__INTEL_COMPILER) // the preprocessors don't type a message location
-#    define _PSTL_PRAGMA_LOCATION __FILE__ ":" _PSTL_STRING(__LINE__) ": [Parallel STL message]: "
+#if defined(_MSC_VER) ||                                                       \
+    defined(                                                                   \
+        __INTEL_COMPILER) // the preprocessors don't type a message location
+#define _PSTL_PRAGMA_LOCATION                                                  \
+  __FILE__ ":" _PSTL_STRING(__LINE__) ": [Parallel STL message]: "
 #else
-#    define _PSTL_PRAGMA_LOCATION " [Parallel STL message]: "
+#define _PSTL_PRAGMA_LOCATION " [Parallel STL message]: "
 #endif
 
-#define _PSTL_PRAGMA_MESSAGE_IMPL(x) _PSTL_PRAGMA(message(_PSTL_STRING_CONCAT(_PSTL_PRAGMA_LOCATION, x)))
+#define _PSTL_PRAGMA_MESSAGE_IMPL(x)                                           \
+  _PSTL_PRAGMA(message(_PSTL_STRING_CONCAT(_PSTL_PRAGMA_LOCATION, x)))
 
 #if defined(_PSTL_USAGE_WARNINGS)
-#    define _PSTL_PRAGMA_MESSAGE(x) _PSTL_PRAGMA_MESSAGE_IMPL(x)
-#    define _PSTL_PRAGMA_MESSAGE_POLICIES(x) _PSTL_PRAGMA_MESSAGE_IMPL(x)
+#define _PSTL_PRAGMA_MESSAGE(x) _PSTL_PRAGMA_MESSAGE_IMPL(x)
+#define _PSTL_PRAGMA_MESSAGE_POLICIES(x) _PSTL_PRAGMA_MESSAGE_IMPL(x)
 #else
-#    define _PSTL_PRAGMA_MESSAGE(x)
-#    define _PSTL_PRAGMA_MESSAGE_POLICIES(x)
+#define _PSTL_PRAGMA_MESSAGE(x)
+#define _PSTL_PRAGMA_MESSAGE_POLICIES(x)
 #endif
 
 // broken macros
-#if (defined(__GLIBCXX__) && __GLIBCXX__ < 20150716) || \
+#if (defined(__GLIBCXX__) && __GLIBCXX__ < 20150716) ||                        \
     (defined(_MSC_VER) && _MSC_VER < 1800)
-#   define _PSTL_CPP11_STD_ROTATE_BROKEN
+#define _PSTL_CPP11_STD_ROTATE_BROKEN
 #endif
 
 #if defined(__INTEL_COMPILER) && __INTEL_COMPILER == 1800
-#   define _PSTL_ICC_18_OMP_SIMD_BROKEN
+#define _PSTL_ICC_18_OMP_SIMD_BROKEN
 #endif
 
 #endif /* _PSTL_CONFIG_H */