diff --git a/devito/ir/clusters/analysis.py b/devito/ir/clusters/analysis.py index 3f6b3099c7..4778f2b2b9 100644 --- a/devito/ir/clusters/analysis.py +++ b/devito/ir/clusters/analysis.py @@ -100,7 +100,7 @@ def _callback(self, clusters, d, prefix): # False alarm, the dependence is over a locally-defined symbol continue - if dep.is_reduction and not (d.is_Custom and d.is_Derived): + if dep.is_reduction: is_parallel_atomic = True continue diff --git a/devito/ir/support/basic.py b/devito/ir/support/basic.py index 4e26c50935..2f7a47c0e9 100644 --- a/devito/ir/support/basic.py +++ b/devito/ir/support/basic.py @@ -660,9 +660,9 @@ def is_const(self, dim): """ True if a constant dependence, that is no Dimensions involved, False otherwise. """ - return (self.source.aindices[dim] is None and - self.sink.aindices[dim] is None and - self.distance_mapper[dim] == 0) + return (self.source.aindices.get(dim, None) is None and + self.sink.aindices.get(dim, None) is None and + self.distance_mapper.get(dim, 0) == 0) @memoized_meth def is_carried(self, dim=None): diff --git a/devito/operations/interpolators.py b/devito/operations/interpolators.py index 116cf1699d..4550b5b143 100644 --- a/devito/operations/interpolators.py +++ b/devito/operations/interpolators.py @@ -3,9 +3,10 @@ import sympy from cached_property import cached_property +from devito.finite_differences.differentiable import Mul from devito.finite_differences.elementary import floor from devito.symbolics import retrieve_function_carriers, INT -from devito.tools import as_tuple, flatten, prod +from devito.tools import as_tuple, flatten from devito.types import (ConditionalDimension, Eq, Inc, Evaluable, Symbol, CustomDimension) from devito.types.utils import DimensionTuple @@ -328,11 +329,9 @@ class LinearInterpolator(WeightedInterpolator): """ @property def _weights(self): - # (1 - p) * (1 - rd) + rd * p - # simplified for better arithmetic - c = [1 - p + rd * (2*p - 1) - for (p, d, rd) in zip(self._point_symbols, self._gdim, self._rdim)] - return prod(c) + c = [(1 - p) * (1 - r) + p * r + for (p, d, r) in zip(self._point_symbols, self._gdim, self._rdim)] + return Mul(*c) @cached_property def _point_symbols(self): @@ -375,5 +374,5 @@ def interpolation_coeffs(self): @property def _weights(self): ddim, cdim = self.interpolation_coeffs.dimensions[1:] - return prod([self.interpolation_coeffs.subs({ddim: ri, cdim: rd-rd.symbolic_min}) + return Mul(*[self.interpolation_coeffs.subs({ddim: ri, cdim: rd-rd.symbolic_min}) for (ri, rd) in enumerate(self._rdim)]) diff --git a/devito/passes/iet/mpi.py b/devito/passes/iet/mpi.py index 34170bcc2e..1343a33b8a 100644 --- a/devito/passes/iet/mpi.py +++ b/devito/passes/iet/mpi.py @@ -47,9 +47,9 @@ def _drop_halospots(iet): # If a HaloSpot is outside any iteration it is not needed for iters, halo_spots in MapNodes(Iteration, HaloSpot, 'groupby').visit(iet).items(): - if not iters and halo_spots: - for hs in halo_spots: - for f in hs.fmapper: + for hs in halo_spots: + for f, v in hs.fmapper.items(): + if not iters and v.loc_indices: mapper[hs].add(f) # Transform the IET introducing the "reduced" HaloSpots diff --git a/devito/passes/iet/parpragma.py b/devito/passes/iet/parpragma.py index 4fdcd36a34..c0b6f2d155 100644 --- a/devito/passes/iet/parpragma.py +++ b/devito/passes/iet/parpragma.py @@ -285,6 +285,10 @@ def _select_candidates(self, candidates): if i.is_Vectorized: break + # Also, we do not want to collapse small atomic reductions + if i.is_ParallelAtomic and i.dim.is_Custom: + break + # Would there be enough work per parallel iteration? nested = candidates[n+1:] if nested: @@ -422,7 +426,7 @@ def _make_guard(self, parregion): def _make_nested_partree(self, partree): # Apply heuristic - if self.nhyperthreads <= self.nested: + if self.nhyperthreads <= self.nested or partree.root.is_ParallelAtomic: return partree # Note: there might be multiple sub-trees amenable to nested parallelism, diff --git a/tests/test_buffering.py b/tests/test_buffering.py index 74eea9764b..b7f59e61a5 100644 --- a/tests/test_buffering.py +++ b/tests/test_buffering.py @@ -272,7 +272,7 @@ def test_over_injection(): # Check generated code assert len(retrieve_iteration_tree(op1)) == \ - 7 + bool(configuration['language'] != 'C') + 8 + bool(configuration['language'] != 'C') buffers = [i for i in FindSymbols().visit(op1) if i.is_Array] assert len(buffers) == 1 diff --git a/tests/test_dimension.py b/tests/test_dimension.py index 0dfe89a647..b9fb7ea2f5 100644 --- a/tests/test_dimension.py +++ b/tests/test_dimension.py @@ -10,7 +10,6 @@ Dimension, DefaultDimension, SubDimension, switchconfig, SubDomain, Lt, Le, Gt, Ge, Ne, Buffer, sin, SpaceDimension, CustomDimension, dimensions, configuration) -from devito.arch.compiler import IntelCompiler, OneapiCompiler from devito.ir.iet import (Conditional, Expression, Iteration, FindNodes, FindSymbols, retrieve_iteration_tree) from devito.symbolics import indexify, retrieve_functions, IntDiv, INT @@ -1417,8 +1416,7 @@ def test_affiness(self): iterations = [i for i in FindNodes(Iteration).visit(op) if i.dim is not time] assert all(i.is_Affine for i in iterations) - @switchconfig(condition=isinstance(configuration['compiler'], - (IntelCompiler, OneapiCompiler)), safe_math=True) + @switchconfig(safe_math=True) def test_sparse_time_function(self): nt = 20 @@ -1452,6 +1450,7 @@ def test_sparse_time_function(self): # Note the endpoint of the range is 12 because we inject at p.forward for i in range(1, 12): assert p.data[i].sum() == i - 1 + print(p.data[i, 10, 10, 10]) assert p.data[i, 10, 10, 10] == i - 1 for i in range(12, 20): assert np.all(p.data[i] == 0) diff --git a/tests/test_dse.py b/tests/test_dse.py index 0bb040982c..2904597128 100644 --- a/tests/test_dse.py +++ b/tests/test_dse.py @@ -2629,7 +2629,7 @@ def test_issue_2163(self): def test_dtype_aliases(self): a = np.arange(64).reshape((8, 8)) - grid = Grid(shape=a.shape, extent=(8, 8)) + grid = Grid(shape=a.shape, extent=(7, 7)) so = 2 f = Function(name='f', grid=grid, space_order=so, dtype=np.int32) @@ -2640,7 +2640,7 @@ def test_dtype_aliases(self): op.apply() assert FindNodes(Expression).visit(op)[0].dtype == np.float32 - assert np.all(fo.data[:-1, :-1] == 6) + assert np.all(fo.data[:-1, :-1] == 8) class TestIsoAcoustic(object): @@ -2685,13 +2685,13 @@ def test_fullopt(self): bns, _ = assert_blocking(op1, {'x0_blk0'}) # due to loop blocking assert summary0[('section0', None)].ops == 50 - assert summary0[('section1', None)].ops == 50 + assert summary0[('section1', None)].ops == 44 assert np.isclose(summary0[('section0', None)].oi, 2.851, atol=0.001) assert summary1[('section0', None)].ops == 9 assert summary1[('section1', None)].ops == 9 assert summary1[('section2', None)].ops == 31 - assert summary1[('section3', None)].ops == 32 + assert summary1[('section3', None)].ops == 26 assert np.isclose(summary1[('section2', None)].oi, 1.767, atol=0.001) assert np.allclose(u0.data, u1.data, atol=10e-5) diff --git a/tests/test_gradient.py b/tests/test_gradient.py index 9c91138c84..5624c5d461 100644 --- a/tests/test_gradient.py +++ b/tests/test_gradient.py @@ -15,6 +15,7 @@ class TestGradient(object): @skipif(['chkpnt', 'cpu64-icc']) + @switchconfig(safe_math=True) @pytest.mark.parametrize('dtype', [np.float32, np.float64]) @pytest.mark.parametrize('opt', [('advanced', {'openmp': True}), ('noop', {'openmp': True})]) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 186ab6cffd..dca94c8f40 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -84,25 +84,31 @@ def custom_points(grid, ranges, npoints, name='points'): return points -def precompute_linear_interpolation(points, grid, origin): - """ Sample precompute function that, given point and grid information - precomputes gridpoints and interpolation coefficients according to a linear - scheme to be used in PrecomputedSparseFunction. +def precompute_linear_interpolation(points, grid, origin, r=2): + """ + Sample precompute function that, given point and grid information + precomputes gridpoints and interpolation coefficients according to a linear + scheme to be used in PrecomputedSparseFunction. + + Allow larger radius with zero weights for testing. """ gridpoints = [tuple(floor((point[i]-origin[i])/grid.spacing[i]) for i in range(len(point))) for point in points] - interpolation_coeffs = np.zeros((len(points), 2, 2)) + interpolation_coeffs = np.zeros((len(points), grid.dim, r)) + rs = r // 2 - 1 for i, point in enumerate(points): for d in range(grid.dim): - interpolation_coeffs[i, d, 0] = ((gridpoints[i][d] + 1)*grid.spacing[d] - - point[d])/grid.spacing[d] - interpolation_coeffs[i, d, 1] = (point[d]-gridpoints[i][d]*grid.spacing[d])\ + gd = gridpoints[i][d] + interpolation_coeffs[i, d, rs] = ((gd + 1)*grid.spacing[d] - + point[d])/grid.spacing[d] + interpolation_coeffs[i, d, rs+1] = (point[d]-gd*grid.spacing[d])\ / grid.spacing[d] return gridpoints, interpolation_coeffs -def test_precomputed_interpolation(): +@pytest.mark.parametrize('r', [2, 4, 6]) +def test_precomputed_interpolation(r): """ Test interpolation with PrecomputedSparseFunction which accepts precomputed values for interpolation coefficients """ @@ -123,7 +129,8 @@ def init(data): m = Function(name='m', grid=grid, initializer=init, space_order=0) gridpoints, interpolation_coeffs = precompute_linear_interpolation(points, - grid, origin) + grid, origin, + r=r) sf = PrecomputedSparseFunction(name='s', grid=grid, r=r, npoint=len(points), gridpoints=gridpoints, @@ -136,7 +143,8 @@ def init(data): assert(all(np.isclose(sf.data, expected_values, rtol=1e-6))) -def test_precomputed_interpolation_time(): +@pytest.mark.parametrize('r', [2, 4, 6]) +def test_precomputed_interpolation_time(r): """ Test interpolation with PrecomputedSparseFunction which accepts precomputed values for interpolation coefficients, but this time with a TimeFunction @@ -154,7 +162,8 @@ def test_precomputed_interpolation_time(): u.data[it, :] = it gridpoints, interpolation_coeffs = precompute_linear_interpolation(points, - grid, origin) + grid, origin, + r=r) sf = PrecomputedSparseTimeFunction(name='s', grid=grid, r=r, npoint=len(points), nt=5, gridpoints=gridpoints, @@ -171,7 +180,8 @@ def test_precomputed_interpolation_time(): assert np.allclose(sf.data[it, :], it) -def test_precomputed_injection(): +@pytest.mark.parametrize('r', [2, 4, 6]) +def test_precomputed_injection(r): """Test injection with PrecomputedSparseFunction which accepts precomputed values for interpolation coefficients """ @@ -188,7 +198,8 @@ def test_precomputed_injection(): m.data[:] = 0. gridpoints, interpolation_coeffs = precompute_linear_interpolation(coords, - m.grid, origin) + m.grid, origin, + r=r) sf = PrecomputedSparseFunction(name='s', grid=m.grid, r=r, npoint=len(coords), gridpoints=gridpoints, @@ -206,7 +217,8 @@ def test_precomputed_injection(): assert np.allclose(m.data[indices], result, rtol=1.e-5) -def test_precomputed_injection_time(): +@pytest.mark.parametrize('r', [2, 4, 6]) +def test_precomputed_injection_time(r): """Test injection with PrecomputedSparseFunction which accepts precomputed values for interpolation coefficients """ @@ -224,7 +236,8 @@ def test_precomputed_injection_time(): m.data[:] = 0. gridpoints, interpolation_coeffs = precompute_linear_interpolation(coords, - m.grid, origin) + m.grid, origin, + r=r) sf = PrecomputedSparseTimeFunction(name='s', grid=m.grid, r=r, npoint=len(coords), gridpoints=gridpoints, nt=nt, @@ -718,7 +731,7 @@ class SparseFirst(SparseFunction): # No time dependence so need the implicit dim rec = s.interpolate(expr=s+fs, implicit_dims=grid.stepping_dim) op = Operator(eqs + rec) - print(op) + op(time_M=10) expected = 10*11/2 # n (n+1)/2 assert np.allclose(s.data, expected)