3D/4D/5D parallelization functions with 1D or no tiling

6 files changed, 3628 insertions, 0 deletions

diff --git a/include/pthreadpool.h b/include/pthreadpool.h
index de4016b..6a7d61f 100644
--- a/include/pthreadpool.h
+++ b/include/pthreadpool.h
@@ -11,8 +11,14 @@ typedef void (*pthreadpool_task_1d_tile_1d_t)(void*, size_t, size_t);
 typedef void (*pthreadpool_task_2d_t)(void*, size_t, size_t);
 typedef void (*pthreadpool_task_2d_tile_1d_t)(void*, size_t, size_t, size_t);
 typedef void (*pthreadpool_task_2d_tile_2d_t)(void*, size_t, size_t, size_t, size_t);
+typedef void (*pthreadpool_task_3d_t)(void*, size_t, size_t, size_t);
+typedef void (*pthreadpool_task_3d_tile_1d_t)(void*, size_t, size_t, size_t, size_t);
 typedef void (*pthreadpool_task_3d_tile_2d_t)(void*, size_t, size_t, size_t, size_t, size_t);
+typedef void (*pthreadpool_task_4d_t)(void*, size_t, size_t, size_t, size_t);
+typedef void (*pthreadpool_task_4d_tile_1d_t)(void*, size_t, size_t, size_t, size_t, size_t);
 typedef void (*pthreadpool_task_4d_tile_2d_t)(void*, size_t, size_t, size_t, size_t, size_t, size_t);
+typedef void (*pthreadpool_task_5d_t)(void*, size_t, size_t, size_t, size_t, size_t);
+typedef void (*pthreadpool_task_5d_tile_1d_t)(void*, size_t, size_t, size_t, size_t, size_t, size_t);
 typedef void (*pthreadpool_task_5d_tile_2d_t)(void*, size_t, size_t, size_t, size_t, size_t, size_t, size_t);
 typedef void (*pthreadpool_task_6d_tile_2d_t)(void*, size_t, size_t, size_t, size_t, size_t, size_t, size_t, size_t);
 
@@ -360,6 +366,86 @@ void pthreadpool_parallelize_2d_tile_2d_with_uarch(
 	uint32_t flags);
 
 /**
+ * Process items on a 3D grid.
+ *
+ * The function implements a parallel version of the following snippet:
+ *
+ *   for (size_t i = 0; i < range_i; i++)
+ *     for (size_t j = 0; j < range_j; j++)
+ *       for (size_t k = 0; k < range_k; k++)
+ *         function(context, i, j, k);
+ *
+ * When the function returns, all items have been processed and the thread pool
+ * is ready for a new task.
+ *
+ * @note If multiple threads call this function with the same thread pool, the
+ *    calls are serialized.
+ *
+ * @param threadpool  the thread pool to use for parallelisation. If threadpool
+ *    is NULL, all items are processed serially on the calling thread.
+ * @param function    the function to call for each tile.
+ * @param context     the first argument passed to the specified function.
+ * @param range_i     the number of items to process along the first dimension
+ *    of the 3D grid.
+ * @param range_j     the number of items to process along the second dimension
+ *    of the 3D grid.
+ * @param range_k     the number of items to process along the third dimension
+ *    of the 3D grid.
+ * @param flags       a bitwise combination of zero or more optional flags
+ *    (PTHREADPOOL_FLAG_DISABLE_DENORMALS or PTHREADPOOL_FLAG_YIELD_WORKERS)
+ */
+void pthreadpool_parallelize_3d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_3d_t function,
+	void* context,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	uint32_t flags);
+
+/**
+ * Process items on a 3D grid with the specified maximum tile size along the
+ * last grid dimension.
+ *
+ * The function implements a parallel version of the following snippet:
+ *
+ *   for (size_t i = 0; i < range_i; i++)
+ *     for (size_t j = 0; j < range_j; j++)
+ *       for (size_t k = 0; k < range_k; k += tile_k)
+ *         function(context, i, j, k, min(range_k - k, tile_k));
+ *
+ * When the function returns, all items have been processed and the thread pool
+ * is ready for a new task.
+ *
+ * @note If multiple threads call this function with the same thread pool, the
+ *    calls are serialized.
+ *
+ * @param threadpool  the thread pool to use for parallelisation. If threadpool
+ *    is NULL, all items are processed serially on the calling thread.
+ * @param function    the function to call for each tile.
+ * @param context     the first argument passed to the specified function.
+ * @param range_i     the number of items to process along the first dimension
+ *    of the 3D grid.
+ * @param range_j     the number of items to process along the second dimension
+ *    of the 3D grid.
+ * @param range_k     the number of items to process along the third dimension
+ *    of the 3D grid.
+ * @param tile_k      the maximum number of items along the third dimension of
+ *    the 3D grid to process in one function call.
+ * @param flags       a bitwise combination of zero or more optional flags
+ *    (PTHREADPOOL_FLAG_DISABLE_DENORMALS or PTHREADPOOL_FLAG_YIELD_WORKERS)
+ */
+void pthreadpool_parallelize_3d_tile_1d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_3d_tile_1d_t function,
+	void* context,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t tile_k,
+	uint32_t flags);
+
+/**
  * Process items on a 3D grid with the specified maximum tile size along the
  * last two grid dimensions.
  *
@@ -468,6 +554,94 @@ void pthreadpool_parallelize_3d_tile_2d_with_uarch(
 	uint32_t flags);
 
 /**
+ * Process items on a 4D grid.
+ *
+ * The function implements a parallel version of the following snippet:
+ *
+ *   for (size_t i = 0; i < range_i; i++)
+ *     for (size_t j = 0; j < range_j; j++)
+ *       for (size_t k = 0; k < range_k; k++)
+ *         for (size_t l = 0; l < range_l; l++)
+ *           function(context, i, j, k, l);
+ *
+ * When the function returns, all items have been processed and the thread pool
+ * is ready for a new task.
+ *
+ * @note If multiple threads call this function with the same thread pool, the
+ *    calls are serialized.
+ *
+ * @param threadpool  the thread pool to use for parallelisation. If threadpool
+ *    is NULL, all items are processed serially on the calling thread.
+ * @param function    the function to call for each tile.
+ * @param context     the first argument passed to the specified function.
+ * @param range_i     the number of items to process along the first dimension
+ *    of the 4D grid.
+ * @param range_j     the number of items to process along the second dimension
+ *    of the 4D grid.
+ * @param range_k     the number of items to process along the third dimension
+ *    of the 4D grid.
+ * @param range_l     the number of items to process along the fourth dimension
+ *    of the 4D grid.
+ * @param flags       a bitwise combination of zero or more optional flags
+ *    (PTHREADPOOL_FLAG_DISABLE_DENORMALS or PTHREADPOOL_FLAG_YIELD_WORKERS)
+ */
+void pthreadpool_parallelize_4d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_4d_t function,
+	void* context,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	uint32_t flags);
+
+/**
+ * Process items on a 4D grid with the specified maximum tile size along the
+ * last grid dimension.
+ *
+ * The function implements a parallel version of the following snippet:
+ *
+ *   for (size_t i = 0; i < range_i; i++)
+ *     for (size_t j = 0; j < range_j; j++)
+ *       for (size_t k = 0; k < range_k; k++)
+ *         for (size_t l = 0; l < range_l; l += tile_l)
+ *           function(context, i, j, k, l, min(range_l - l, tile_l));
+ *
+ * When the function returns, all items have been processed and the thread pool
+ * is ready for a new task.
+ *
+ * @note If multiple threads call this function with the same thread pool, the
+ *    calls are serialized.
+ *
+ * @param threadpool  the thread pool to use for parallelisation. If threadpool
+ *    is NULL, all items are processed serially on the calling thread.
+ * @param function    the function to call for each tile.
+ * @param context     the first argument passed to the specified function.
+ * @param range_i     the number of items to process along the first dimension
+ *    of the 4D grid.
+ * @param range_j     the number of items to process along the second dimension
+ *    of the 4D grid.
+ * @param range_k     the number of items to process along the third dimension
+ *    of the 4D grid.
+ * @param range_l     the number of items to process along the fourth dimension
+ *    of the 4D grid.
+ * @param tile_l      the maximum number of items along the fourth dimension of
+ *    the 4D grid to process in one function call.
+ * @param flags       a bitwise combination of zero or more optional flags
+ *    (PTHREADPOOL_FLAG_DISABLE_DENORMALS or PTHREADPOOL_FLAG_YIELD_WORKERS)
+ */
+void pthreadpool_parallelize_4d_tile_1d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_4d_tile_1d_t function,
+	void* context,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	size_t tile_l,
+	uint32_t flags);
+
+/**
  * Process items on a 4D grid with the specified maximum tile size along the
  * last two grid dimensions.
  *
@@ -584,6 +758,102 @@ void pthreadpool_parallelize_4d_tile_2d_with_uarch(
 	uint32_t flags);
 
 /**
+ * Process items on a 5D grid.
+ *
+ * The function implements a parallel version of the following snippet:
+ *
+ *   for (size_t i = 0; i < range_i; i++)
+ *     for (size_t j = 0; j < range_j; j++)
+ *       for (size_t k = 0; k < range_k; k++)
+ *         for (size_t l = 0; l < range_l; l++)
+ *           for (size_t m = 0; m < range_m; m++)
+ *             function(context, i, j, k, l, m);
+ *
+ * When the function returns, all items have been processed and the thread pool
+ * is ready for a new task.
+ *
+ * @note If multiple threads call this function with the same thread pool, the
+ *    calls are serialized.
+ *
+ * @param threadpool  the thread pool to use for parallelisation. If threadpool
+ *    is NULL, all items are processed serially on the calling thread.
+ * @param function    the function to call for each tile.
+ * @param context     the first argument passed to the specified function.
+ * @param range_i     the number of items to process along the first dimension
+ *    of the 5D grid.
+ * @param range_j     the number of items to process along the second dimension
+ *    of the 5D grid.
+ * @param range_k     the number of items to process along the third dimension
+ *    of the 5D grid.
+ * @param range_l     the number of items to process along the fourth dimension
+ *    of the 5D grid.
+ * @param range_m     the number of items to process along the fifth dimension
+ *    of the 5D grid.
+ * @param flags       a bitwise combination of zero or more optional flags
+ *    (PTHREADPOOL_FLAG_DISABLE_DENORMALS or PTHREADPOOL_FLAG_YIELD_WORKERS)
+ */
+void pthreadpool_parallelize_5d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_5d_t function,
+	void* context,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	size_t range_m,
+	uint32_t flags);
+
+/**
+ * Process items on a 5D grid with the specified maximum tile size along the
+ * last grid dimension.
+ *
+ * The function implements a parallel version of the following snippet:
+ *
+ *   for (size_t i = 0; i < range_i; i++)
+ *     for (size_t j = 0; j < range_j; j++)
+ *       for (size_t k = 0; k < range_k; k++)
+ *         for (size_t l = 0; l < range_l; l++)
+ *           for (size_t m = 0; m < range_m; m += tile_m)
+ *             function(context, i, j, k, l, m, min(range_m - m, tile_m));
+ *
+ * When the function returns, all items have been processed and the thread pool
+ * is ready for a new task.
+ *
+ * @note If multiple threads call this function with the same thread pool, the
+ *    calls are serialized.
+ *
+ * @param threadpool  the thread pool to use for parallelisation. If threadpool
+ *    is NULL, all items are processed serially on the calling thread.
+ * @param function    the function to call for each tile.
+ * @param context     the first argument passed to the specified function.
+ * @param range_i     the number of items to process along the first dimension
+ *    of the 5D grid.
+ * @param range_j     the number of items to process along the second dimension
+ *    of the 5D grid.
+ * @param range_k     the number of items to process along the third dimension
+ *    of the 5D grid.
+ * @param range_l     the number of items to process along the fourth dimension
+ *    of the 5D grid.
+ * @param range_m     the number of items to process along the fifth dimension
+ *    of the 5D grid.
+ * @param tile_m      the maximum number of items along the fifth dimension of
+ *    the 5D grid to process in one function call.
+ * @param flags       a bitwise combination of zero or more optional flags
+ *    (PTHREADPOOL_FLAG_DISABLE_DENORMALS or PTHREADPOOL_FLAG_YIELD_WORKERS)
+ */
+void pthreadpool_parallelize_5d_tile_1d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_5d_tile_1d_t function,
+	void* context,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	size_t range_m,
+	size_t tile_m,
+	uint32_t flags);
+
+/**
  * Process items on a 5D grid with the specified maximum tile size along the
  * last two grid dimensions.
  *
diff --git a/src/fastpath.c b/src/fastpath.c
index 1a5066a..b4e40c5 100644
--- a/src/fastpath.c
+++ b/src/fastpath.c
@@ -356,6 +356,116 @@ PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_2d_tile_2d_with_uarch_f
 	pthreadpool_fence_release();
 }
 
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_3d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread)
+{
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_3d_t task = (pthreadpool_task_3d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	const size_t threads_count = threadpool->threads_count.value;
+	const size_t range_threshold = -threads_count;
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t range_k = threadpool->params.parallelize_3d.range_k;
+	const struct fxdiv_result_size_t index_ij_k = fxdiv_divide_size_t(range_start, range_k);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_3d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_k.quotient, range_j);
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t k = index_ij_k.remainder;
+
+	while (pthreadpool_decrement_fetch_relaxed_size_t(&thread->range_length) < range_threshold) {
+		task(argument, i, j, k);
+		if (++k == range_k.value) {
+			k = 0;
+			if (++j == range_j.value) {
+				j = 0;
+				i += 1;
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_length) < range_threshold) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t index_ij_k = fxdiv_divide_size_t(linear_index, range_k);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_k.quotient, range_j);
+			task(argument, index_i_j.quotient, index_i_j.remainder, index_ij_k.remainder);
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_3d_tile_1d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread)
+{
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_3d_tile_1d_t task = (pthreadpool_task_3d_tile_1d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	const size_t threads_count = threadpool->threads_count.value;
+	const size_t range_threshold = -threads_count;
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t tile_range_k = threadpool->params.parallelize_3d_tile_1d.tile_range_k;
+	const struct fxdiv_result_size_t tile_index_ij_k = fxdiv_divide_size_t(range_start, tile_range_k);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_3d_tile_1d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(tile_index_ij_k.quotient, range_j);
+	const size_t tile_k = threadpool->params.parallelize_3d_tile_1d.tile_k;
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t start_k = tile_index_ij_k.remainder * tile_k;
+
+	const size_t range_k = threadpool->params.parallelize_3d_tile_1d.range_k;
+	while (pthreadpool_decrement_fetch_relaxed_size_t(&thread->range_length) < range_threshold) {
+		task(argument, i, j, start_k, min(range_k - start_k, tile_k));
+		start_k += tile_k;
+		if (start_k >= range_k) {
+			start_k = 0;
+			if (++j == range_j.value) {
+				j = 0;
+				i += 1;
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_length) < range_threshold) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t tile_index_ij_k = fxdiv_divide_size_t(linear_index, tile_range_k);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(tile_index_ij_k.quotient, range_j);
+			const size_t start_k = tile_index_ij_k.remainder * tile_k;
+			task(argument, index_i_j.quotient, index_i_j.remainder, start_k, min(range_k - start_k, tile_k));
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
 PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_3d_tile_2d_fastpath(
 	struct pthreadpool* threadpool,
 	struct thread_info* thread)
@@ -487,6 +597,132 @@ PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_3d_tile_2d_with_uarch_f
 	pthreadpool_fence_release();
 }
 
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_4d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread)
+{
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_4d_t task = (pthreadpool_task_4d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	const size_t threads_count = threadpool->threads_count.value;
+	const size_t range_threshold = -threads_count;
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t range_kl = threadpool->params.parallelize_4d.range_kl;
+	const struct fxdiv_result_size_t index_ij_kl = fxdiv_divide_size_t(range_start, range_kl);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_4d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_kl.quotient, range_j);
+	const struct fxdiv_divisor_size_t range_l = threadpool->params.parallelize_4d.range_l;
+	const struct fxdiv_result_size_t index_k_l = fxdiv_divide_size_t(index_ij_kl.remainder, range_l);
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t k = index_k_l.quotient;
+	size_t l = index_k_l.remainder;
+
+	const size_t range_k = threadpool->params.parallelize_4d.range_k;
+	while (pthreadpool_decrement_fetch_relaxed_size_t(&thread->range_length) < range_threshold) {
+		task(argument, i, j, k, l);
+		if (++l == range_l.value) {
+			l = 0;
+			if (++k == range_k) {
+				k = 0;
+				if (++j == range_j.value) {
+					j = 0;
+					i += 1;
+				}
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_length) < range_threshold) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t index_ij_kl = fxdiv_divide_size_t(linear_index, range_kl);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_kl.quotient, range_j);
+			const struct fxdiv_result_size_t index_k_l = fxdiv_divide_size_t(index_ij_kl.remainder, range_l);
+			task(argument, index_i_j.quotient, index_i_j.remainder, index_k_l.quotient, index_k_l.remainder);
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_4d_tile_1d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread)
+{
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_4d_tile_1d_t task = (pthreadpool_task_4d_tile_1d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	const size_t threads_count = threadpool->threads_count.value;
+	const size_t range_threshold = -threads_count;
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t tile_range_kl = threadpool->params.parallelize_4d_tile_1d.tile_range_kl;
+	const struct fxdiv_result_size_t tile_index_ij_kl = fxdiv_divide_size_t(range_start, tile_range_kl);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_4d_tile_1d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(tile_index_ij_kl.quotient, range_j);
+	const struct fxdiv_divisor_size_t tile_range_l = threadpool->params.parallelize_4d_tile_1d.tile_range_l;
+	const struct fxdiv_result_size_t tile_index_k_l = fxdiv_divide_size_t(tile_index_ij_kl.remainder, tile_range_l);
+	const size_t tile_l = threadpool->params.parallelize_4d_tile_1d.tile_l;
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t k = tile_index_k_l.quotient;
+	size_t start_l = tile_index_k_l.remainder * tile_l;
+
+	const size_t range_l = threadpool->params.parallelize_4d_tile_1d.range_l;
+	const size_t range_k = threadpool->params.parallelize_4d_tile_1d.range_k;
+	while (pthreadpool_decrement_fetch_relaxed_size_t(&thread->range_length) < range_threshold) {
+		task(argument, i, j, k, start_l, min(range_l - start_l, tile_l));
+		start_l += tile_l;
+		if (start_l >= range_l) {
+			start_l = 0;
+			if (++k == range_k) {
+				k = 0;
+				if (++j == range_j.value) {
+					j = 0;
+					i += 1;
+				}
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_length) < range_threshold) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t tile_index_ij_kl = fxdiv_divide_size_t(linear_index, tile_range_kl);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(tile_index_ij_kl.quotient, range_j);
+			const struct fxdiv_result_size_t tile_index_k_l = fxdiv_divide_size_t(tile_index_ij_kl.remainder, tile_range_l);
+			const size_t start_l = tile_index_k_l.remainder * tile_l;
+			task(argument, index_i_j.quotient, index_i_j.remainder, tile_index_k_l.quotient, start_l, min(range_l - start_l, tile_l));
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
 PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_4d_tile_2d_fastpath(
 	struct pthreadpool* threadpool,
 	struct thread_info* thread)
@@ -632,6 +868,147 @@ PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_4d_tile_2d_with_uarch_f
 	pthreadpool_fence_release();
 }
 
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_5d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread)
+{
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_5d_t task = (pthreadpool_task_5d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	const size_t threads_count = threadpool->threads_count.value;
+	const size_t range_threshold = -threads_count;
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t range_lm = threadpool->params.parallelize_5d.range_lm;
+	const struct fxdiv_result_size_t index_ijk_lm = fxdiv_divide_size_t(range_start, range_lm);
+	const struct fxdiv_divisor_size_t range_k = threadpool->params.parallelize_5d.range_k;
+	const struct fxdiv_result_size_t index_ij_k = fxdiv_divide_size_t(index_ijk_lm.quotient, range_k);
+	const struct fxdiv_divisor_size_t range_m = threadpool->params.parallelize_5d.range_m;
+	const struct fxdiv_result_size_t index_l_m = fxdiv_divide_size_t(index_ijk_lm.remainder, range_m);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_5d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_k.quotient, range_j);
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t k = index_ij_k.remainder;
+	size_t l = index_l_m.quotient;
+	size_t m = index_l_m.remainder;
+
+	const size_t range_l = threadpool->params.parallelize_5d.range_l;
+	while (pthreadpool_decrement_fetch_relaxed_size_t(&thread->range_length) < range_threshold) {
+		task(argument, i, j, k, l, m);
+		if (++m == range_m.value) {
+			m = 0;
+			if (++l == range_l) {
+				l = 0;
+				if (++k == range_k.value) {
+					k = 0;
+					if (++j == range_j.value) {
+						j = 0;
+						i += 1;
+					}
+				}
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_length) < range_threshold) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t index_ijk_lm = fxdiv_divide_size_t(linear_index, range_lm);
+			const struct fxdiv_result_size_t index_ij_k = fxdiv_divide_size_t(index_ijk_lm.quotient, range_k);
+			const struct fxdiv_result_size_t index_l_m = fxdiv_divide_size_t(index_ijk_lm.remainder, range_m);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_k.quotient, range_j);
+			task(argument, index_i_j.quotient, index_i_j.remainder, index_ij_k.remainder, index_l_m.quotient, index_l_m.remainder);
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_5d_tile_1d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread)
+{
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_5d_tile_1d_t task = (pthreadpool_task_5d_tile_1d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	const size_t threads_count = threadpool->threads_count.value;
+	const size_t range_threshold = -threads_count;
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t tile_range_m = threadpool->params.parallelize_5d_tile_1d.tile_range_m;
+	const struct fxdiv_result_size_t tile_index_ijkl_m = fxdiv_divide_size_t(range_start, tile_range_m);
+	const struct fxdiv_divisor_size_t range_kl = threadpool->params.parallelize_5d_tile_1d.range_kl;
+	const struct fxdiv_result_size_t index_ij_kl = fxdiv_divide_size_t(tile_index_ijkl_m.quotient, range_kl);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_5d_tile_1d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_kl.quotient, range_j);
+	const struct fxdiv_divisor_size_t range_l = threadpool->params.parallelize_5d_tile_1d.range_l;
+	const struct fxdiv_result_size_t index_k_l = fxdiv_divide_size_t(index_ij_kl.remainder, range_l);
+	const size_t tile_m = threadpool->params.parallelize_5d_tile_1d.tile_m;
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t k = index_k_l.quotient;
+	size_t l = index_k_l.remainder;
+	size_t start_m = tile_index_ijkl_m.remainder * tile_m;
+
+	const size_t range_m = threadpool->params.parallelize_5d_tile_1d.range_m;
+	const size_t range_k = threadpool->params.parallelize_5d_tile_1d.range_k;
+	while (pthreadpool_decrement_fetch_relaxed_size_t(&thread->range_length) < range_threshold) {
+		task(argument, i, j, k, l, start_m, min(range_m - start_m, tile_m));
+		start_m += tile_m;
+		if (start_m >= range_m) {
+			start_m = 0;
+			if (++l == range_l.value) {
+				l = 0;
+				if (++k == range_k) {
+					k = 0;
+					if (++j == range_j.value) {
+						j = 0;
+						i += 1;
+					}
+				}
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_length) < range_threshold) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t tile_index_ijkl_m = fxdiv_divide_size_t(linear_index, tile_range_m);
+			const struct fxdiv_result_size_t index_ij_kl = fxdiv_divide_size_t(tile_index_ijkl_m.quotient, range_kl);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_kl.quotient, range_j);
+			const struct fxdiv_result_size_t index_k_l = fxdiv_divide_size_t(index_ij_kl.remainder, range_l);
+			size_t start_m = tile_index_ijkl_m.remainder * tile_m;
+			task(argument, index_i_j.quotient, index_i_j.remainder, index_k_l.quotient, index_k_l.remainder, start_m,
+				min(range_m - start_m, tile_m));
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
 PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_5d_tile_2d_fastpath(
 	struct pthreadpool* threadpool,
 	struct thread_info* thread)
diff --git a/src/portable-api.c b/src/portable-api.c
index 6b16674..6a8ccf2 100644
--- a/src/portable-api.c
+++ b/src/portable-api.c
@@ -328,6 +328,106 @@ static void thread_parallelize_2d_tile_2d_with_uarch(struct pthreadpool* threadp
 	pthreadpool_fence_release();
 }
 
+static void thread_parallelize_3d(struct pthreadpool* threadpool, struct thread_info* thread) {
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_3d_t task = (pthreadpool_task_3d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t range_k = threadpool->params.parallelize_3d.range_k;
+	const struct fxdiv_result_size_t index_ij_k = fxdiv_divide_size_t(range_start, range_k);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_3d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_k.quotient, range_j);
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t k = index_ij_k.remainder;
+
+	while (pthreadpool_try_decrement_relaxed_size_t(&thread->range_length)) {
+		task(argument, i, j, k);
+		if (++k == range_k.value) {
+			k = 0;
+			if (++j == range_j.value) {
+				j = 0;
+				i += 1;
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	const size_t threads_count = threadpool->threads_count.value;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_try_decrement_relaxed_size_t(&other_thread->range_length)) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t index_ij_k = fxdiv_divide_size_t(linear_index, range_k);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_k.quotient, range_j);
+			task(argument, index_i_j.quotient, index_i_j.remainder, index_ij_k.remainder);
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
+static void thread_parallelize_3d_tile_1d(struct pthreadpool* threadpool, struct thread_info* thread) {
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_3d_tile_1d_t task = (pthreadpool_task_3d_tile_1d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t tile_range_k = threadpool->params.parallelize_3d_tile_1d.tile_range_k;
+	const struct fxdiv_result_size_t tile_index_ij_k = fxdiv_divide_size_t(range_start, tile_range_k);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_3d_tile_1d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(tile_index_ij_k.quotient, range_j);
+	const size_t tile_k = threadpool->params.parallelize_3d_tile_1d.tile_k;
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t start_k = tile_index_ij_k.remainder * tile_k;
+
+	const size_t range_k = threadpool->params.parallelize_3d_tile_1d.range_k;
+	while (pthreadpool_try_decrement_relaxed_size_t(&thread->range_length)) {
+		task(argument, i, j, start_k, min(range_k - start_k, tile_k));
+		start_k += tile_k;
+		if (start_k >= range_k) {
+			start_k = 0;
+			if (++j == range_j.value) {
+				j = 0;
+				i += 1;
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	const size_t threads_count = threadpool->threads_count.value;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_try_decrement_relaxed_size_t(&other_thread->range_length)) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t tile_index_ij_k = fxdiv_divide_size_t(linear_index, tile_range_k);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(tile_index_ij_k.quotient, range_j);
+			const size_t start_k = tile_index_ij_k.remainder * tile_k;
+			task(argument, index_i_j.quotient, index_i_j.remainder, start_k, min(range_k - start_k, tile_k));
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
 static void thread_parallelize_3d_tile_2d(struct pthreadpool* threadpool, struct thread_info* thread) {
 	assert(threadpool != NULL);
 	assert(thread != NULL);
@@ -449,6 +549,122 @@ static void thread_parallelize_3d_tile_2d_with_uarch(struct pthreadpool* threadp
 	pthreadpool_fence_release();
 }
 
+static void thread_parallelize_4d(struct pthreadpool* threadpool, struct thread_info* thread) {
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_4d_t task = (pthreadpool_task_4d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t range_kl = threadpool->params.parallelize_4d.range_kl;
+	const struct fxdiv_result_size_t index_ij_kl = fxdiv_divide_size_t(range_start, range_kl);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_4d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_kl.quotient, range_j);
+	const struct fxdiv_divisor_size_t range_l = threadpool->params.parallelize_4d.range_l;
+	const struct fxdiv_result_size_t index_k_l = fxdiv_divide_size_t(index_ij_kl.remainder, range_l);
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t k = index_k_l.quotient;
+	size_t l = index_k_l.remainder;
+
+	const size_t range_k = threadpool->params.parallelize_4d.range_k;
+	while (pthreadpool_try_decrement_relaxed_size_t(&thread->range_length)) {
+		task(argument, i, j, k, l);
+		if (++l == range_l.value) {
+			l = 0;
+			if (++k == range_k) {
+				k = 0;
+				if (++j == range_j.value) {
+					j = 0;
+					i += 1;
+				}
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	const size_t threads_count = threadpool->threads_count.value;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_try_decrement_relaxed_size_t(&other_thread->range_length)) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t index_ij_kl = fxdiv_divide_size_t(linear_index, range_kl);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_kl.quotient, range_j);
+			const struct fxdiv_result_size_t index_k_l = fxdiv_divide_size_t(index_ij_kl.remainder, range_l);
+			task(argument, index_i_j.quotient, index_i_j.remainder, index_k_l.quotient, index_k_l.remainder);
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
+static void thread_parallelize_4d_tile_1d(struct pthreadpool* threadpool, struct thread_info* thread) {
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_4d_tile_1d_t task = (pthreadpool_task_4d_tile_1d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t tile_range_kl = threadpool->params.parallelize_4d_tile_1d.tile_range_kl;
+	const struct fxdiv_result_size_t tile_index_ij_kl = fxdiv_divide_size_t(range_start, tile_range_kl);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_4d_tile_1d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(tile_index_ij_kl.quotient, range_j);
+	const struct fxdiv_divisor_size_t tile_range_l = threadpool->params.parallelize_4d_tile_1d.tile_range_l;
+	const struct fxdiv_result_size_t tile_index_k_l = fxdiv_divide_size_t(tile_index_ij_kl.remainder, tile_range_l);
+	const size_t tile_l = threadpool->params.parallelize_4d_tile_1d.tile_l;
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t k = tile_index_k_l.quotient;
+	size_t start_l = tile_index_k_l.remainder * tile_l;
+
+	const size_t range_k = threadpool->params.parallelize_4d_tile_1d.range_k;
+	const size_t range_l = threadpool->params.parallelize_4d_tile_1d.range_l;
+	while (pthreadpool_try_decrement_relaxed_size_t(&thread->range_length)) {
+		task(argument, i, j, k, start_l, min(range_l - start_l, tile_l));
+		start_l += tile_l;
+		if (start_l >= range_l) {
+			start_l = 0;
+			if (++k == range_k) {
+				k = 0;
+				if (++j == range_j.value) {
+					j = 0;
+					i += 1;
+				}
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	const size_t threads_count = threadpool->threads_count.value;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_try_decrement_relaxed_size_t(&other_thread->range_length)) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t tile_index_ij_kl = fxdiv_divide_size_t(linear_index, tile_range_kl);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(tile_index_ij_kl.quotient, range_j);
+			const struct fxdiv_result_size_t tile_index_k_l = fxdiv_divide_size_t(tile_index_ij_kl.remainder, tile_range_l);
+			const size_t start_l = tile_index_k_l.remainder * tile_l;
+			task(argument, index_i_j.quotient, index_i_j.remainder, tile_index_k_l.quotient, start_l, min(range_l - start_l, tile_l));
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
 static void thread_parallelize_4d_tile_2d(struct pthreadpool* threadpool, struct thread_info* thread) {
 	assert(threadpool != NULL);
 	assert(thread != NULL);
@@ -584,6 +800,137 @@ static void thread_parallelize_4d_tile_2d_with_uarch(struct pthreadpool* threadp
 	pthreadpool_fence_release();
 }
 
+static void thread_parallelize_5d(struct pthreadpool* threadpool, struct thread_info* thread) {
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_5d_t task = (pthreadpool_task_5d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t range_lm = threadpool->params.parallelize_5d.range_lm;
+	const struct fxdiv_result_size_t index_ijk_lm = fxdiv_divide_size_t(range_start, range_lm);
+	const struct fxdiv_divisor_size_t range_k = threadpool->params.parallelize_5d.range_k;
+	const struct fxdiv_result_size_t index_ij_k = fxdiv_divide_size_t(index_ijk_lm.quotient, range_k);
+	const struct fxdiv_divisor_size_t range_m = threadpool->params.parallelize_5d.range_m;
+	const struct fxdiv_result_size_t index_l_m = fxdiv_divide_size_t(index_ijk_lm.remainder, range_m);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_5d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_k.quotient, range_j);
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t k = index_ij_k.remainder;
+	size_t l = index_l_m.quotient;
+	size_t m = index_l_m.remainder;
+
+	const size_t range_l = threadpool->params.parallelize_5d.range_l;
+	while (pthreadpool_try_decrement_relaxed_size_t(&thread->range_length)) {
+		task(argument, i, j, k, l, m);
+		if (++m == range_m.value) {
+			m = 0;
+			if (++l == range_l) {
+				l = 0;
+				if (++k == range_k.value) {
+					k = 0;
+					if (++j == range_j.value) {
+						j = 0;
+						i += 1;
+					}
+				}
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	const size_t threads_count = threadpool->threads_count.value;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_try_decrement_relaxed_size_t(&other_thread->range_length)) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t index_ijk_lm = fxdiv_divide_size_t(linear_index, range_lm);
+			const struct fxdiv_result_size_t index_ij_k = fxdiv_divide_size_t(index_ijk_lm.quotient, range_k);
+			const struct fxdiv_result_size_t index_l_m = fxdiv_divide_size_t(index_ijk_lm.remainder, range_m);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_k.quotient, range_j);
+			task(argument, index_i_j.quotient, index_i_j.remainder, index_ij_k.remainder, index_l_m.quotient, index_l_m.remainder);
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
+static void thread_parallelize_5d_tile_1d(struct pthreadpool* threadpool, struct thread_info* thread) {
+	assert(threadpool != NULL);
+	assert(thread != NULL);
+
+	const pthreadpool_task_5d_tile_1d_t task = (pthreadpool_task_5d_tile_1d_t) pthreadpool_load_relaxed_void_p(&threadpool->task);
+	void *const argument = pthreadpool_load_relaxed_void_p(&threadpool->argument);
+
+	/* Process thread's own range of items */
+	const size_t range_start = pthreadpool_load_relaxed_size_t(&thread->range_start);
+	const struct fxdiv_divisor_size_t tile_range_m = threadpool->params.parallelize_5d_tile_1d.tile_range_m;
+	const struct fxdiv_result_size_t tile_index_ijkl_m = fxdiv_divide_size_t(range_start, tile_range_m);
+	const struct fxdiv_divisor_size_t range_kl = threadpool->params.parallelize_5d_tile_1d.range_kl;
+	const struct fxdiv_result_size_t index_ij_kl = fxdiv_divide_size_t(tile_index_ijkl_m.quotient, range_kl);
+	const struct fxdiv_divisor_size_t range_j = threadpool->params.parallelize_5d_tile_1d.range_j;
+	const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_kl.quotient, range_j);
+	const struct fxdiv_divisor_size_t range_l = threadpool->params.parallelize_5d_tile_1d.range_l;
+	const struct fxdiv_result_size_t index_k_l = fxdiv_divide_size_t(index_ij_kl.remainder, range_l);
+	const size_t tile_m = threadpool->params.parallelize_5d_tile_1d.tile_m;
+	size_t i = index_i_j.quotient;
+	size_t j = index_i_j.remainder;
+	size_t k = index_k_l.quotient;
+	size_t l = index_k_l.remainder;
+	size_t start_m = tile_index_ijkl_m.remainder * tile_m;
+
+	const size_t range_m = threadpool->params.parallelize_5d_tile_1d.range_m;
+	const size_t range_k = threadpool->params.parallelize_5d_tile_1d.range_k;
+	while (pthreadpool_try_decrement_relaxed_size_t(&thread->range_length)) {
+		task(argument, i, j, k, l, start_m, min(range_m - start_m, tile_m));
+		start_m += tile_m;
+		if (start_m >= range_m) {
+			start_m = 0;
+			if (++l == range_l.value) {
+				l = 0;
+				if (++k == range_k) {
+					k = 0;
+					if (++j == range_j.value) {
+						j = 0;
+						i += 1;
+					}
+				}
+			}
+		}
+	}
+
+	/* There still may be other threads with work */
+	const size_t thread_number = thread->thread_number;
+	const size_t threads_count = threadpool->threads_count.value;
+	for (size_t tid = modulo_decrement(thread_number, threads_count);
+		tid != thread_number;
+		tid = modulo_decrement(tid, threads_count))
+	{
+		struct thread_info* other_thread = &threadpool->threads[tid];
+		while (pthreadpool_try_decrement_relaxed_size_t(&other_thread->range_length)) {
+			const size_t linear_index = pthreadpool_decrement_fetch_relaxed_size_t(&other_thread->range_end);
+			const struct fxdiv_result_size_t tile_index_ijkl_m = fxdiv_divide_size_t(linear_index, tile_range_m);
+			const struct fxdiv_result_size_t index_ij_kl = fxdiv_divide_size_t(tile_index_ijkl_m.quotient, range_kl);
+			const struct fxdiv_result_size_t index_i_j = fxdiv_divide_size_t(index_ij_kl.quotient, range_j);
+			const struct fxdiv_result_size_t index_k_l = fxdiv_divide_size_t(index_ij_kl.remainder, range_l);
+			size_t start_m = tile_index_ijkl_m.remainder * tile_m;
+			task(argument, index_i_j.quotient, index_i_j.remainder, index_k_l.quotient, index_k_l.remainder, start_m,
+				min(range_m - start_m, tile_m));
+		}
+	}
+
+	/* Make changes by this thread visible to other threads */
+	pthreadpool_fence_release();
+}
+
 static void thread_parallelize_5d_tile_2d(struct pthreadpool* threadpool, struct thread_info* thread) {
 	assert(threadpool != NULL);
 	assert(thread != NULL);
@@ -1061,6 +1408,102 @@ void pthreadpool_parallelize_2d_tile_2d_with_uarch(
 	}
 }
 
+void pthreadpool_parallelize_3d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_3d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	uint32_t flags)
+{
+	size_t threads_count;
+	if (threadpool == NULL || (threads_count = threadpool->threads_count.value) <= 1 || (range_i | range_j | range_k) <= 1) {
+		/* No thread pool used: execute task sequentially on the calling thread */
+		struct fpu_state saved_fpu_state = { 0 };
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			saved_fpu_state = get_fpu_state();
+			disable_fpu_denormals();
+		}
+		for (size_t i = 0; i < range_i; i++) {
+			for (size_t j = 0; j < range_j; j++) {
+				for (size_t k = 0; k < range_k; k++) {
+					task(argument, i, j, k);
+				}
+			}
+		}
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			set_fpu_state(saved_fpu_state);
+		}
+	} else {
+		const size_t range = range_i * range_j * range_k;
+		const struct pthreadpool_3d_params params = {
+			.range_j = fxdiv_init_size_t(range_j),
+			.range_k = fxdiv_init_size_t(range_k),
+		};
+		thread_function_t parallelize_3d = &thread_parallelize_3d;
+		#if PTHREADPOOL_USE_FASTPATH
+			const size_t range_threshold = -threads_count;
+			if (range < range_threshold) {
+				parallelize_3d = &pthreadpool_thread_parallelize_3d_fastpath;
+			}
+		#endif
+		pthreadpool_parallelize(
+			threadpool, parallelize_3d, &params, sizeof(params),
+			task, argument, range, flags);
+	}
+}
+
+void pthreadpool_parallelize_3d_tile_1d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_3d_tile_1d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t tile_k,
+	uint32_t flags)
+{
+	size_t threads_count;
+	if (threadpool == NULL || (threads_count = threadpool->threads_count.value) <= 1 || ((range_i | range_j) <= 1 && range_k <= tile_k)) {
+		/* No thread pool used: execute task sequentially on the calling thread */
+		struct fpu_state saved_fpu_state = { 0 };
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			saved_fpu_state = get_fpu_state();
+			disable_fpu_denormals();
+		}
+		for (size_t i = 0; i < range_i; i++) {
+			for (size_t j = 0; j < range_j; j++) {
+				for (size_t k = 0; k < range_k; k += tile_k) {
+					task(argument, i, j, k, min(range_k - k, tile_k));
+				}
+			}
+		}
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			set_fpu_state(saved_fpu_state);
+		}
+	} else {
+		const size_t tile_range_k = divide_round_up(range_k, tile_k);
+		const size_t tile_range = range_i * range_j * tile_range_k;
+		const struct pthreadpool_3d_tile_1d_params params = {
+			.range_k = range_k,
+			.tile_k = tile_k,
+			.range_j = fxdiv_init_size_t(range_j),
+			.tile_range_k = fxdiv_init_size_t(tile_range_k),
+		};
+		thread_function_t parallelize_3d_tile_1d = &thread_parallelize_3d_tile_1d;
+		#if PTHREADPOOL_USE_FASTPATH
+			const size_t range_threshold = -threads_count;
+			if (tile_range < range_threshold) {
+				parallelize_3d_tile_1d = &pthreadpool_thread_parallelize_3d_tile_1d_fastpath;
+			}
+		#endif
+		pthreadpool_parallelize(
+			threadpool, parallelize_3d_tile_1d, &params, sizeof(params),
+			task, argument, tile_range, flags);
+	}
+}
+
 void pthreadpool_parallelize_3d_tile_2d(
 	pthreadpool_t threadpool,
 	pthreadpool_task_3d_tile_2d_t task,
@@ -1182,6 +1625,114 @@ void pthreadpool_parallelize_3d_tile_2d_with_uarch(
 	}
 }
 
+void pthreadpool_parallelize_4d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_4d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	uint32_t flags)
+{
+	size_t threads_count;
+	if (threadpool == NULL || (threads_count = threadpool->threads_count.value) <= 1 || (range_i | range_j | range_k | range_l) <= 1) {
+		/* No thread pool used: execute task sequentially on the calling thread */
+		struct fpu_state saved_fpu_state = { 0 };
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			saved_fpu_state = get_fpu_state();
+			disable_fpu_denormals();
+		}
+		for (size_t i = 0; i < range_i; i++) {
+			for (size_t j = 0; j < range_j; j++) {
+				for (size_t k = 0; k < range_k; k++) {
+					for (size_t l = 0; l < range_l; l++) {
+						task(argument, i, j, k, l);
+					}
+				}
+			}
+		}
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			set_fpu_state(saved_fpu_state);
+		}
+	} else {
+		const size_t range_kl = range_k * range_l;
+		const size_t range = range_i * range_j * range_kl;
+		const struct pthreadpool_4d_params params = {
+			.range_k = range_k,
+			.range_j = fxdiv_init_size_t(range_j),
+			.range_kl = fxdiv_init_size_t(range_kl),
+			.range_l = fxdiv_init_size_t(range_l),
+		};
+		thread_function_t parallelize_4d = &thread_parallelize_4d;
+		#if PTHREADPOOL_USE_FASTPATH
+			const size_t range_threshold = -threads_count;
+			if (range < range_threshold) {
+				parallelize_4d = &pthreadpool_thread_parallelize_4d_fastpath;
+			}
+		#endif
+		pthreadpool_parallelize(
+			threadpool, parallelize_4d, &params, sizeof(params),
+			task, argument, range, flags);
+	}
+}
+
+void pthreadpool_parallelize_4d_tile_1d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_4d_tile_1d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	size_t tile_l,
+	uint32_t flags)
+{
+	size_t threads_count;
+	if (threadpool == NULL || (threads_count = threadpool->threads_count.value) <= 1 || ((range_i | range_j | range_k) <= 1 && range_l <= tile_l)) {
+		/* No thread pool used: execute task sequentially on the calling thread */
+		struct fpu_state saved_fpu_state = { 0 };
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			saved_fpu_state = get_fpu_state();
+			disable_fpu_denormals();
+		}
+		for (size_t i = 0; i < range_i; i++) {
+			for (size_t j = 0; j < range_j; j++) {
+				for (size_t k = 0; k < range_k; k++) {
+					for (size_t l = 0; l < range_l; l += tile_l) {
+						task(argument, i, j, k, l, min(range_l - l, tile_l));
+					}
+				}
+			}
+		}
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			set_fpu_state(saved_fpu_state);
+		}
+	} else {
+		const size_t tile_range_l = divide_round_up(range_l, tile_l);
+		const size_t tile_range_kl = range_k * tile_range_l;
+		const size_t tile_range = range_i * range_j * tile_range_kl;
+		const struct pthreadpool_4d_tile_1d_params params = {
+			.range_k = range_k,
+			.range_l = range_l,
+			.tile_l = tile_l,
+			.range_j = fxdiv_init_size_t(range_j),
+			.tile_range_kl = fxdiv_init_size_t(tile_range_kl),
+			.tile_range_l = fxdiv_init_size_t(tile_range_l),
+		};
+		thread_function_t parallelize_4d_tile_1d = &thread_parallelize_4d_tile_1d;
+		#if PTHREADPOOL_USE_FASTPATH
+			const size_t range_threshold = -threads_count;
+			if (tile_range < range_threshold) {
+				parallelize_4d_tile_1d = &pthreadpool_thread_parallelize_4d_tile_1d_fastpath;
+			}
+		#endif
+		pthreadpool_parallelize(
+			threadpool, parallelize_4d_tile_1d, &params, sizeof(params),
+			task, argument, tile_range, flags);
+	}
+}
+
 void pthreadpool_parallelize_4d_tile_2d(
 	pthreadpool_t threadpool,
 	pthreadpool_task_4d_tile_2d_t task,
@@ -1313,6 +1864,122 @@ void pthreadpool_parallelize_4d_tile_2d_with_uarch(
 	}
 }
 
+void pthreadpool_parallelize_5d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_5d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	size_t range_m,
+	uint32_t flags)
+{
+	size_t threads_count;
+	if (threadpool == NULL || (threads_count = threadpool->threads_count.value) <= 1 || (range_i | range_j | range_k | range_l | range_m) <= 1) {
+		/* No thread pool used: execute task sequentially on the calling thread */
+		struct fpu_state saved_fpu_state = { 0 };
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			saved_fpu_state = get_fpu_state();
+			disable_fpu_denormals();
+		}
+		for (size_t i = 0; i < range_i; i++) {
+			for (size_t j = 0; j < range_j; j++) {
+				for (size_t k = 0; k < range_k; k++) {
+					for (size_t l = 0; l < range_l; l++) {
+						for (size_t m = 0; m < range_m; m++) {
+							task(argument, i, j, k, l, m);
+						}
+					}
+				}
+			}
+		}
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			set_fpu_state(saved_fpu_state);
+		}
+	} else {
+		const size_t range_lm = range_l * range_m;
+		const size_t range = range_i * range_j * range_k * range_lm;
+		const struct pthreadpool_5d_params params = {
+			.range_l = range_l,
+			.range_j = fxdiv_init_size_t(range_j),
+			.range_k = fxdiv_init_size_t(range_k),
+			.range_lm = fxdiv_init_size_t(range_lm),
+			.range_m = fxdiv_init_size_t(range_m),
+		};
+		thread_function_t parallelize_5d = &thread_parallelize_5d;
+		#if PTHREADPOOL_USE_FASTPATH
+			const size_t range_threshold = -threads_count;
+			if (range < range_threshold) {
+				parallelize_5d = &pthreadpool_thread_parallelize_5d_fastpath;
+			}
+		#endif
+		pthreadpool_parallelize(
+			threadpool, parallelize_5d, &params, sizeof(params),
+			task, argument, range, flags);
+	}
+}
+
+void pthreadpool_parallelize_5d_tile_1d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_5d_tile_1d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	size_t range_m,
+	size_t tile_m,
+	uint32_t flags)
+{
+	size_t threads_count;
+	if (threadpool == NULL || (threads_count = threadpool->threads_count.value) <= 1 || ((range_i | range_j | range_k | range_l) <= 1 && range_m <= tile_m)) {
+		/* No thread pool used: execute task sequentially on the calling thread */
+		struct fpu_state saved_fpu_state = { 0 };
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			saved_fpu_state = get_fpu_state();
+			disable_fpu_denormals();
+		}
+		for (size_t i = 0; i < range_i; i++) {
+			for (size_t j = 0; j < range_j; j++) {
+				for (size_t k = 0; k < range_k; k++) {
+					for (size_t l = 0; l < range_l; l++) {
+						for (size_t m = 0; m < range_m; m += tile_m) {
+							task(argument, i, j, k, l, m, min(range_m - m, tile_m));
+						}
+					}
+				}
+			}
+		}
+		if (flags & PTHREADPOOL_FLAG_DISABLE_DENORMALS) {
+			set_fpu_state(saved_fpu_state);
+		}
+	} else {
+		const size_t tile_range_m = divide_round_up(range_m, tile_m);
+		const size_t range_kl = range_k * range_l;
+		const size_t tile_range = range_i * range_j * range_kl * tile_range_m;
+		const struct pthreadpool_5d_tile_1d_params params = {
+			.range_k = range_k,
+			.range_m = range_m,
+			.tile_m = tile_m,
+			.range_j = fxdiv_init_size_t(range_j),
+			.range_kl = fxdiv_init_size_t(range_kl),
+			.range_l = fxdiv_init_size_t(range_l),
+			.tile_range_m = fxdiv_init_size_t(tile_range_m),
+		};
+		thread_function_t parallelize_5d_tile_1d = &thread_parallelize_5d_tile_1d;
+		#if PTHREADPOOL_USE_FASTPATH
+			const size_t range_threshold = -threads_count;
+			if (tile_range < range_threshold) {
+				parallelize_5d_tile_1d = &pthreadpool_thread_parallelize_5d_tile_1d_fastpath;
+			}
+		#endif
+		pthreadpool_parallelize(
+			threadpool, parallelize_5d_tile_1d, &params, sizeof(params),
+			task, argument, tile_range, flags);
+	}
+}
+
 void pthreadpool_parallelize_5d_tile_2d(
 	pthreadpool_t threadpool,
 	pthreadpool_task_5d_tile_2d_t task,
diff --git a/src/shim.c b/src/shim.c
index 7bf378c..e90ac45 100644
--- a/src/shim.c
+++ b/src/shim.c
@@ -133,6 +133,43 @@ void pthreadpool_parallelize_2d_tile_2d_with_uarch(
 	}
 }
 
+void pthreadpool_parallelize_3d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_3d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	uint32_t flags)
+{
+	for (size_t i = 0; i < range_i; i++) {
+		for (size_t j = 0; j < range_j; j++) {
+			for (size_t k = 0; k < range_k; k++) {
+				task(argument, i, j, k);
+			}
+		}
+	}
+}
+
+void pthreadpool_parallelize_3d_tile_1d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_3d_tile_1d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t tile_k,
+	uint32_t flags)
+{
+	for (size_t i = 0; i < range_i; i++) {
+		for (size_t j = 0; j < range_j; j++) {
+			for (size_t k = 0; k < range_k; k += tile_k) {
+				task(argument, i, j, k, min(range_k - k, tile_k));
+			}
+		}
+	}
+}
+
 void pthreadpool_parallelize_3d_tile_2d(
 	pthreadpool_t threadpool,
 	pthreadpool_task_3d_tile_2d_t task,
@@ -177,6 +214,49 @@ void pthreadpool_parallelize_3d_tile_2d_with_uarch(
 	}
 }
 
+void pthreadpool_parallelize_4d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_4d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	uint32_t flags)
+{
+	for (size_t i = 0; i < range_i; i++) {
+		for (size_t j = 0; j < range_j; j++) {
+			for (size_t k = 0; k < range_k; k++) {
+				for (size_t l = 0; l < range_l; l++) {
+					task(argument, i, j, k, l);
+				}
+			}
+		}
+	}
+}
+
+void pthreadpool_parallelize_4d_tile_1d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_4d_tile_1d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	size_t tile_l,
+	uint32_t flags)
+{
+	for (size_t i = 0; i < range_i; i++) {
+		for (size_t j = 0; j < range_j; j++) {
+			for (size_t k = 0; k < range_k; k++) {
+				for (size_t l = 0; l < range_l; l += tile_l) {
+					task(argument, i, j, k, l, min(range_l - l, tile_l));
+				}
+			}
+		}
+	}
+}
+
 void pthreadpool_parallelize_4d_tile_2d(
 	pthreadpool_t threadpool,
 	pthreadpool_task_4d_tile_2d_t task,
@@ -227,6 +307,55 @@ void pthreadpool_parallelize_4d_tile_2d_with_uarch(
 	}
 }
 
+void pthreadpool_parallelize_5d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_5d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	size_t range_m,
+	uint32_t flags)
+{
+	for (size_t i = 0; i < range_i; i++) {
+		for (size_t j = 0; j < range_j; j++) {
+			for (size_t k = 0; k < range_k; k++) {
+				for (size_t l = 0; l < range_l; l++) {
+					for (size_t m = 0; m < range_m; m++) {
+						task(argument, i, j, k, l, m);
+					}
+				}
+			}
+		}
+	}
+}
+
+void pthreadpool_parallelize_5d_tile_1d(
+	pthreadpool_t threadpool,
+	pthreadpool_task_5d_tile_1d_t task,
+	void* argument,
+	size_t range_i,
+	size_t range_j,
+	size_t range_k,
+	size_t range_l,
+	size_t range_m,
+	size_t tile_m,
+	uint32_t flags)
+{
+	for (size_t i = 0; i < range_i; i++) {
+		for (size_t j = 0; j < range_j; j++) {
+			for (size_t k = 0; k < range_k; k++) {
+				for (size_t l = 0; l < range_l; l++) {
+					for (size_t m = 0; m < range_m; m += tile_m) {
+						task(argument, i, j, k, l, m, min(range_m - m, tile_m));
+					}
+				}
+			}
+		}
+	}
+}
+
 void pthreadpool_parallelize_5d_tile_2d(
 	pthreadpool_t threadpool,
 	pthreadpool_task_5d_tile_2d_t task,
diff --git a/src/threadpool-object.h b/src/threadpool-object.h
index 7b643c6..9870e8a 100644
--- a/src/threadpool-object.h
+++ b/src/threadpool-object.h
@@ -179,6 +179,36 @@ struct pthreadpool_2d_tile_2d_with_uarch_params {
 	struct fxdiv_divisor_size_t tile_range_j;
 };
 
+struct pthreadpool_3d_params {
+	/**
+	 * FXdiv divisor for the range_j argument passed to the pthreadpool_parallelize_3d function.
+	 */
+	struct fxdiv_divisor_size_t range_j;
+	/**
+	 * FXdiv divisor for the range_k argument passed to the pthreadpool_parallelize_3d function.
+	 */
+	struct fxdiv_divisor_size_t range_k;
+};
+
+struct pthreadpool_3d_tile_1d_params {
+	/**
+	 * Copy of the range_k argument passed to the pthreadpool_parallelize_3d_tile_1d function.
+	 */
+	size_t range_k;
+	/**
+	 * Copy of the tile_k argument passed to the pthreadpool_parallelize_3d_tile_1d function.
+	 */
+	size_t tile_k;
+	/**
+	 * FXdiv divisor for the range_j argument passed to the pthreadpool_parallelize_3d_tile_1d function.
+	 */
+	struct fxdiv_divisor_size_t range_j;
+	/**
+	 * FXdiv divisor for the divide_round_up(range_k, tile_k) value.
+	 */
+	struct fxdiv_divisor_size_t tile_range_k;
+};
+
 struct pthreadpool_3d_tile_2d_params {
 	/**
 	 * Copy of the range_j argument passed to the pthreadpool_parallelize_3d_tile_2d function.
@@ -241,6 +271,52 @@ struct pthreadpool_3d_tile_2d_with_uarch_params {
 	struct fxdiv_divisor_size_t tile_range_k;
 };
 
+struct pthreadpool_4d_params {
+	/**
+	 * Copy of the range_k argument passed to the pthreadpool_parallelize_4d function.
+	 */
+	size_t range_k;
+	/**
+	 * FXdiv divisor for the range_j argument passed to the pthreadpool_parallelize_4d function.
+	 */
+	struct fxdiv_divisor_size_t range_j;
+	/**
+	 * FXdiv divisor for the range_k * range_l value.
+	 */
+	struct fxdiv_divisor_size_t range_kl;
+	/**
+	 * FXdiv divisor for the range_l argument passed to the pthreadpool_parallelize_4d function.
+	 */
+	struct fxdiv_divisor_size_t range_l;
+};
+
+struct pthreadpool_4d_tile_1d_params {
+	/**
+	 * Copy of the range_k argument passed to the pthreadpool_parallelize_4d_tile_1d function.
+	 */
+	size_t range_k;
+	/**
+	 * Copy of the range_l argument passed to the pthreadpool_parallelize_4d_tile_1d function.
+	 */
+	size_t range_l;
+	/**
+	 * Copy of the tile_l argument passed to the pthreadpool_parallelize_4d_tile_1d function.
+	 */
+	size_t tile_l;
+	/**
+	 * FXdiv divisor for the range_j argument passed to the pthreadpool_parallelize_4d_tile_1d function.
+	 */
+	struct fxdiv_divisor_size_t range_j;
+	/**
+	 * FXdiv divisor for the range_k * divide_round_up(range_l, tile_l) value.
+	 */
+	struct fxdiv_divisor_size_t tile_range_kl;
+	/**
+	 * FXdiv divisor for the divide_round_up(range_l, tile_l) value.
+	 */
+	struct fxdiv_divisor_size_t tile_range_l;
+};
+
 struct pthreadpool_4d_tile_2d_params {
 	/**
 	 * Copy of the range_k argument passed to the pthreadpool_parallelize_4d_tile_2d function.
@@ -311,6 +387,60 @@ struct pthreadpool_4d_tile_2d_with_uarch_params {
 	struct fxdiv_divisor_size_t tile_range_l;
 };
 
+struct pthreadpool_5d_params {
+	/**
+	 * Copy of the range_l argument passed to the pthreadpool_parallelize_5d function.
+	 */
+	size_t range_l;
+	/**
+	 * FXdiv divisor for the range_j argument passed to the pthreadpool_parallelize_5d function.
+	 */
+	struct fxdiv_divisor_size_t range_j;
+	/**
+	 * FXdiv divisor for the range_k argument passed to the pthreadpool_parallelize_5d function.
+	 */
+	struct fxdiv_divisor_size_t range_k;
+	/**
+	 * FXdiv divisor for the range_l * range_m value.
+	 */
+	struct fxdiv_divisor_size_t range_lm;
+	/**
+	 * FXdiv divisor for the range_m argument passed to the pthreadpool_parallelize_5d function.
+	 */
+	struct fxdiv_divisor_size_t range_m;
+};
+
+struct pthreadpool_5d_tile_1d_params {
+	/**
+	 * Copy of the range_k argument passed to the pthreadpool_parallelize_5d_tile_1d function.
+	 */
+	size_t range_k;
+	/**
+	 * Copy of the range_m argument passed to the pthreadpool_parallelize_5d_tile_1d function.
+	 */
+	size_t range_m;
+	/**
+	 * Copy of the tile_m argument passed to the pthreadpool_parallelize_5d_tile_1d function.
+	 */
+	size_t tile_m;
+	/**
+	 * FXdiv divisor for the range_j argument passed to the pthreadpool_parallelize_5d_tile_1d function.
+	 */
+	struct fxdiv_divisor_size_t range_j;
+	/**
+	 * FXdiv divisor for the range_k * range_l value.
+	 */
+	struct fxdiv_divisor_size_t range_kl;
+	/**
+	 * FXdiv divisor for the range_l argument passed to the pthreadpool_parallelize_5d_tile_1d function.
+	 */
+	struct fxdiv_divisor_size_t range_l;
+	/**
+	 * FXdiv divisor for the divide_round_up(range_m, tile_m) value.
+	 */
+	struct fxdiv_divisor_size_t tile_range_m;
+};
+
 struct pthreadpool_5d_tile_2d_params {
 	/**
 	 * Copy of the range_l argument passed to the pthreadpool_parallelize_5d_tile_2d function.
@@ -434,10 +564,16 @@ struct PTHREADPOOL_CACHELINE_ALIGNED pthreadpool {
 		struct pthreadpool_2d_tile_1d_params parallelize_2d_tile_1d;
 		struct pthreadpool_2d_tile_2d_params parallelize_2d_tile_2d;
 		struct pthreadpool_2d_tile_2d_with_uarch_params parallelize_2d_tile_2d_with_uarch;
+		struct pthreadpool_3d_params parallelize_3d;
+		struct pthreadpool_3d_tile_1d_params parallelize_3d_tile_1d;
 		struct pthreadpool_3d_tile_2d_params parallelize_3d_tile_2d;
 		struct pthreadpool_3d_tile_2d_with_uarch_params parallelize_3d_tile_2d_with_uarch;
+		struct pthreadpool_4d_params parallelize_4d;
+		struct pthreadpool_4d_tile_1d_params parallelize_4d_tile_1d;
 		struct pthreadpool_4d_tile_2d_params parallelize_4d_tile_2d;
 		struct pthreadpool_4d_tile_2d_with_uarch_params parallelize_4d_tile_2d_with_uarch;
+		struct pthreadpool_5d_params parallelize_5d;
+		struct pthreadpool_5d_tile_1d_params parallelize_5d_tile_1d;
 		struct pthreadpool_5d_tile_2d_params parallelize_5d_tile_2d;
 		struct pthreadpool_6d_tile_2d_params parallelize_6d_tile_2d;
 	} params;
@@ -555,6 +691,14 @@ PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_2d_tile_2d_with_uarch_f
 	struct pthreadpool* threadpool,
 	struct thread_info* thread);
 
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_3d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread);
+
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_3d_tile_1d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread);
+
 PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_3d_tile_2d_fastpath(
 	struct pthreadpool* threadpool,
 	struct thread_info* thread);
@@ -563,6 +707,14 @@ PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_3d_tile_2d_with_uarch_f
 	struct pthreadpool* threadpool,
 	struct thread_info* thread);
 
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_4d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread);
+
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_4d_tile_1d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread);
+
 PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_4d_tile_2d_fastpath(
 	struct pthreadpool* threadpool,
 	struct thread_info* thread);
@@ -571,6 +723,14 @@ PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_4d_tile_2d_with_uarch_f
 	struct pthreadpool* threadpool,
 	struct thread_info* thread);
 
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_5d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread);
+
+PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_5d_tile_1d_fastpath(
+	struct pthreadpool* threadpool,
+	struct thread_info* thread);
+
 PTHREADPOOL_INTERNAL void pthreadpool_thread_parallelize_5d_tile_2d_fastpath(
 	struct pthreadpool* threadpool,
 	struct thread_info* thread);
diff --git a/test/pthreadpool.cc b/test/pthreadpool.cc
index b8a6803..f822506 100644
--- a/test/pthreadpool.cc
+++ b/test/pthreadpool.cc
@@ -23,17 +23,44 @@ const size_t kParallelize2DTile2DRangeI = 53;
 const size_t kParallelize2DTile2DRangeJ = 59;
 const size_t kParallelize2DTile2DTileI = 5;
 const size_t kParallelize2DTile2DTileJ = 7;
+const size_t kParallelize3DRangeI = 13;
+const size_t kParallelize3DRangeJ = 17;
+const size_t kParallelize3DRangeK = 19;
+const size_t kParallelize3DTile1DRangeI = 17;
+const size_t kParallelize3DTile1DRangeJ = 19;
+const size_t kParallelize3DTile1DRangeK = 23;
+const size_t kParallelize3DTile1DTileK = 5;
 const size_t kParallelize3DTile2DRangeI = 19;
 const size_t kParallelize3DTile2DRangeJ = 23;
 const size_t kParallelize3DTile2DRangeK = 29;
 const size_t kParallelize3DTile2DTileJ = 2;
 const size_t kParallelize3DTile2DTileK = 3;
+const size_t kParallelize4DRangeI = 11;
+const size_t kParallelize4DRangeJ = 13;
+const size_t kParallelize4DRangeK = 17;
+const size_t kParallelize4DRangeL = 19;
+const size_t kParallelize4DTile1DRangeI = 13;
+const size_t kParallelize4DTile1DRangeJ = 17;
+const size_t kParallelize4DTile1DRangeK = 19;
+const size_t kParallelize4DTile1DRangeL = 23;
+const size_t kParallelize4DTile1DTileL = 5;
 const size_t kParallelize4DTile2DRangeI = 17;
 const size_t kParallelize4DTile2DRangeJ = 19;
 const size_t kParallelize4DTile2DRangeK = 23;
 const size_t kParallelize4DTile2DRangeL = 29;
 const size_t kParallelize4DTile2DTileK = 2;
 const size_t kParallelize4DTile2DTileL = 3;
+const size_t kParallelize5DRangeI = 7;
+const size_t kParallelize5DRangeJ = 11;
+const size_t kParallelize5DRangeK = 13;
+const size_t kParallelize5DRangeL = 17;
+const size_t kParallelize5DRangeM = 19;
+const size_t kParallelize5DTile1DRangeI = 11;
+const size_t kParallelize5DTile1DRangeJ = 13;
+const size_t kParallelize5DTile1DRangeK = 17;
+const size_t kParallelize5DTile1DRangeL = 19;
+const size_t kParallelize5DTile1DRangeM = 23;
+const size_t kParallelize5DTile1DTileM = 5;
 const size_t kParallelize5DTile2DRangeI = 13;
 const size_t kParallelize5DTile2DRangeJ = 17;
 const size_t kParallelize5DTile2DRangeK = 19;
@@ -2286,6 +2313,646 @@ TEST(Parallelize2DTile2DWithUArch, MultiThreadPoolWorkStealing) {
 	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize2DTile2DRangeI * kParallelize2DTile2DRangeJ);
 }
 
+static void ComputeNothing3D(void*, size_t, size_t, size_t) {
+}
+
+TEST(Parallelize3D, SingleThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_3d(threadpool.get(),
+		ComputeNothing3D,
+		nullptr,
+		kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+		0 /* flags */);
+}
+
+TEST(Parallelize3D, MultiThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d(
+		threadpool.get(),
+		ComputeNothing3D,
+		nullptr,
+		kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+		0 /* flags */);
+}
+
+static void CheckBounds3D(void*, size_t i, size_t j, size_t k) {
+	EXPECT_LT(i, kParallelize3DRangeI);
+	EXPECT_LT(j, kParallelize3DRangeJ);
+	EXPECT_LT(k, kParallelize3DRangeK);
+}
+
+TEST(Parallelize3D, SingleThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_3d(
+		threadpool.get(),
+		CheckBounds3D,
+		nullptr,
+		kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+		0 /* flags */);
+}
+
+TEST(Parallelize3D, MultiThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d(
+		threadpool.get(),
+		CheckBounds3D,
+		nullptr,
+		kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+		0 /* flags */);
+}
+
+static void SetTrue3D(std::atomic_bool* processed_indicators, size_t i, size_t j, size_t k) {
+	const size_t linear_idx = (i * kParallelize3DRangeJ + j) * kParallelize3DRangeK + k;
+	processed_indicators[linear_idx].store(true, std::memory_order_relaxed);
+}
+
+TEST(Parallelize3D, SingleThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize3DRangeI * kParallelize3DRangeJ * kParallelize3DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_3d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_t>(SetTrue3D),
+		static_cast<void*>(indicators.data()),
+		kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize3DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DRangeJ + j) * kParallelize3DRangeK + k;
+				EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+					<< "Element (" << i << ", " << j << ", " << k << ") not processed";
+			}
+		}
+	}
+}
+
+TEST(Parallelize3D, MultiThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize3DRangeI * kParallelize3DRangeJ * kParallelize3DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_t>(SetTrue3D),
+		static_cast<void*>(indicators.data()),
+		kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize3DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DRangeJ + j) * kParallelize3DRangeK + k;
+				EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+					<< "Element (" << i << ", " << j << ", " << k << ") not processed";
+			}
+		}
+	}
+}
+
+static void Increment3D(std::atomic_int* processed_counters, size_t i, size_t j, size_t k) {
+	const size_t linear_idx = (i * kParallelize3DRangeJ + j) * kParallelize3DRangeK + k;
+	processed_counters[linear_idx].fetch_add(1, std::memory_order_relaxed);
+}
+
+TEST(Parallelize3D, SingleThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize3DRangeI * kParallelize3DRangeJ * kParallelize3DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_3d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_t>(Increment3D),
+		static_cast<void*>(counters.data()),
+		kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize3DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DRangeJ + j) * kParallelize3DRangeK + k;
+				EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+					<< "Element (" << i << ", " << j << ", " << k << ") was processed "
+					<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+			}
+		}
+	}
+}
+
+TEST(Parallelize3D, MultiThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize3DRangeI * kParallelize3DRangeJ * kParallelize3DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_t>(Increment3D),
+		static_cast<void*>(counters.data()),
+		kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize3DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DRangeJ + j) * kParallelize3DRangeK + k;
+				EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+					<< "Element (" << i << ", " << j << ", " << k << ") was processed "
+					<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+			}
+		}
+	}
+}
+
+TEST(Parallelize3D, SingleThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize3DRangeI * kParallelize3DRangeJ * kParallelize3DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	for (size_t iteration = 0; iteration < kIncrementIterations; iteration++) {
+		pthreadpool_parallelize_3d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_3d_t>(Increment3D),
+			static_cast<void*>(counters.data()),
+			kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+				0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize3DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DRangeJ + j) * kParallelize3DRangeK + k;
+				EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations)
+					<< "Element (" << i << ", " << j << ", " << k << ") was processed "
+					<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+					<< "(expected: " << kIncrementIterations << ")";
+			}
+		}
+	}
+}
+
+TEST(Parallelize3D, MultiThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize3DRangeI * kParallelize3DRangeJ * kParallelize3DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	for (size_t iteration = 0; iteration < kIncrementIterations; iteration++) {
+		pthreadpool_parallelize_3d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_3d_t>(Increment3D),
+			static_cast<void*>(counters.data()),
+			kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+				0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize3DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DRangeJ + j) * kParallelize3DRangeK + k;
+				EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations)
+					<< "Element (" << i << ", " << j << ", " << k << ") was processed "
+					<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+					<< "(expected: " << kIncrementIterations << ")";
+			}
+		}
+	}
+}
+
+static void IncrementSame3D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t k) {
+	num_processed_items->fetch_add(1, std::memory_order_relaxed);
+}
+
+TEST(Parallelize3D, MultiThreadPoolHighContention) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_t>(IncrementSame3D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize3DRangeI * kParallelize3DRangeJ * kParallelize3DRangeK);
+}
+
+static void WorkImbalance3D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t k) {
+	num_processed_items->fetch_add(1, std::memory_order_relaxed);
+	if (i == 0 && j == 0 && k == 0) {
+		/* Spin-wait until all items are computed */
+		while (num_processed_items->load(std::memory_order_relaxed) != kParallelize3DRangeI * kParallelize3DRangeJ * kParallelize3DRangeK) {
+			std::atomic_thread_fence(std::memory_order_acquire);
+		}
+	}
+}
+
+TEST(Parallelize3D, MultiThreadPoolWorkStealing) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_t>(WorkImbalance3D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize3DRangeI, kParallelize3DRangeJ, kParallelize3DRangeK,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize3DRangeI * kParallelize3DRangeJ * kParallelize3DRangeK);
+}
+
+static void ComputeNothing3DTile1D(void*, size_t, size_t, size_t, size_t) {
+}
+
+TEST(Parallelize3DTile1D, SingleThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_3d_tile_1d(threadpool.get(),
+		ComputeNothing3DTile1D,
+		nullptr,
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+}
+
+TEST(Parallelize3DTile1D, MultiThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		ComputeNothing3DTile1D,
+		nullptr,
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+}
+
+static void CheckBounds3DTile1D(void*, size_t i, size_t j, size_t start_k, size_t tile_k) {
+	EXPECT_LT(i, kParallelize3DTile1DRangeI);
+	EXPECT_LT(j, kParallelize3DTile1DRangeJ);
+	EXPECT_LT(start_k, kParallelize3DTile1DRangeK);
+	EXPECT_LE(start_k + tile_k, kParallelize3DTile1DRangeK);
+}
+
+TEST(Parallelize3DTile1D, SingleThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		CheckBounds3DTile1D,
+		nullptr,
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+}
+
+TEST(Parallelize3DTile1D, MultiThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		CheckBounds3DTile1D,
+		nullptr,
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+}
+
+static void CheckTiling3DTile1D(void*, size_t i, size_t j, size_t start_k, size_t tile_k) {
+	EXPECT_GT(tile_k, 0);
+	EXPECT_LE(tile_k, kParallelize3DTile1DTileK);
+	EXPECT_EQ(start_k % kParallelize3DTile1DTileK, 0);
+	EXPECT_EQ(tile_k, std::min<size_t>(kParallelize3DTile1DTileK, kParallelize3DTile1DRangeK - start_k));
+}
+
+TEST(Parallelize3DTile1D, SingleThreadPoolUniformTiling) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		CheckTiling3DTile1D,
+		nullptr,
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+}
+
+TEST(Parallelize3DTile1D, MultiThreadPoolUniformTiling) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		CheckTiling3DTile1D,
+		nullptr,
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+}
+
+static void SetTrue3DTile1D(std::atomic_bool* processed_indicators, size_t i, size_t j, size_t start_k, size_t tile_k) {
+	for (size_t k = start_k; k < start_k + tile_k; k++) {
+		const size_t linear_idx = (i * kParallelize3DTile1DRangeJ + j) * kParallelize3DTile1DRangeK + k;
+		processed_indicators[linear_idx].store(true, std::memory_order_relaxed);
+	}
+}
+
+TEST(Parallelize3DTile1D, SingleThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize3DTile1DRangeI * kParallelize3DTile1DRangeJ * kParallelize3DTile1DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_tile_1d_t>(SetTrue3DTile1D),
+		static_cast<void*>(indicators.data()),
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize3DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DTile1DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DTile1DRangeJ + j) * kParallelize3DTile1DRangeK + k;
+				EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+					<< "Element (" << i << ", " << j << ", " << k << ") not processed";
+			}
+		}
+	}
+}
+
+TEST(Parallelize3DTile1D, MultiThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize3DTile1DRangeI * kParallelize3DTile1DRangeJ * kParallelize3DTile1DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_tile_1d_t>(SetTrue3DTile1D),
+		static_cast<void*>(indicators.data()),
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize3DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DTile1DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DTile1DRangeJ + j) * kParallelize3DTile1DRangeK + k;
+				EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+					<< "Element (" << i << ", " << j << ", " << k << ") not processed";
+			}
+		}
+	}
+}
+
+static void Increment3DTile1D(std::atomic_int* processed_counters, size_t i, size_t j, size_t start_k, size_t tile_k) {
+	for (size_t k = start_k; k < start_k + tile_k; k++) {
+		const size_t linear_idx = (i * kParallelize3DTile1DRangeJ + j) * kParallelize3DTile1DRangeK + k;
+		processed_counters[linear_idx].fetch_add(1, std::memory_order_relaxed);
+	}
+}
+
+TEST(Parallelize3DTile1D, SingleThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize3DTile1DRangeI * kParallelize3DTile1DRangeJ * kParallelize3DTile1DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_tile_1d_t>(Increment3DTile1D),
+		static_cast<void*>(counters.data()),
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize3DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DTile1DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DTile1DRangeJ + j) * kParallelize3DTile1DRangeK + k;
+				EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+					<< "Element (" << i << ", " << j << ", " << k << ") was processed "
+					<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+			}
+		}
+	}
+}
+
+TEST(Parallelize3DTile1D, MultiThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize3DTile1DRangeI * kParallelize3DTile1DRangeJ * kParallelize3DTile1DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_tile_1d_t>(Increment3DTile1D),
+		static_cast<void*>(counters.data()),
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize3DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DTile1DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DTile1DRangeJ + j) * kParallelize3DTile1DRangeK + k;
+				EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+					<< "Element (" << i << ", " << j << ", " << k << ") was processed "
+					<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+			}
+		}
+	}
+}
+
+TEST(Parallelize3DTile1D, SingleThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize3DTile1DRangeI * kParallelize3DTile1DRangeJ * kParallelize3DTile1DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	for (size_t iteration = 0; iteration < kIncrementIterations; iteration++) {
+		pthreadpool_parallelize_3d_tile_1d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_3d_tile_1d_t>(Increment3DTile1D),
+			static_cast<void*>(counters.data()),
+			kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+			kParallelize3DTile1DTileK,
+			0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize3DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DTile1DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DTile1DRangeJ + j) * kParallelize3DTile1DRangeK + k;
+				EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations)
+					<< "Element (" << i << ", " << j << ", " << k << ") was processed "
+					<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+					<< "(expected: " << kIncrementIterations << ")";
+			}
+		}
+	}
+}
+
+TEST(Parallelize3DTile1D, MultiThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize3DTile1DRangeI * kParallelize3DTile1DRangeJ * kParallelize3DTile1DRangeK);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	for (size_t iteration = 0; iteration < kIncrementIterations; iteration++) {
+		pthreadpool_parallelize_3d_tile_1d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_3d_tile_1d_t>(Increment3DTile1D),
+			static_cast<void*>(counters.data()),
+			kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+			kParallelize3DTile1DTileK,
+			0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize3DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize3DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize3DTile1DRangeK; k++) {
+				const size_t linear_idx = (i * kParallelize3DTile1DRangeJ + j) * kParallelize3DTile1DRangeK + k;
+				EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations)
+					<< "Element (" << i << ", " << j << ", " << k << ") was processed "
+					<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+					<< "(expected: " << kIncrementIterations << ")";
+			}
+		}
+	}
+}
+
+static void IncrementSame3DTile1D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t start_k, size_t tile_k) {
+	for (size_t k = start_k; k < start_k + tile_k; k++) {
+		num_processed_items->fetch_add(1, std::memory_order_relaxed);
+	}
+}
+
+TEST(Parallelize3DTile1D, MultiThreadPoolHighContention) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_tile_1d_t>(IncrementSame3DTile1D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize3DTile1DRangeI * kParallelize3DTile1DRangeJ * kParallelize3DTile1DRangeK);
+}
+
+static void WorkImbalance3DTile1D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t start_k, size_t tile_k) {
+	num_processed_items->fetch_add(tile_k, std::memory_order_relaxed);
+	if (i == 0 && j == 0 && start_k == 0) {
+		/* Spin-wait until all items are computed */
+		while (num_processed_items->load(std::memory_order_relaxed) != kParallelize3DTile1DRangeI * kParallelize3DTile1DRangeJ * kParallelize3DTile1DRangeK) {
+			std::atomic_thread_fence(std::memory_order_acquire);
+		}
+	}
+}
+
+TEST(Parallelize3DTile1D, MultiThreadPoolWorkStealing) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_3d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_3d_tile_1d_t>(WorkImbalance3DTile1D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize3DTile1DRangeI, kParallelize3DTile1DRangeJ, kParallelize3DTile1DRangeK,
+		kParallelize3DTile1DTileK,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize3DTile1DRangeI * kParallelize3DTile1DRangeJ * kParallelize3DTile1DRangeK);
+}
+
 static void ComputeNothing3DTile2D(void*, size_t, size_t, size_t, size_t, size_t) {
 }
 
@@ -3058,6 +3725,672 @@ TEST(Parallelize3DTile2DWithUArch, MultiThreadPoolWorkStealing) {
 	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize3DTile2DRangeI * kParallelize3DTile2DRangeJ * kParallelize3DTile2DRangeK);
 }
 
+static void ComputeNothing4D(void*, size_t, size_t, size_t, size_t) {
+}
+
+TEST(Parallelize4D, SingleThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_4d(threadpool.get(),
+		ComputeNothing4D,
+		nullptr,
+		kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+		0 /* flags */);
+}
+
+TEST(Parallelize4D, MultiThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d(
+		threadpool.get(),
+		ComputeNothing4D,
+		nullptr,
+		kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+		0 /* flags */);
+}
+
+static void CheckBounds4D(void*, size_t i, size_t j, size_t k, size_t l) {
+	EXPECT_LT(i, kParallelize4DRangeI);
+	EXPECT_LT(j, kParallelize4DRangeJ);
+	EXPECT_LT(k, kParallelize4DRangeK);
+	EXPECT_LT(l, kParallelize4DRangeL);
+}
+
+TEST(Parallelize4D, SingleThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_4d(
+		threadpool.get(),
+		CheckBounds4D,
+		nullptr,
+		kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+		0 /* flags */);
+}
+
+TEST(Parallelize4D, MultiThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d(
+		threadpool.get(),
+		CheckBounds4D,
+		nullptr,
+		kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+		0 /* flags */);
+}
+
+static void SetTrue4D(std::atomic_bool* processed_indicators, size_t i, size_t j, size_t k, size_t l) {
+	const size_t linear_idx = ((i * kParallelize4DRangeJ + j) * kParallelize4DRangeK + k) * kParallelize4DRangeL + l;
+	processed_indicators[linear_idx].store(true, std::memory_order_relaxed);
+}
+
+TEST(Parallelize4D, SingleThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize4DRangeI * kParallelize4DRangeJ * kParallelize4DRangeK * kParallelize4DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_4d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_t>(SetTrue4D),
+		static_cast<void*>(indicators.data()),
+		kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize4DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DRangeJ + j) * kParallelize4DRangeK + k) * kParallelize4DRangeL + l;
+					EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") not processed";
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize4D, MultiThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize4DRangeI * kParallelize4DRangeJ * kParallelize4DRangeK * kParallelize4DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_t>(SetTrue4D),
+		static_cast<void*>(indicators.data()),
+		kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize4DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DRangeJ + j) * kParallelize4DRangeK + k) * kParallelize4DRangeL + l;
+					EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") not processed";
+				}
+			}
+		}
+	}
+}
+
+static void Increment4D(std::atomic_int* processed_counters, size_t i, size_t j, size_t k, size_t l) {
+	const size_t linear_idx = ((i * kParallelize4DRangeJ + j) * kParallelize4DRangeK + k) * kParallelize4DRangeL + l;
+	processed_counters[linear_idx].fetch_add(1, std::memory_order_relaxed);
+}
+
+TEST(Parallelize4D, SingleThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize4DRangeI * kParallelize4DRangeJ * kParallelize4DRangeK * kParallelize4DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_4d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_t>(Increment4D),
+		static_cast<void*>(counters.data()),
+		kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize4DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DRangeJ + j) * kParallelize4DRangeK + k) * kParallelize4DRangeL + l;
+					EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") was processed "
+						<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize4D, MultiThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize4DRangeI * kParallelize4DRangeJ * kParallelize4DRangeK * kParallelize4DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_t>(Increment4D),
+		static_cast<void*>(counters.data()),
+		kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize4DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DRangeJ + j) * kParallelize4DRangeK + k) * kParallelize4DRangeL + l;
+					EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") was processed "
+						<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize4D, SingleThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize4DRangeI * kParallelize4DRangeJ * kParallelize4DRangeK * kParallelize4DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	for (size_t iteration = 0; iteration < kIncrementIterations; iteration++) {
+		pthreadpool_parallelize_4d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_4d_t>(Increment4D),
+			static_cast<void*>(counters.data()),
+			kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+				0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize4DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DRangeJ + j) * kParallelize4DRangeK + k) * kParallelize4DRangeL + l;
+					EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations)
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") was processed "
+						<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+						<< "(expected: " << kIncrementIterations << ")";
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize4D, MultiThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize4DRangeI * kParallelize4DRangeJ * kParallelize4DRangeK * kParallelize4DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	for (size_t iteration = 0; iteration < kIncrementIterations; iteration++) {
+		pthreadpool_parallelize_4d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_4d_t>(Increment4D),
+			static_cast<void*>(counters.data()),
+			kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+				0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize4DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DRangeJ + j) * kParallelize4DRangeK + k) * kParallelize4DRangeL + l;
+					EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations)
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") was processed "
+						<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+						<< "(expected: " << kIncrementIterations << ")";
+				}
+			}
+		}
+	}
+}
+
+static void IncrementSame4D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t k, size_t l) {
+	num_processed_items->fetch_add(1, std::memory_order_relaxed);
+}
+
+TEST(Parallelize4D, MultiThreadPoolHighContention) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_t>(IncrementSame4D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize4DRangeI * kParallelize4DRangeJ * kParallelize4DRangeK * kParallelize4DRangeL);
+}
+
+static void WorkImbalance4D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t k, size_t l) {
+	num_processed_items->fetch_add(1, std::memory_order_relaxed);
+	if (i == 0 && j == 0 && k == 0 && l == 0) {
+		/* Spin-wait until all items are computed */
+		while (num_processed_items->load(std::memory_order_relaxed) != kParallelize4DRangeI * kParallelize4DRangeJ * kParallelize4DRangeK * kParallelize4DRangeL) {
+			std::atomic_thread_fence(std::memory_order_acquire);
+		}
+	}
+}
+
+TEST(Parallelize4D, MultiThreadPoolWorkStealing) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_t>(WorkImbalance4D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize4DRangeI, kParallelize4DRangeJ, kParallelize4DRangeK, kParallelize4DRangeL,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize4DRangeI * kParallelize4DRangeJ * kParallelize4DRangeK * kParallelize4DRangeL);
+}
+
+static void ComputeNothing4DTile1D(void*, size_t, size_t, size_t, size_t, size_t) {
+}
+
+TEST(Parallelize4DTile1D, SingleThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_4d_tile_1d(threadpool.get(),
+		ComputeNothing4DTile1D,
+		nullptr,
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+}
+
+TEST(Parallelize4DTile1D, MultiThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		ComputeNothing4DTile1D,
+		nullptr,
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+}
+
+static void CheckBounds4DTile1D(void*, size_t i, size_t j, size_t k, size_t start_l, size_t tile_l) {
+	EXPECT_LT(i, kParallelize4DTile1DRangeI);
+	EXPECT_LT(j, kParallelize4DTile1DRangeJ);
+	EXPECT_LT(k, kParallelize4DTile1DRangeK);
+	EXPECT_LT(start_l, kParallelize4DTile1DRangeL);
+	EXPECT_LE(start_l + tile_l, kParallelize4DTile1DRangeL);
+}
+
+TEST(Parallelize4DTile1D, SingleThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		CheckBounds4DTile1D,
+		nullptr,
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+}
+
+TEST(Parallelize4DTile1D, MultiThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		CheckBounds4DTile1D,
+		nullptr,
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+}
+
+static void CheckTiling4DTile1D(void*, size_t i, size_t j, size_t k, size_t start_l, size_t tile_l) {
+	EXPECT_GT(tile_l, 0);
+	EXPECT_LE(tile_l, kParallelize4DTile1DTileL);
+	EXPECT_EQ(start_l % kParallelize4DTile1DTileL, 0);
+	EXPECT_EQ(tile_l, std::min<size_t>(kParallelize4DTile1DTileL, kParallelize4DTile1DRangeL - start_l));
+}
+
+TEST(Parallelize4DTile1D, SingleThreadPoolUniformTiling) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		CheckTiling4DTile1D,
+		nullptr,
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+}
+
+TEST(Parallelize4DTile1D, MultiThreadPoolUniformTiling) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		CheckTiling4DTile1D,
+		nullptr,
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+}
+
+static void SetTrue4DTile1D(std::atomic_bool* processed_indicators, size_t i, size_t j, size_t k, size_t start_l, size_t tile_l) {
+	for (size_t l = start_l; l < start_l + tile_l; l++) {
+		const size_t linear_idx = ((i * kParallelize4DTile1DRangeJ + j) * kParallelize4DTile1DRangeK + k) * kParallelize4DTile1DRangeL + l;
+		processed_indicators[linear_idx].store(true, std::memory_order_relaxed);
+	}
+}
+
+TEST(Parallelize4DTile1D, SingleThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize4DTile1DRangeI * kParallelize4DTile1DRangeJ * kParallelize4DTile1DRangeK * kParallelize4DTile1DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_tile_1d_t>(SetTrue4DTile1D),
+		static_cast<void*>(indicators.data()),
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize4DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DTile1DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DTile1DRangeJ + j) * kParallelize4DTile1DRangeK + k) * kParallelize4DTile1DRangeL + l;
+					EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") not processed";
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize4DTile1D, MultiThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize4DTile1DRangeI * kParallelize4DTile1DRangeJ * kParallelize4DTile1DRangeK * kParallelize4DTile1DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_tile_1d_t>(SetTrue4DTile1D),
+		static_cast<void*>(indicators.data()),
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize4DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DTile1DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DTile1DRangeJ + j) * kParallelize4DTile1DRangeK + k) * kParallelize4DTile1DRangeL + l;
+					EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") not processed";
+				}
+			}
+		}
+	}
+}
+
+static void Increment4DTile1D(std::atomic_int* processed_counters, size_t i, size_t j, size_t k, size_t start_l, size_t tile_l) {
+	for (size_t l = start_l; l < start_l + tile_l; l++) {
+		const size_t linear_idx = ((i * kParallelize4DTile1DRangeJ + j) * kParallelize4DTile1DRangeK + k) * kParallelize4DTile1DRangeL + l;
+		processed_counters[linear_idx].fetch_add(1, std::memory_order_relaxed);
+	}
+}
+
+TEST(Parallelize4DTile1D, SingleThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize4DTile1DRangeI * kParallelize4DTile1DRangeJ * kParallelize4DTile1DRangeK * kParallelize4DTile1DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_tile_1d_t>(Increment4DTile1D),
+		static_cast<void*>(counters.data()),
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize4DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DTile1DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DTile1DRangeJ + j) * kParallelize4DTile1DRangeK + k) * kParallelize4DTile1DRangeL + l;
+					EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") was processed "
+						<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize4DTile1D, MultiThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize4DTile1DRangeI * kParallelize4DTile1DRangeJ * kParallelize4DTile1DRangeK * kParallelize4DTile1DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_tile_1d_t>(Increment4DTile1D),
+		static_cast<void*>(counters.data()),
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize4DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DTile1DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DTile1DRangeJ + j) * kParallelize4DTile1DRangeK + k) * kParallelize4DTile1DRangeL + l;
+					EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") was processed "
+						<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize4DTile1D, SingleThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize4DTile1DRangeI * kParallelize4DTile1DRangeJ * kParallelize4DTile1DRangeK * kParallelize4DTile1DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	for (size_t iteration = 0; iteration < kIncrementIterations; iteration++) {
+		pthreadpool_parallelize_4d_tile_1d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_4d_tile_1d_t>(Increment4DTile1D),
+			static_cast<void*>(counters.data()),
+			kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+			kParallelize4DTile1DTileL,
+			0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize4DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DTile1DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DTile1DRangeJ + j) * kParallelize4DTile1DRangeK + k) * kParallelize4DTile1DRangeL + l;
+					EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations)
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") was processed "
+						<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+						<< "(expected: " << kIncrementIterations << ")";
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize4DTile1D, MultiThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize4DTile1DRangeI * kParallelize4DTile1DRangeJ * kParallelize4DTile1DRangeK * kParallelize4DTile1DRangeL);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	for (size_t iteration = 0; iteration < kIncrementIterations; iteration++) {
+		pthreadpool_parallelize_4d_tile_1d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_4d_tile_1d_t>(Increment4DTile1D),
+			static_cast<void*>(counters.data()),
+			kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+			kParallelize4DTile1DTileL,
+			0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize4DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize4DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize4DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize4DTile1DRangeL; l++) {
+					const size_t linear_idx = ((i * kParallelize4DTile1DRangeJ + j) * kParallelize4DTile1DRangeK + k) * kParallelize4DTile1DRangeL + l;
+					EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations)
+						<< "Element (" << i << ", " << j << ", " << k << ", " << l << ") was processed "
+						<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+						<< "(expected: " << kIncrementIterations << ")";
+				}
+			}
+		}
+	}
+}
+
+static void IncrementSame4DTile1D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t k, size_t start_l, size_t tile_l) {
+	for (size_t l = start_l; l < start_l + tile_l; l++) {
+		num_processed_items->fetch_add(1, std::memory_order_relaxed);
+	}
+}
+
+TEST(Parallelize4DTile1D, MultiThreadPoolHighContention) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_tile_1d_t>(IncrementSame4DTile1D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize4DTile1DRangeI * kParallelize4DTile1DRangeJ * kParallelize4DTile1DRangeK * kParallelize4DTile1DRangeL);
+}
+
+static void WorkImbalance4DTile1D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t k, size_t start_l, size_t tile_l) {
+	num_processed_items->fetch_add(tile_l, std::memory_order_relaxed);
+	if (i == 0 && j == 0 && k == 0 && start_l == 0) {
+		/* Spin-wait until all items are computed */
+		while (num_processed_items->load(std::memory_order_relaxed) != kParallelize4DTile1DRangeI * kParallelize4DTile1DRangeJ * kParallelize4DTile1DRangeK * kParallelize4DTile1DRangeL) {
+			std::atomic_thread_fence(std::memory_order_acquire);
+		}
+	}
+}
+
+TEST(Parallelize4DTile1D, MultiThreadPoolWorkStealing) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_4d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_4d_tile_1d_t>(WorkImbalance4DTile1D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize4DTile1DRangeI, kParallelize4DTile1DRangeJ, kParallelize4DTile1DRangeK, kParallelize4DTile1DRangeL,
+		kParallelize4DTile1DTileL,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize4DTile1DRangeI * kParallelize4DTile1DRangeJ * kParallelize4DTile1DRangeK * kParallelize4DTile1DRangeL);
+}
+
 static void ComputeNothing4DTile2D(void*, size_t, size_t, size_t, size_t, size_t, size_t) {
 }
 
@@ -3856,6 +5189,698 @@ TEST(Parallelize4DTile2DWithUArch, MultiThreadPoolWorkStealing) {
 	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize4DTile2DRangeI * kParallelize4DTile2DRangeJ * kParallelize4DTile2DRangeK * kParallelize4DTile2DRangeL);
 }
 
+static void ComputeNothing5D(void*, size_t, size_t, size_t, size_t, size_t) {
+}
+
+TEST(Parallelize5D, SingleThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_5d(threadpool.get(),
+		ComputeNothing5D,
+		nullptr,
+		kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+		0 /* flags */);
+}
+
+TEST(Parallelize5D, MultiThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d(
+		threadpool.get(),
+		ComputeNothing5D,
+		nullptr,
+		kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+		0 /* flags */);
+}
+
+static void CheckBounds5D(void*, size_t i, size_t j, size_t k, size_t l, size_t m) {
+	EXPECT_LT(i, kParallelize5DRangeI);
+	EXPECT_LT(j, kParallelize5DRangeJ);
+	EXPECT_LT(k, kParallelize5DRangeK);
+	EXPECT_LT(l, kParallelize5DRangeL);
+	EXPECT_LT(m, kParallelize5DRangeM);
+}
+
+TEST(Parallelize5D, SingleThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_5d(
+		threadpool.get(),
+		CheckBounds5D,
+		nullptr,
+		kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+		0 /* flags */);
+}
+
+TEST(Parallelize5D, MultiThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d(
+		threadpool.get(),
+		CheckBounds5D,
+		nullptr,
+		kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+		0 /* flags */);
+}
+
+static void SetTrue5D(std::atomic_bool* processed_indicators, size_t i, size_t j, size_t k, size_t l, size_t m) {
+	const size_t linear_idx = (((i * kParallelize5DRangeJ + j) * kParallelize5DRangeK + k) * kParallelize5DRangeL + l) * kParallelize5DRangeM + m;
+	processed_indicators[linear_idx].store(true, std::memory_order_relaxed);
+}
+
+TEST(Parallelize5D, SingleThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize5DRangeI * kParallelize5DRangeJ * kParallelize5DRangeK * kParallelize5DRangeL * kParallelize5DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_5d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_t>(SetTrue5D),
+		static_cast<void*>(indicators.data()),
+		kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize5DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DRangeJ + j) * kParallelize5DRangeK + k) * kParallelize5DRangeL + l) * kParallelize5DRangeM + m;
+						EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") not processed";
+					}
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize5D, MultiThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize5DRangeI * kParallelize5DRangeJ * kParallelize5DRangeK * kParallelize5DRangeL * kParallelize5DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_t>(SetTrue5D),
+		static_cast<void*>(indicators.data()),
+		kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize5DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DRangeJ + j) * kParallelize5DRangeK + k) * kParallelize5DRangeL + l) * kParallelize5DRangeM + m;
+						EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") not processed";
+					}
+				}
+			}
+		}
+	}
+}
+
+static void Increment5D(std::atomic_int* processed_counters, size_t i, size_t j, size_t k, size_t l, size_t m) {
+	const size_t linear_idx = (((i * kParallelize5DRangeJ + j) * kParallelize5DRangeK + k) * kParallelize5DRangeL + l) * kParallelize5DRangeM + m;
+	processed_counters[linear_idx].fetch_add(1, std::memory_order_relaxed);
+}
+
+TEST(Parallelize5D, SingleThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize5DRangeI * kParallelize5DRangeJ * kParallelize5DRangeK * kParallelize5DRangeL * kParallelize5DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_5d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_t>(Increment5D),
+		static_cast<void*>(counters.data()),
+		kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize5DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DRangeJ + j) * kParallelize5DRangeK + k) * kParallelize5DRangeL + l) * kParallelize5DRangeM + m;
+						EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") was processed "
+							<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+					}
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize5D, MultiThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize5DRangeI * kParallelize5DRangeJ * kParallelize5DRangeK * kParallelize5DRangeL * kParallelize5DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_t>(Increment5D),
+		static_cast<void*>(counters.data()),
+		kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize5DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DRangeJ + j) * kParallelize5DRangeK + k) * kParallelize5DRangeL + l) * kParallelize5DRangeM + m;
+						EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") was processed "
+							<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+					}
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize5D, SingleThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize5DRangeI * kParallelize5DRangeJ * kParallelize5DRangeK * kParallelize5DRangeL * kParallelize5DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	for (size_t iteration = 0; iteration < kIncrementIterations5D; iteration++) {
+		pthreadpool_parallelize_5d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_5d_t>(Increment5D),
+			static_cast<void*>(counters.data()),
+			kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+				0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize5DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DRangeJ + j) * kParallelize5DRangeK + k) * kParallelize5DRangeL + l) * kParallelize5DRangeM + m;
+						EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations5D)
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") was processed "
+							<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+							<< "(expected: " << kIncrementIterations5D << ")";
+					}
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize5D, MultiThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize5DRangeI * kParallelize5DRangeJ * kParallelize5DRangeK * kParallelize5DRangeL * kParallelize5DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	for (size_t iteration = 0; iteration < kIncrementIterations5D; iteration++) {
+		pthreadpool_parallelize_5d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_5d_t>(Increment5D),
+			static_cast<void*>(counters.data()),
+			kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+				0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize5DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DRangeJ + j) * kParallelize5DRangeK + k) * kParallelize5DRangeL + l) * kParallelize5DRangeM + m;
+						EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations5D)
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") was processed "
+							<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+							<< "(expected: " << kIncrementIterations5D << ")";
+					}
+				}
+			}
+		}
+	}
+}
+
+static void IncrementSame5D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t k, size_t l, size_t m) {
+	num_processed_items->fetch_add(1, std::memory_order_relaxed);
+}
+
+TEST(Parallelize5D, MultiThreadPoolHighContention) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_t>(IncrementSame5D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize5DRangeI * kParallelize5DRangeJ * kParallelize5DRangeK * kParallelize5DRangeL * kParallelize5DRangeM);
+}
+
+static void WorkImbalance5D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t k, size_t l, size_t m) {
+	num_processed_items->fetch_add(1, std::memory_order_relaxed);
+	if (i == 0 && j == 0 && k == 0 && l == 0 && m == 0) {
+		/* Spin-wait until all items are computed */
+		while (num_processed_items->load(std::memory_order_relaxed) != kParallelize5DRangeI * kParallelize5DRangeJ * kParallelize5DRangeK * kParallelize5DRangeL * kParallelize5DRangeM) {
+			std::atomic_thread_fence(std::memory_order_acquire);
+		}
+	}
+}
+
+TEST(Parallelize5D, MultiThreadPoolWorkStealing) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_t>(WorkImbalance5D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize5DRangeI, kParallelize5DRangeJ, kParallelize5DRangeK, kParallelize5DRangeL, kParallelize5DRangeM,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize5DRangeI * kParallelize5DRangeJ * kParallelize5DRangeK * kParallelize5DRangeL * kParallelize5DRangeM);
+}
+
+static void ComputeNothing5DTile1D(void*, size_t, size_t, size_t, size_t, size_t, size_t) {
+}
+
+TEST(Parallelize5DTile1D, SingleThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_5d_tile_1d(threadpool.get(),
+		ComputeNothing5DTile1D,
+		nullptr,
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+}
+
+TEST(Parallelize5DTile1D, MultiThreadPoolCompletes) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		ComputeNothing5DTile1D,
+		nullptr,
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+}
+
+static void CheckBounds5DTile1D(void*, size_t i, size_t j, size_t k, size_t l, size_t start_m, size_t tile_m) {
+	EXPECT_LT(i, kParallelize5DTile1DRangeI);
+	EXPECT_LT(j, kParallelize5DTile1DRangeJ);
+	EXPECT_LT(k, kParallelize5DTile1DRangeK);
+	EXPECT_LT(l, kParallelize5DTile1DRangeL);
+	EXPECT_LT(start_m, kParallelize5DTile1DRangeM);
+	EXPECT_LE(start_m + tile_m, kParallelize5DTile1DRangeM);
+}
+
+TEST(Parallelize5DTile1D, SingleThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		CheckBounds5DTile1D,
+		nullptr,
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+}
+
+TEST(Parallelize5DTile1D, MultiThreadPoolAllItemsInBounds) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		CheckBounds5DTile1D,
+		nullptr,
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+}
+
+static void CheckTiling5DTile1D(void*, size_t i, size_t j, size_t k, size_t l, size_t start_m, size_t tile_m) {
+	EXPECT_GT(tile_m, 0);
+	EXPECT_LE(tile_m, kParallelize5DTile1DTileM);
+	EXPECT_EQ(start_m % kParallelize5DTile1DTileM, 0);
+	EXPECT_EQ(tile_m, std::min<size_t>(kParallelize5DTile1DTileM, kParallelize5DTile1DRangeM - start_m));
+}
+
+TEST(Parallelize5DTile1D, SingleThreadPoolUniformTiling) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		CheckTiling5DTile1D,
+		nullptr,
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+}
+
+TEST(Parallelize5DTile1D, MultiThreadPoolUniformTiling) {
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		CheckTiling5DTile1D,
+		nullptr,
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+}
+
+static void SetTrue5DTile1D(std::atomic_bool* processed_indicators, size_t i, size_t j, size_t k, size_t l, size_t start_m, size_t tile_m) {
+	for (size_t m = start_m; m < start_m + tile_m; m++) {
+		const size_t linear_idx = (((i * kParallelize5DTile1DRangeJ + j) * kParallelize5DTile1DRangeK + k) * kParallelize5DTile1DRangeL + l) * kParallelize5DTile1DRangeM + m;
+		processed_indicators[linear_idx].store(true, std::memory_order_relaxed);
+	}
+}
+
+TEST(Parallelize5DTile1D, SingleThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize5DTile1DRangeI * kParallelize5DTile1DRangeJ * kParallelize5DTile1DRangeK * kParallelize5DTile1DRangeL * kParallelize5DTile1DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_tile_1d_t>(SetTrue5DTile1D),
+		static_cast<void*>(indicators.data()),
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize5DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DTile1DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DTile1DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DTile1DRangeJ + j) * kParallelize5DTile1DRangeK + k) * kParallelize5DTile1DRangeL + l) * kParallelize5DTile1DRangeM + m;
+						EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") not processed";
+					}
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize5DTile1D, MultiThreadPoolAllItemsProcessed) {
+	std::vector<std::atomic_bool> indicators(kParallelize5DTile1DRangeI * kParallelize5DTile1DRangeJ * kParallelize5DTile1DRangeK * kParallelize5DTile1DRangeL * kParallelize5DTile1DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_tile_1d_t>(SetTrue5DTile1D),
+		static_cast<void*>(indicators.data()),
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize5DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DTile1DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DTile1DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DTile1DRangeJ + j) * kParallelize5DTile1DRangeK + k) * kParallelize5DTile1DRangeL + l) * kParallelize5DTile1DRangeM + m;
+						EXPECT_TRUE(indicators[linear_idx].load(std::memory_order_relaxed))
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") not processed";
+					}
+				}
+			}
+		}
+	}
+}
+
+static void Increment5DTile1D(std::atomic_int* processed_counters, size_t i, size_t j, size_t k, size_t l, size_t start_m, size_t tile_m) {
+	for (size_t m = start_m; m < start_m + tile_m; m++) {
+		const size_t linear_idx = (((i * kParallelize5DTile1DRangeJ + j) * kParallelize5DTile1DRangeK + k) * kParallelize5DTile1DRangeL + l) * kParallelize5DTile1DRangeM + m;
+		processed_counters[linear_idx].fetch_add(1, std::memory_order_relaxed);
+	}
+}
+
+TEST(Parallelize5DTile1D, SingleThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize5DTile1DRangeI * kParallelize5DTile1DRangeJ * kParallelize5DTile1DRangeK * kParallelize5DTile1DRangeL * kParallelize5DTile1DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_tile_1d_t>(Increment5DTile1D),
+		static_cast<void*>(counters.data()),
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize5DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DTile1DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DTile1DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DTile1DRangeJ + j) * kParallelize5DTile1DRangeK + k) * kParallelize5DTile1DRangeL + l) * kParallelize5DTile1DRangeM + m;
+						EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") was processed "
+							<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+					}
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize5DTile1D, MultiThreadPoolEachItemProcessedOnce) {
+	std::vector<std::atomic_int> counters(kParallelize5DTile1DRangeI * kParallelize5DTile1DRangeJ * kParallelize5DTile1DRangeK * kParallelize5DTile1DRangeL * kParallelize5DTile1DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_tile_1d_t>(Increment5DTile1D),
+		static_cast<void*>(counters.data()),
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+
+	for (size_t i = 0; i < kParallelize5DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DTile1DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DTile1DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DTile1DRangeJ + j) * kParallelize5DTile1DRangeK + k) * kParallelize5DTile1DRangeL + l) * kParallelize5DTile1DRangeM + m;
+						EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), 1)
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") was processed "
+							<< counters[linear_idx].load(std::memory_order_relaxed) << " times (expected: 1)";
+					}
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize5DTile1D, SingleThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize5DTile1DRangeI * kParallelize5DTile1DRangeJ * kParallelize5DTile1DRangeK * kParallelize5DTile1DRangeL * kParallelize5DTile1DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(1), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	for (size_t iteration = 0; iteration < kIncrementIterations5D; iteration++) {
+		pthreadpool_parallelize_5d_tile_1d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_5d_tile_1d_t>(Increment5DTile1D),
+			static_cast<void*>(counters.data()),
+			kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+			kParallelize5DTile1DTileM,
+			0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize5DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DTile1DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DTile1DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DTile1DRangeJ + j) * kParallelize5DTile1DRangeK + k) * kParallelize5DTile1DRangeL + l) * kParallelize5DTile1DRangeM + m;
+						EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations5D)
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") was processed "
+							<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+							<< "(expected: " << kIncrementIterations5D << ")";
+					}
+				}
+			}
+		}
+	}
+}
+
+TEST(Parallelize5DTile1D, MultiThreadPoolEachItemProcessedMultipleTimes) {
+	std::vector<std::atomic_int> counters(kParallelize5DTile1DRangeI * kParallelize5DTile1DRangeJ * kParallelize5DTile1DRangeK * kParallelize5DTile1DRangeL * kParallelize5DTile1DRangeM);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	for (size_t iteration = 0; iteration < kIncrementIterations5D; iteration++) {
+		pthreadpool_parallelize_5d_tile_1d(
+			threadpool.get(),
+			reinterpret_cast<pthreadpool_task_5d_tile_1d_t>(Increment5DTile1D),
+			static_cast<void*>(counters.data()),
+			kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+			kParallelize5DTile1DTileM,
+			0 /* flags */);
+	}
+
+	for (size_t i = 0; i < kParallelize5DTile1DRangeI; i++) {
+		for (size_t j = 0; j < kParallelize5DTile1DRangeJ; j++) {
+			for (size_t k = 0; k < kParallelize5DTile1DRangeK; k++) {
+				for (size_t l = 0; l < kParallelize5DTile1DRangeL; l++) {
+					for (size_t m = 0; m < kParallelize5DTile1DRangeM; m++) {
+						const size_t linear_idx = (((i * kParallelize5DTile1DRangeJ + j) * kParallelize5DTile1DRangeK + k) * kParallelize5DTile1DRangeL + l) * kParallelize5DTile1DRangeM + m;
+						EXPECT_EQ(counters[linear_idx].load(std::memory_order_relaxed), kIncrementIterations5D)
+							<< "Element (" << i << ", " << j << ", " << k << ", " << l << ", " << m << ") was processed "
+							<< counters[linear_idx].load(std::memory_order_relaxed) << " times "
+							<< "(expected: " << kIncrementIterations5D << ")";
+					}
+				}
+			}
+		}
+	}
+}
+
+static void IncrementSame5DTile1D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t k, size_t l, size_t start_m, size_t tile_m) {
+	for (size_t m = start_m; m < start_m + tile_m; m++) {
+		num_processed_items->fetch_add(1, std::memory_order_relaxed);
+	}
+}
+
+TEST(Parallelize5DTile1D, MultiThreadPoolHighContention) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_tile_1d_t>(IncrementSame5DTile1D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize5DTile1DRangeI * kParallelize5DTile1DRangeJ * kParallelize5DTile1DRangeK * kParallelize5DTile1DRangeL * kParallelize5DTile1DRangeM);
+}
+
+static void WorkImbalance5DTile1D(std::atomic_int* num_processed_items, size_t i, size_t j, size_t k, size_t l, size_t start_m, size_t tile_m) {
+	num_processed_items->fetch_add(tile_m, std::memory_order_relaxed);
+	if (i == 0 && j == 0 && k == 0 && l == 0 && start_m == 0) {
+		/* Spin-wait until all items are computed */
+		while (num_processed_items->load(std::memory_order_relaxed) != kParallelize5DTile1DRangeI * kParallelize5DTile1DRangeJ * kParallelize5DTile1DRangeK * kParallelize5DTile1DRangeL * kParallelize5DTile1DRangeM) {
+			std::atomic_thread_fence(std::memory_order_acquire);
+		}
+	}
+}
+
+TEST(Parallelize5DTile1D, MultiThreadPoolWorkStealing) {
+	std::atomic_int num_processed_items = ATOMIC_VAR_INIT(0);
+
+	auto_pthreadpool_t threadpool(pthreadpool_create(0), pthreadpool_destroy);
+	ASSERT_TRUE(threadpool.get());
+
+	if (pthreadpool_get_threads_count(threadpool.get()) <= 1) {
+		GTEST_SKIP();
+	}
+
+	pthreadpool_parallelize_5d_tile_1d(
+		threadpool.get(),
+		reinterpret_cast<pthreadpool_task_5d_tile_1d_t>(WorkImbalance5DTile1D),
+		static_cast<void*>(&num_processed_items),
+		kParallelize5DTile1DRangeI, kParallelize5DTile1DRangeJ, kParallelize5DTile1DRangeK, kParallelize5DTile1DRangeL, kParallelize5DTile1DRangeM,
+		kParallelize5DTile1DTileM,
+		0 /* flags */);
+	EXPECT_EQ(num_processed_items.load(std::memory_order_relaxed), kParallelize5DTile1DRangeI * kParallelize5DTile1DRangeJ * kParallelize5DTile1DRangeK * kParallelize5DTile1DRangeL * kParallelize5DTile1DRangeM);
+}
+
 static void ComputeNothing5DTile2D(void*, size_t, size_t, size_t, size_t, size_t, size_t, size_t) {
 }