Missing barrier on stencil code after parallelizing

[jdsteve2]

This 5-point stencil kernel produces the following OpenCL code without a barrier, which would be required for correctness.

import loopy as lp
import numpy as np
knl = lp.make_kernel(
    "[nx,nt] -> {[x, t]: 0<=x<nx and 0<=t<nt}",
    """
    u[t+2,x+1] = 2*u[t+1,x+1] + dt**2/dx**2 * (u[t+1,x+2] - 2*u[t+1,x+1] + u[t+1,x]) - u[t,x+1]
    """
    ) 
knl = lp.add_dtypes(knl, {"u": np.float32, "dx": np.float32, "dt": np.float32})
knl = lp.split_iname(knl, "x", 14)
knl = lp.assume(knl, "nx % 14 = 0 and nt >= 1 and nx >= 1")
knl = lp.tag_inames(knl, "x_outer:g.0, x_inner:l.0")
print(lp.generate_code_v2(knl).device_code())

Generated code:

#define lid(N) ((int) get_local_id(N))
#define gid(N) ((int) get_group_id(N))
__kernel void __attribute__ ((reqd_work_group_size(14, 1, 1))) loopy_kernel(float const dt, float const dx, int const nt, int const nx, __global float *__restrict__ u)
{
  for (int t = 0; t <= -1 + nt; ++t)
    u[(2 + nx) * (2 + t) + 1 + 14 * gid(0) + lid(0)] = 2.0f * u[(2 + nx) * (1 + t) + 1 + 14 * gid(0) + lid(0)] + ((dt * dt) / (dx * dx)) * (u[(2 + nx) * (1 + t) + 2 + 14 * gid(0) + lid(0)] + -1.0f * 2.0f * u[(2 + nx) * (1 + t) + 1 + 14 * gid(0) + lid(0)] + u[(2 + nx) * (1 + t) + 14 * gid(0) + lid(0)]) + -1.0f * u[(2 + nx) * t + 1 + 14 * gid(0) + lid(0)];
}

[kaushikcfd]

Currently the dependency checking in the barrier placement code does not take into account that an instruction in a loop might induce WW, RW, WR races on its own. That's an easy fix here:

loopy/loopy/schedule/__init__.py

Lines 1674 to 1679 in f7a5a0f

	target_deps = self.kernel.recursive_insn_dep_map()[target.id]
	if source.id in target_deps:
	if self.reverse:
	dep_descr = "{tgt} rev-depends on {src}"
	else:
	dep_descr = "{tgt} depends on {src}"

But the catch is that currently loopy.symbolic.AccessRangeOverlapChecker is quite imprecise and would lead to false positives even for cases like out[i] = 2*in[i] when i goes along the device's work item loop.

[inducer]

As @kaushikcfd points out, a big part of the issue is the absence of self-dependencies in the current, simplistic dependency model. If there are two statements like so:

import loopy as lp
import numpy as np
knl = lp.make_kernel(
    "[nx,nt] -> {[x, t]: 0<=x<nx and 0<=t<nt}",
    """
    u[t+2,x+1] = 2*u[t+1,x+1] + dt**2/dx**2 * (u[t+1,x+2] - 2*u[t+1,x+1] + u[t+1,x]) - u[t,x+1] {id=a}
    u[t+2,x+1] = 2*u[t+1,x+1] + dt**2/dx**2 * (u[t+1,x+2] - 2*u[t+1,x+1] + u[t+1,x]) - u[t,x+1] {dep=a}
    """
    ) 
knl = lp.add_dtypes(knl, {"u": np.float32, "dx": np.float32, "dt": np.float32})
knl = lp.split_iname(knl, "x", 14)
knl = lp.assume(knl, "nx % 14 = 0 and nt >= 1 and nx >= 1")
knl = lp.tag_inames(knl, "x_outer:g.0, x_inner:l.0")
print(lp.generate_code_v2(knl).device_code())

Loopy will notice that something is off:

loopy.diagnostic.MissingBarrierError: loopy_kernel: Dependency 'insn depends on a' (for variable 'u') requires synchronizatio
n by a global barrier (add a 'no_sync_with' instruction option to state that no synchronization is needed)

@jdsteve2's work towards better dependency tracking (#168) should help address this.

[kaushikcfd]

Applying

diff --git a/loopy/schedule/__init__.py b/loopy/schedule/__init__.py
index 12fdf6b0..496502a8 100644
--- a/loopy/schedule/__init__.py
+++ b/loopy/schedule/__init__.py
@@ -1672,6 +1672,10 @@ class DependencyTracker:
             source, target = target, source
 
         target_deps = self.kernel.recursive_insn_dep_map()[target.id]
+
+        if self.reverse:
+            target_deps = target_deps | {target.id}
+
         if source.id in target_deps:
             if self.reverse:
                 dep_descr = "{tgt} rev-depends on {src}"

to #224 and transforming the kernel as : knl = lp.set_options(knl, "insert_additional_gbarriers"), we obtain the following linearization:

LINEARIZATION:
   0: for t
   1:     ... gbarrier
   2:     CALL KERNEL loopy_kernel(extra_args=[], extra_inames=['t'])
   3:         u[t + 2, 1 + x_inner + x_outer*14] = 2*u[t + 1, 1 + x_inner + x_outer*14] + (dt**2 / dx**2)*(u[t + 1, 2 + x_inner + x_outer*14] + (-1)*2*u[t + 1, 1 + x_inner + x_outer*14] + u[t + 1, x_inner + x_outer*14]) + (-1)*u[t, 1 + x_inner + x_outer*14]  {id=insn}
   4:     RETURN FROM KERNEL loopy_kernel
   5: end t