Hi! Before OpenMP 5.0, all OpenMP loop nests had to be rectangular (and OpenMP has various other restrictions that make implementation easier), so that it was very easy to compute the number of iterations of the collapsed loop by computing number of iterations of each of the loops in the loop nest and multiplying them; then either inline code or runtime library from that total number of iterations, the thread number, number of threads in the team and scheduling policies can determine an interval of logical iterations which the thread will handle and then it was fairly cheap to compute from that single 0 .. num_iterations-1 iterator what iterator values it starts with.
For non-rectangular loops, this is getting harder, both how to compute the number of iterations and how to cheaply if possible from the single logical number iteration start compute the different iterator values. Below is a proof of concept of what I came up so far. I have tried to list the relevant restrictions in OpenMP in the first big comment in bar function, foo function has #if 0 what I'm trying to implement and #else a version that does it serially in a single thread only, and then bar is an attempt to write in C what GCC could roughly emit for it. Note, the test uses (mostly) int types, but in reality it can be other integral types, signed or unsigned, and OpenMP just adds assurances that the loops will not wrap around and otherwise behave well. For both number of iterations computation and the logical iteration to actual iterator values computations I have two versions, one lame fallback which worst case let's each thread basically run the whole loop as is except for the body. This one is because I don't want to spend months on it and deal with Bernoulli constants for 120 nested loops etc., so something that will be standard compliant and for loops that have costly body will be beneficial too. And then an optimized version, for now limited to triangular loops (can be used also if there are rectangular loops around those), where the first part tries to compute total number of iterations using Faulhaber's formula and the second part attempts to compute quadratic equation root using integer square root. The proof of concept right now uses the fallback even if it sees that the inner loop will not have at least a single iteration for all values of the outer loop iterator, will try to change that to just artificially change the a and b values (bounds of the outer iterator) for the purposes of total number of iterations computation and for the purpose of transforming the quadratic equation root into the actual iterator values (and keep the original ones for the purpose of computing lastprivate iterator values). Any thoughts on how to simplify this, what to do differently, whether e.g. using floating point math instead would be beneficial etc.? Any input appreciated. /* Proof of concept for OpenMP non-rectangular worksharing-loop implementation. Copyright (C) 2020 Free Software Foundation, Inc. GCC is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3, or (at your option) any later version. GCC is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with GCC; see the file COPYING3. If not see <http://www.gnu.org/licenses/>. */ #include <stdlib.h> #include <stdio.h> #ifdef _OPENMP #include <omp.h> #endif int x, i, j, nitersv, niterscnt; #ifdef DEBUG #define DPRINTF(...) printf (__VA_ARGS__) #else #define DPRINTF(...) do {} while (0) #endif #ifndef _OPENMP #define omp_get_num_threads() 1 #define omp_get_thread_num() 0 #endif void foo (int a, int b, int c, int d, int e, int f, int g, int h) { #ifdef DEBUG #pragma omp single nitersv = 0; #endif #if 0 #pragma omp for collapse(2) lastprivate (i, j, x) #else #pragma omp single #endif for (i = a; i < b; i += c) for (j = d * i + e; j < f * i + g; j += h) { x = i * 1024 + (j & 1023); #ifdef DEBUG nitersv++; #endif DPRINTF ("%d %d %d %d\n", i, j, x, omp_get_thread_num ()); } #ifdef DEBUG #pragma omp single DPRINTF ("niters = %d\n", nitersv); #endif } void bar (int a, int b, int c, int d, int e, int f, int g, int h) { /* Proposed implementation of: #pragma omp for collapse(2) lastprivate (i, j, x) for (i = a; i < b; i += c) for (j = d * i + e; j < f * i + g; j += h) */ /* OpenMP requires that ((f - d) * c) % h == 0 and that either the initializer and condition expressions are outermost loop invariant, or have syntactic forms that can be represented as integral a1 * var-outer + a2 where var-outer is some outer loop iterator with compatible type and a1 and a2 are integral expressions (and have compatible type too), which are outermost loop invariant. Comparisons can be <, <=, >, >= or != but in the last case the step is required to be compile time constant so that one can determine iteration direction and for the others the step has to match the iteration direction of the comparison operator. And also a requirement that essentially says that there is no wrap-around in any of the iterators and that the number of iterations can be computed without risks of overflows/wrap-arounds. Any number of loops can be collapsed and all but the outermost can be non-rectangular (or at least potentially one, basically where the expressions refer to the outer loop iterator). The step expressions must be always outermost loop invariant. */ /* First try to calculate the total number of iterations. Can be simplified by computing all outermost rectangular loops whose iterator vars are not referenced in the non-rectangular loops separately, and similarly all innermost rectangular loops separately. */ int niters = 0; if (!(a < b)) { /* No iterations at all, only i defined after loop. */ i = a; goto end; } /* If the (middle) non-rectangular loops are triangular (or perhaps in some more cases using Faulhaber's formula?), check if for all the iterators the inner loop will have at least one iteration. If all of a, b, c, d, e, f, g, h are compile time constants, we want to compute niters at compile time obviously. If only some of them are constant, let the normal optimizations simplify the expressions correspondingly. */ int t4 = (b + (c - 1) - a) / c; int t5 = a + ((t4 - 1) * c); int t8 = d * a + e; int t9 = f * a + g; int t10, t11; if (t8 < t9 && d * t5 + e < f * t5 + g) { t10 = ((t9 + (h - 1) - t8) / h); t11 = ((f - d) * c / h); niters = t4 * t10 + t11 * (((t4 - 1) * t4) / 2); } else { t10 = t11 = 0; /* Fallback implementation, if it above gets too ugly/hard. Repeat all loops except the innermost, hope loop optimizations optimize at least something. */ for (int t1 = a; t1 < b; t1 += c) { int t2 = d * t1 + e; int t3 = f * t1 + g; if (t2 < t3) niters += (t3 + (h - 1) - t2) / h; } } /* Second step, the usual GCC OpenMP schedule(static) computation of which iterations should the current thread take. */ int num_threads = omp_get_num_threads (); int thr = omp_get_thread_num (); int t6 = niters / num_threads; int t7 = niters % num_threads; if (thr < t7) { t7 = 0; t6++; } int start = t6 * thr + t7; int end = start + t6; if (!(start < end)) goto end; /* Now, start contains the first logical iteration that this thread should handle (counted from 0) and end should be the first one it should already not handle. */ /* Third step, from start try to determine the initial values of the loop iteration variables. */ int ipriv, jpriv; if (t10) { /* We want to find maximum x such that start >= x * t10 + t11 * (((x - 1) * x) / 2) and from that x compute both indexes: ipriv = a + x * c; jpriv = d * ipriv + e + (start - (x * t10 + t11 * (((x - 1) * x) / 2))) * h. */ /* x = (isqrt((t10-t11/2)*(t10-t11/2)+start)-(t10-t11/2)) / t11 would be a rough guess that needs to be verified. */ int t12 = t10 - t11 / 2; /* Quick overflow check for the ^2. */ if (__builtin_expect (t12 + 45000U > 90000U, 0)) goto fallback; unsigned t13 = t12 * t12 + (unsigned) start; if (__builtin_expect (t13 == 0 || t11 == 0, 0)) goto fallback; /* Compute isqrt(t13). */ unsigned isqrtb = (1U << (__SIZEOF_INT__ * __CHAR_BIT__ + 1 - __builtin_clz (t13)) /2) - 1; unsigned isqrta = (isqrtb + 3) / 2; do { unsigned isqrtm = (isqrta + isqrtb) >> 1; if (isqrtm * isqrtm > t13) isqrtb = isqrtm - 1; else isqrta = isqrtm + 1; } while (isqrtb >= isqrta); unsigned isqrt = isqrta - 1; unsigned t14 = (isqrt - t12) / t11; unsigned t15 = t14 * t10 + t11 * (((t14 - 1) * t14) / 2); if (__builtin_expect (start >= t15, 1)) { unsigned t16 = t15 + t10 + t14; if (__builtin_expect (start >= t16, 0)) goto fallback; } else { unsigned t16 = t15 - t10 - t14; if (__builtin_expect (start < t16, 0)) goto fallback; t14--; t15 = t16; } ipriv = a + (int) t14 * c; jpriv = d * ipriv + e + (start - (int) t15) * h; } else { /* Fallback implementation, if it above gets too ugly/hard. Repeat all loops except the innermost, hope loop optimizations optimize at least something. */ fallback:; int cnt = 0; for (int t1 = a; t1 < b; t1 += c) { int t2 = d * t1 + e; int t3 = f * t1 + g; if (t2 < t3) { int t8 = (t3 + (h - 1) - t2) / h; if (cnt + t8 > start) { ipriv = t1; jpriv = t2 + (start - cnt) * h; goto done; } else cnt += t8; } } done:; } int jmax = f * ipriv + g; int xpriv; do { /* Now the body, with the privatized copies of the loop iterators as well as other privatized variables as usual in OpenMP. */ { xpriv = ipriv * 1024 + (jpriv & 1023); DPRINTF ("%d %d %d %d\n", ipriv, jpriv, xpriv, omp_get_thread_num ()); } #ifdef DEBUG if (ipriv < a || ipriv >= b) abort (); if ((ipriv - a) % c) abort (); if (jpriv < d * ipriv + e || jpriv >= f * ipriv + g) abort (); if ((jpriv - (d * ipriv + e)) % h) abort (); #pragma omp atomic ++niterscnt; #endif /* Use start as the logical iteration counter. */ start++; if (!(start < end)) break; /* Now bump the innermost iterator. */ jpriv += h; if (jpriv < jmax) continue; /* The outermost iterator doesn't need condition checking, we have done that already through the start < end check. */ ipriv += c; /* Or precompute earlier how to bump jmax and jmin less expensively? */ jpriv = d * ipriv + e; jmax = f * ipriv + g; } while (1); /* Lastprivate handling. */ if (start == niters) { /* The thread that has been assigned the last iteration will handle this. */ /* The variables other than iterators are very easy. */ x = xpriv; /* The iterators can be harder, at least in cases where the innermost loop is not (or might not) be executed at all for some of the outer loop iterator values. */ /* Try to do something smarter for the cases where the first phase proved that is not the case? */ /* As fallback, continue iterating with empty bodies the outer loops until all the conditions fail. */ jpriv += h; do { ipriv += c; if (!(ipriv < b)) break; jpriv = d * ipriv + e; } while (1); /* And assign those to the shared variables. */ i = ipriv; j = jpriv; } end:; DPRINTF ("niters %d\n", niters); } volatile int v; int main () { x = i = j = -1; #pragma omp parallel num_threads(15) foo (v + 4, v + 10, v + 1, v + 2, v - 9, v + 1, v, v + 1); DPRINTF ("last %d %d %d\n", i, j, x); if (i != 10 || j != 9 || x != 8 * 1024 + 7) abort (); DPRINTF ("===\n"); x = i = j = -1; #pragma omp parallel num_threads(15) bar (v + 4, v + 10, v + 1, v + 2, v - 9, v + 1, v, v + 1); DPRINTF ("last %d %d %d\n", i, j, x); if (i != 10 || j != 9 || x != 8 * 1024 + 7) abort (); DPRINTF ("===\n"); x = i = j = -1; #pragma omp parallel num_threads(15) foo (v + 1, v + 10, v + 2, v + 0, v + 1, v + 1, v + 1, v + 1); DPRINTF ("last %d %d %d\n", i, j, x); if (i != 11 || j != 10 || x != 9 * 1024 + 9) abort (); DPRINTF ("===\n"); x = i = j = -1; #pragma omp parallel num_threads(15) bar (v + 1, v + 10, v + 2, v + 0, v + 1, v + 1, v + 1, v + 1); DPRINTF ("last %d %d %d\n", i, j, x); if (i != 11 || j != 10 || x != 9 * 1024 + 9) abort (); DPRINTF ("===\n"); x = i = j = -1; #pragma omp parallel num_threads(15) foo (v + 4, v + 8, v + 12, v - 8, v - 9, v - 3, v + 6, v + 15); DPRINTF ("last %d %d %d\n", i, j, x); if (i != 16 || j != 4 || x != 5 * 1024 - 11) abort (); DPRINTF ("===\n"); x = i = j = -1; #pragma omp parallel num_threads(15) bar (v + 4, v + 8, v + 12, v - 8, v - 9, v - 3, v + 6, v + 15); DPRINTF ("last %d %d %d\n", i, j, x); if (i != 16 || j != 4 || x != 5 * 1024 - 11) abort (); DPRINTF ("===\n"); x = i = j = -1; #pragma omp parallel num_threads(15) foo (v - 13, v + 7, v + 12, v + 3, v + 5, v, v - 6, v + 1); DPRINTF ("last %d %d %d\n", i, j, x); if (i != 11 || j != 2 || x != -12 * 1024 - 7) abort (); DPRINTF ("===\n"); x = i = j = -1; #pragma omp parallel num_threads(15) bar (v - 13, v + 7, v + 12, v + 3, v + 5, v, v - 6, v + 1); DPRINTF ("last %d %d %d\n", i, j, x); if (i != 11 || j != 2 || x != -12 * 1024 - 7) abort (); DPRINTF ("===\n"); for (int idx = 0; idx < 16384 * 1024; idx++) { int a = (random () & 31) - 16; int b = (random () & 31) - 16; int c = (random () & 15) + 1; int d = (random () & 31) - 16; int e = (random () & 31) - 16; int f = (random () & 31) - 16; int g = (random () & 31) - 16; int h = (random () & 15) + 1; while (((f - d) * c % h) != 0) h = (random () & 15) + 1; x = i = j = -1; #pragma omp parallel num_threads(15) foo (a, b, c, d, e, f, g, h); int xs = x; int is = i; int js = j; x = i = j = -1; niterscnt = 0; #pragma omp parallel num_threads(15) bar (a, b, c, d, e, f, g, h); #ifdef DEBUG if (nitersv != niterscnt) abort (); #endif } return 0; } Jakub
