Skip to content

Commit c6a97c4

Browse files
authored
Merge pull request #3080 from devitocodes/subgrid-strided-subsampling
API: Subgrid (integer) support
2 parents 6d2ad82 + ad24b75 commit c6a97c4

8 files changed

Lines changed: 3191 additions & 10 deletions

File tree

‎devito/ir/support/basic.py‎

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -869,8 +869,15 @@ def is_storage_related(self, dims=None):
869869
"""
870870
for d in self.findices:
871871
if d._defines & set(as_tuple(dims)): # noqa: SIM102
872-
if any(i.is_NonlinearDerived for i in d._defines) or \
873-
self.is_const(d):
872+
# An injective Dimension maps distinct iterations to distinct
873+
# locations by construction -- e.g. a factor-based
874+
# ConditionalDimension, whose guard admits only
875+
# `parent % factor == 0` and whose `i -> i//factor` is
876+
# one-to-one on those. Unlike, say, a SteppingDimension, where
877+
# `time % 2` deliberately revisits the same slot
878+
nonlinear = [i for i in d._defines
879+
if i.is_NonlinearDerived and not i.is_Injective]
880+
if nonlinear or self.is_const(d):
874881
return True
875882
return False
876883

‎devito/operator/operator.py‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -34,7 +34,7 @@
3434
from devito.passes import (
3535
Graph, error_mapper, finalize_args, generate_implicit, generate_macros, is_on_device,
3636
lower_dtypes, lower_index_derivatives, lower_local_sums, minimize_symbols,
37-
optimize_pows, unevaluate
37+
optimize_pows, subsampling, unevaluate
3838
)
3939
from devito.symbolics import estimate_cost, subs_op_args
4040
from devito.tools import (
@@ -423,6 +423,11 @@ def _lower_clusters(cls, expressions, profiler=None, **kwargs):
423423
# Generate implicit Clusters from higher level abstractions
424424
clusters = generate_implicit(clusters)
425425

426+
# Strengthen factor-based ConditionalDimension guards into strided
427+
# iteration. Before `lower_index_derivatives` so that the Dimension
428+
# swap is seen by everything downstream
429+
clusters = subsampling(clusters)
430+
426431
# Lower all remaining high order symbolic objects
427432
clusters = lower_local_sums(clusters, **kwargs)
428433
clusters = lower_index_derivatives(clusters, **kwargs)

‎devito/passes/clusters/__init__.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
from .factorization import * # noqa
66
from .fusion import * # noqa
77
from .blocking import * # noqa
8+
from .subsampling import * # noqa
89
from .asynchrony import * # noqa
910
from .implicit import * # noqa
1011
from .misc import * # noqa
Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
from devito.ir.clusters import Queue
2+
from devito.ir.support import GuardFactor
3+
from devito.symbolics import IntDiv, uxreplace
4+
from devito.tools import timed_pass
5+
6+
__all__ = ['subsampling']
7+
8+
9+
@timed_pass()
10+
def subsampling(clusters):
11+
"""
12+
Strengthen factor-based ConditionalDimension guards into strided iteration.
13+
14+
A Function subsampled by `factor` along `d` is written under a
15+
`d % factor == 0` guard, so the enclosing loop visits every point of `d`
16+
and discards all but one in `factor` of them. Where legal, drop the guard
17+
and step the loop `factor` points at a time instead
18+
19+
for (x = x_m; x <= x_M; x += 1) for (x = ...; x <= x_M; x += xf)
20+
if (x%xf == 0) -> g[x/xf] = ...
21+
g[x/xf] = ...
22+
23+
The loop is taken over the subsampled Dimension `cd` itself, with a *unit*
24+
stride, and every access to the parent `d` is rescaled to `cd*factor`:
25+
26+
for (x = x_m; x <= x_M; x += 1) for (xs = xs_m; xs <= xs_M; xs++)
27+
if (x%xf == 0) -> g[xs] = u[xs*xf]
28+
g[x/xf] = u[x]
29+
30+
A non-unitary increment is deliberately avoided: it would escape the
31+
Lamport-style dependence analysis the compiler relies on. The loop bounds
32+
are `cd`'s own runtime values, which `ConditionalDimension._arg_values`
33+
derives as `ceil(x_m/factor) .. floor(x_M/factor)`, so no anchoring of the
34+
loop minimum is needed either.
35+
36+
Parameters
37+
----------
38+
clusters : tuple of Clusters
39+
Input Clusters, subject of the optimization pass.
40+
"""
41+
return Subsampling().process(clusters)
42+
43+
44+
class Subsampling(Queue):
45+
46+
"""
47+
Replace uniform factor-based guards with strided iteration.
48+
49+
Clusters are grouped by iteration space prefix alone, deliberately *not*
50+
by guard: a guard may only become a stride if every Cluster sharing the
51+
prefix carries the very same one. Striding a loop nest that also hosts an
52+
unguarded expression would run that expression on a subset of its
53+
iterations -- and dropping the guard is what makes the nest shared in the
54+
first place, since `stree_build` reuses a subtree across Clusters with
55+
matching guards (see `reuse_whole_subtree`).
56+
"""
57+
58+
def callback(self, clusters, prefix):
59+
if not prefix:
60+
return clusters
61+
62+
d = prefix[-1].dim
63+
64+
guards = {c.guards.get(d) for c in clusters}
65+
if len(guards) > 1:
66+
# Non-uniformly guarded, hence `d`'s loop nest is shared with
67+
# expressions that must keep running at every point
68+
return clusters
69+
70+
guard = guards.pop()
71+
if not all(_is_strideable(c, d, guard) for c in clusters):
72+
return clusters
73+
74+
cd = guard.d
75+
factor = cd.symbolic_factor
76+
77+
processed = []
78+
for c in clusters:
79+
# Iterate the subsampled Dimension itself, with unit stride, and
80+
# rescale every parent access `d` to `cd*factor`
81+
ispace = c.ispace.switch(d, cd)
82+
83+
# `IntDiv(d, factor)` is how the subsampled access was lowered;
84+
# with `d -> cd*factor` it collapses to `cd`, so rewrite it
85+
# directly rather than leaving `(cd*factor)/factor` behind
86+
subs = {IntDiv(d, factor): cd, d: cd*factor}
87+
exprs = [uxreplace(e, subs) for e in c.exprs]
88+
89+
# `cd` inherits `d`'s computational properties: it covers the very
90+
# same points, just addressed in subsampled coordinates. Without
91+
# this the loop would lose PARALLEL and so all of its OpenMP
92+
properties = c.properties.drop(d)
93+
properties = properties.add(cd, c.properties[d])
94+
95+
processed.append(c.rebuild(exprs=exprs, ispace=ispace,
96+
guards=c.guards.popany(d),
97+
properties=properties))
98+
99+
return processed
100+
101+
102+
def _is_strideable(cluster, d, guard):
103+
"""
104+
True if iteration over `d` may step `guard`'s subsampling factor points at
105+
a time in place of evaluating `guard` at every point.
106+
"""
107+
# The guard must be exactly a subsampling guard. Anything `And`/`Or`-ed in
108+
# (e.g. a user-supplied `condition`) must still be evaluated pointwise
109+
if not isinstance(guard, GuardFactor):
110+
return False
111+
112+
if d not in cluster.ispace.itdims:
113+
return False
114+
115+
# Skipping points is only sound if no dependence is carried along `d`;
116+
# otherwise the discarded iterations may be doing work that `d`'s
117+
# remaining iterations rely upon
118+
if not cluster.properties.is_parallel_relaxed(d):
119+
return False
120+
121+
# Non-affine accesses may read the points that would be skipped
122+
return cluster.properties.is_affine(d)

‎devito/types/dimension.py‎

Lines changed: 22 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -116,6 +116,7 @@ class Dimension(ArgProvider):
116116
is_Custom = False
117117
is_Derived = False
118118
is_NonlinearDerived = False
119+
is_Injective = False
119120
is_AbstractSub = False
120121
is_Sub = False
121122
is_MultiSub = False
@@ -1090,6 +1091,19 @@ def __init_finalize__(self, name, parent=None, factor=None, condition=None,
10901091
self._indirect = indirect
10911092
self._relation = relation
10921093

1094+
@property
1095+
def is_Injective(self):
1096+
"""
1097+
True if distinct iterations of the parent that pass the guard map to
1098+
distinct indices of `self`, so that no two iterations can access the
1099+
same memory location through it.
1100+
1101+
A pure factor subsampling is injective: the guard admits only
1102+
`parent % factor == 0`, and `i -> i//factor` is one-to-one on those.
1103+
An arbitrary `condition` carries no such promise.
1104+
"""
1105+
return self._factor is not None and self._condition is None
1106+
10931107
@property
10941108
def uses_symbolic_factor(self):
10951109
return self._factor is not None
@@ -1153,13 +1167,18 @@ def _arg_values(self, interval, grid=None, args=None, **kwargs):
11531167
if isinstance(fact, Constant):
11541168
fact = fact.data
11551169

1156-
toint = lambda x: math.ceil(x / fact)
11571170
vals = {}
1171+
1172+
# The first subsampled point at or after the parent's min...
11581173
with suppress(KeyError, TypeError):
1159-
vals[self.min_name] = toint(kwargs.get(self.parent.min_name))
1174+
vals[self.min_name] = math.ceil(kwargs.get(self.parent.min_name) / fact)
11601175

1176+
# ...and the last one at or before its max. Note `floor` rather than
1177+
# `ceil`: with e.g. `factor=3` and `x_M=11` the last subsampled point
1178+
# is `x=9`, hence index 3, while `ceil(11/3) == 4` would overshoot the
1179+
# Function by one
11611180
with suppress(KeyError, TypeError):
1162-
vals[self.max_name] = toint(kwargs.get(self.parent.max_name))
1181+
vals[self.max_name] = kwargs.get(self.parent.max_name) // fact
11631182

11641183
vals[self.symbolic_factor.name] = fact
11651184

0 commit comments

Comments
 (0)