From a3c2426ffd9e4326f6d7dd738dd8d9029f43f427 Mon Sep 17 00:00:00 2001 From: mloubout Date: Wed, 28 Jan 2026 13:12:57 -0500 Subject: [PATCH 1/2] compiler: add rudimentary support for multi-cond buffering --- devito/ir/clusters/algorithms.py | 5 +- devito/ir/equations/equation.py | 2 +- devito/types/dimension.py | 11 +- .../userapi/05_conditional_dimension.ipynb | 187 ++++++++++-------- tests/test_buffering.py | 36 +++- 5 files changed, 153 insertions(+), 88 deletions(-) diff --git a/devito/ir/clusters/algorithms.py b/devito/ir/clusters/algorithms.py index 4223e789e7..788c31060b 100644 --- a/devito/ir/clusters/algorithms.py +++ b/devito/ir/clusters/algorithms.py @@ -254,6 +254,7 @@ def guard(clusters): # Chain together all `cds` conditions from all expressions in `c` guards = {} + mode = sympy.Or for cd in cds: # `BOTTOM` parent implies a guard that lives outside of # any iteration space, which corresponds to the placeholder None @@ -270,6 +271,7 @@ def guard(clusters): # Pull `cd` from any expr condition = guards.setdefault(k, []) + mode = mode and cd.relation for e in exprs: try: condition.append(e.conditionals[cd]) @@ -284,7 +286,8 @@ def guard(clusters): conditionals.pop(cd, None) exprs[i] = e.func(*e.args, conditionals=conditionals) - guards = {d: sympy.And(*v, evaluate=False) for d, v in guards.items()} + # Combination mode is And by default and Or if all conditions are + guards = {d: mode(*v, evaluate=False) for d, v in guards.items()} # Construct a guarded Cluster processed.append(c.rebuild(exprs=exprs, guards=guards)) diff --git a/devito/ir/equations/equation.py b/devito/ir/equations/equation.py index 29945903a9..8f7d35e155 100644 --- a/devito/ir/equations/equation.py +++ b/devito/ir/equations/equation.py @@ -223,7 +223,7 @@ def __new__(cls, *args, **kwargs): else: cond = diff2sympy(lower_exprs(d.condition)) if d._factor is not None: - cond = sympy.And(cond, GuardFactor(d)) + cond = d.relation(cond, GuardFactor(d)) conditionals[d] = cond # Replace dimension with index index = d.index diff --git a/devito/types/dimension.py b/devito/types/dimension.py index fa02ebb32d..ea07f5cff1 100644 --- a/devito/types/dimension.py +++ b/devito/types/dimension.py @@ -860,6 +860,8 @@ class ConditionalDimension(DerivedDimension): If True, use `self`, rather than the parent Dimension, to index into arrays. A typical use case is when arrays are accessed indirectly via the ``condition`` expression. + relation: Or/And, default=And + How this ConditionalDimension will be combined with other ones. Examples -------- @@ -913,10 +915,10 @@ class ConditionalDimension(DerivedDimension): is_Conditional = True __rkwargs__ = DerivedDimension.__rkwargs__ + \ - ('factor', 'condition', 'indirect') + ('factor', 'condition', 'indirect', 'relation') def __init_finalize__(self, name, parent=None, factor=None, condition=None, - indirect=False, **kwargs): + indirect=False, relation=sympy.And, **kwargs): # `parent=None` degenerates to a ConditionalDimension outside of # any iteration space if parent is None: @@ -937,6 +939,7 @@ def __init_finalize__(self, name, parent=None, factor=None, condition=None, self._condition = condition self._indirect = indirect + self._relation = relation @property def uses_symbolic_factor(self): @@ -978,6 +981,10 @@ def condition(self): def indirect(self): return self._indirect + @property + def relation(self): + return self._relation + @cached_property def free_symbols(self): retval = set(super().free_symbols) diff --git a/examples/userapi/05_conditional_dimension.ipynb b/examples/userapi/05_conditional_dimension.ipynb index 7f5e75fed5..e38fa989d8 100644 --- a/examples/userapi/05_conditional_dimension.ipynb +++ b/examples/userapi/05_conditional_dimension.ipynb @@ -71,6 +71,7 @@ "name": "stderr", "output_type": "stream", "text": [ + "NUMA domain count autodetection failed, assuming 1\n", "Operator `Kernel` ran in 0.01 s\n" ] }, @@ -158,75 +159,77 @@ "output_type": "stream", "text": [ "\n", - " Symbol defining a non-convex iteration sub-space derived from a ``parent``\n", - " Dimension, implemented by the compiler generating conditional \"if-then\" code\n", - " within the parent Dimension's iteration space.\n", + "Symbol defining a non-convex iteration sub-space derived from a ``parent``\n", + "Dimension, implemented by the compiler generating conditional \"if-then\" code\n", + "within the parent Dimension's iteration space.\n", "\n", - " Parameters\n", - " ----------\n", - " name : str\n", - " Name of the dimension.\n", - " parent : Dimension\n", - " The parent Dimension.\n", - " factor : int, optional, default=None\n", - " The number of iterations between two executions of the if-branch. If None\n", - " (default), ``condition`` must be provided.\n", - " condition : expr-like, optional, default=None\n", - " An arbitrary SymPy expression, typically involving the ``parent``\n", - " Dimension. When it evaluates to True, the if-branch is executed. If None\n", - " (default), ``factor`` must be provided.\n", - " indirect : bool, optional, default=False\n", - " If True, use `self`, rather than the parent Dimension, to\n", - " index into arrays. A typical use case is when arrays are accessed\n", - " indirectly via the ``condition`` expression.\n", + "Parameters\n", + "----------\n", + "name : str\n", + " Name of the dimension.\n", + "parent : Dimension\n", + " The parent Dimension.\n", + "factor : int, optional, default=None\n", + " The number of iterations between two executions of the if-branch. If None\n", + " (default), ``condition`` must be provided.\n", + "condition : expr-like, optional, default=None\n", + " An arbitrary SymPy expression, typically involving the ``parent``\n", + " Dimension. When it evaluates to True, the if-branch is executed. If None\n", + " (default), ``factor`` must be provided.\n", + "indirect : bool, optional, default=False\n", + " If True, use `self`, rather than the parent Dimension, to\n", + " index into arrays. A typical use case is when arrays are accessed\n", + " indirectly via the ``condition`` expression.\n", + "relation: Or/And, default=And\n", + " How this ConditionalDimension will be combined with other ones.\n", "\n", - " Examples\n", - " --------\n", - " Among the other things, ConditionalDimensions are indicated to implement\n", - " Function subsampling. In the following example, an Operator evaluates the\n", - " Function ``g`` and saves its content into ``f`` every ``factor=4`` iterations.\n", + "Examples\n", + "--------\n", + "Among the other things, ConditionalDimensions are indicated to implement\n", + "Function subsampling. In the following example, an Operator evaluates the\n", + "Function ``g`` and saves its content into ``f`` every ``factor=4`` iterations.\n", "\n", - " >>> from devito import Dimension, ConditionalDimension, Function, Eq, Operator\n", - " >>> size, factor = 16, 4\n", - " >>> i = Dimension(name='i')\n", - " >>> ci = ConditionalDimension(name='ci', parent=i, factor=factor)\n", - " >>> g = Function(name='g', shape=(size,), dimensions=(i,))\n", - " >>> f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))\n", - " >>> op = Operator([Eq(g, 1), Eq(f, g)])\n", + ">>> from devito import Dimension, ConditionalDimension, Function, Eq, Operator\n", + ">>> size, factor = 16, 4\n", + ">>> i = Dimension(name='i')\n", + ">>> ci = ConditionalDimension(name='ci', parent=i, factor=factor)\n", + ">>> g = Function(name='g', shape=(size,), dimensions=(i,))\n", + ">>> f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))\n", + ">>> op = Operator([Eq(g, 1), Eq(f, g)])\n", "\n", - " The Operator generates the following for-loop (pseudocode)\n", + "The Operator generates the following for-loop (pseudocode)\n", "\n", - " .. code-block:: C\n", + ".. code-block:: C\n", "\n", - " for (int i = i_m; i <= i_M; i += 1) {\n", - " g[i] = 1;\n", - " if (i%4 == 0) {\n", - " f[i / 4] = g[i];\n", - " }\n", - " }\n", + " for (int i = i_m; i <= i_M; i += 1) {\n", + " g[i] = 1;\n", + " if (i%4 == 0) {\n", + " f[i / 4] = g[i];\n", + " }\n", + " }\n", "\n", - " Another typical use case is when one needs to constrain the execution of\n", - " loop iterations so that certain conditions are honoured. The following\n", - " artificial example uses ConditionalDimension to guard against out-of-bounds\n", - " accesses in indirectly accessed arrays.\n", + "Another typical use case is when one needs to constrain the execution of\n", + "loop iterations so that certain conditions are honoured. The following\n", + "artificial example uses ConditionalDimension to guard against out-of-bounds\n", + "accesses in indirectly accessed arrays.\n", "\n", - " >>> from sympy import And\n", - " >>> ci = ConditionalDimension(name='ci', parent=i,\n", - " ... condition=And(g[i] > 0, g[i] < 4, evaluate=False))\n", - " >>> f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))\n", - " >>> op = Operator(Eq(f[g[i]], f[g[i]] + 1))\n", + ">>> from sympy import And\n", + ">>> ci = ConditionalDimension(name='ci', parent=i,\n", + "... condition=And(g[i] > 0, g[i] < 4, evaluate=False))\n", + ">>> f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))\n", + ">>> op = Operator(Eq(f[g[i]], f[g[i]] + 1))\n", "\n", - " The Operator generates the following for-loop (pseudocode)\n", + "The Operator generates the following for-loop (pseudocode)\n", "\n", - " .. code-block:: C\n", + ".. code-block:: C\n", "\n", - " for (int i = i_m; i <= i_M; i += 1) {\n", - " if (g[i] > 0 && g[i] < 4) {\n", - " f[g[i]] = f[g[i]] + 1;\n", - " }\n", - " }\n", + " for (int i = i_m; i <= i_M; i += 1) {\n", + " if (g[i] > 0 && g[i] < 4) {\n", + " f[g[i]] = f[g[i]] + 1;\n", + " }\n", + " }\n", "\n", - " \n" + "\n" ] } ], @@ -248,9 +251,7 @@ { "cell_type": "code", "execution_count": 5, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [ { "name": "stderr", @@ -321,14 +322,18 @@ "output_type": "stream", "text": [ "START(section0)\n", - "for (int x = x_m; x <= x_M; x += 1)\n", + "#pragma omp parallel num_threads(nthreads)\n", "{\n", - " #pragma omp simd aligned(f:32)\n", - " for (int y = y_m; y <= y_M; y += 1)\n", + " #pragma omp for schedule(static,1)\n", + " for (int x = x_m; x <= x_M; x += 1)\n", " {\n", - " if (f[x + 1][y + 1] > 0)\n", + " #pragma omp simd aligned(f:16)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " f[x + 1][y + 1] = f[x + 1][y + 1] + 1;\n", + " if (f[x + 1][y + 1] > 0)\n", + " {\n", + " f[x + 1][y + 1] = f[x + 1][y + 1] + 1;\n", + " }\n", " }\n", " }\n", "}\n", @@ -397,14 +402,18 @@ "output_type": "stream", "text": [ "START(section0)\n", - "for (int x = x_m; x <= x_M; x += 1)\n", + "#pragma omp parallel num_threads(nthreads)\n", "{\n", - " #pragma omp simd aligned(f,g:32)\n", - " for (int y = y_m; y <= y_M; y += 1)\n", + " #pragma omp for schedule(static,1)\n", + " for (int x = x_m; x <= x_M; x += 1)\n", " {\n", - " if (y < 5 && g[x + 1][y + 1] != 0)\n", + " #pragma omp simd aligned(f,g:16)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " f[x + 1][y + 1] = f[x + 1][y + 1] + g[x + 1][y + 1];\n", + " if (y < 5 && g[x + 1][y + 1] != 0)\n", + " {\n", + " f[x + 1][y + 1] = f[x + 1][y + 1] + g[x + 1][y + 1];\n", + " }\n", " }\n", " }\n", "}\n", @@ -489,13 +498,17 @@ "output_type": "stream", "text": [ "START(section0)\n", - "for (int x = x_m; x <= x_M; x += 1)\n", + "#pragma omp parallel num_threads(nthreads)\n", "{\n", - " for (int y = y_m; y <= y_M; y += 1)\n", + " #pragma omp for schedule(static,1)\n", + " for (int x = x_m; x <= x_M; x += 1)\n", " {\n", - " if (y < 5 && g[x + 1][y + 1] != 0)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " h[x + 1][y + 1] = g[x + 1][y + 1] + h[x + 1][y + 1];\n", + " if (y < 5 && g[x + 1][y + 1] != 0)\n", + " {\n", + " h[x + 1][y + 1] = g[x + 1][y + 1] + h[x + 1][y + 1];\n", + " }\n", " }\n", " }\n", "}\n", @@ -563,11 +576,15 @@ "output_type": "stream", "text": [ "START(section0)\n", - "for (int i = i_m; i <= i_M; i += 1)\n", + "#pragma omp parallel num_threads(nthreads)\n", "{\n", - " if ((i)%(cif) == 0)\n", + " #pragma omp for schedule(static,1)\n", + " for (int i = i_m; i <= i_M; i += 1)\n", " {\n", - " f[i / cif] = g[i];\n", + " if ((i)%(cif) == 0)\n", + " {\n", + " f[i / cif] = g[i];\n", + " }\n", " }\n", "}\n", "STOP(section0,timers)\n", @@ -678,13 +695,17 @@ "output_type": "stream", "text": [ "START(section0)\n", - "for (int x = x_m; x <= x_M; x += 1)\n", + "#pragma omp parallel num_threads(nthreads)\n", "{\n", - " for (int y = y_m; y <= y_M; y += 1)\n", + " #pragma omp for collapse(2) schedule(static,1) reduction(+:g[1])\n", + " for (int x = x_m; x <= x_M; x += 1)\n", " {\n", - " if (f[x][y] != 0)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " g[1] += 1;\n", + " if (f[x][y] != 0)\n", + " {\n", + " g[1] += 1;\n", + " }\n", " }\n", " }\n", "}\n", @@ -694,7 +715,7 @@ { "data": { "text/plain": [ - "10" + "np.int32(10)" ] }, "execution_count": 11, @@ -831,7 +852,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.2" + "version": "3.13.11" } }, "nbformat": 4, diff --git a/tests/test_buffering.py b/tests/test_buffering.py index 2959796851..5cc5c5b440 100644 --- a/tests/test_buffering.py +++ b/tests/test_buffering.py @@ -1,9 +1,10 @@ import numpy as np import pytest +from sympy import Or from conftest import skipif from devito import ( - ConditionalDimension, Constant, Eq, Grid, Operator, SubDimension, SubDomain, + CondEq, ConditionalDimension, Constant, Eq, Grid, Operator, SubDimension, SubDomain, TimeFunction, configuration, switchconfig ) from devito.arch.archinfo import AppleArm @@ -751,3 +752,36 @@ def test_buffer_reuse(): assert all(np.all(usave.data[i-1] == i) for i in range(1, nt + 1)) assert all(np.all(vsave.data[i-1] == i + 1) for i in range(1, nt + 1)) + + +def test_multi_cond(): + grid = Grid((3, 3)) + nt = 5 + + x, y = grid.dimensions + + factor = 2 + ntmod = (nt - 1) * factor + 1 + + ct1 = ConditionalDimension(name="ct1", parent=grid.time_dim, + factor=factor, relation=Or) + ctend = ConditionalDimension(name="ctend", parent=grid.time_dim, + condition=CondEq(grid.time_dim, ntmod - 2), + relation=Or) + + f = TimeFunction(grid=grid, name='f', time_order=0, + space_order=0, save=nt, time_dim=ct1) + T = TimeFunction(grid=grid, name='T', time_order=0, space_order=0) + + eqs = [Eq(T, grid.time_dim)] + # this to save times from 0 to nt - 2 + eqs.append(Eq(f, T)) + # this to save the last time sample nt - 1 + eqs.append(Eq(f.forward, T+1, implicit_dims=ctend)) + + # run operator with buffering + op = Operator(eqs, opt=('streaming', 'buffering')) + op.apply(time_m=0, time_M=ntmod-2) + + for i in range(nt): + assert np.allclose(f.data[i], i*2) From 8304fdb2897af384826478f85f05fefe98899e4e Mon Sep 17 00:00:00 2001 From: mloubout Date: Thu, 29 Jan 2026 13:00:33 -0500 Subject: [PATCH 2/2] CI: misc docker deploy fixes --- .github/workflows/docker-bases.yaml | 12 +- devito/ir/clusters/algorithms.py | 3 +- devito/types/dimension.py | 5 +- docker/Dockerfile.intel | 4 +- docker/Dockerfile.nvidia | 4 +- examples/seismic/tutorials/02_rtm.ipynb | 2 +- .../userapi/05_conditional_dimension.ipynb | 418 +++++++++++++++--- tests/test_buffering.py | 2 +- 8 files changed, 365 insertions(+), 85 deletions(-) diff --git a/.github/workflows/docker-bases.yaml b/.github/workflows/docker-bases.yaml index 48ada2ba8f..56b6bf0860 100644 --- a/.github/workflows/docker-bases.yaml +++ b/.github/workflows/docker-bases.yaml @@ -40,7 +40,7 @@ jobs: ############## Basic gcc CPU ########################## ####################################################### deploy-cpu-bases: - if: inputs.cpu + if: ${{ github.event_name != 'workflow_dispatch' || inputs.cpu }} name: "cpu-base-${{ matrix.arch }}-gcc${{ matrix.gcc }}" runs-on: ${{ matrix.runner }} env: @@ -88,7 +88,7 @@ jobs: tags: "devitocodes/bases:cpu-gcc${{ matrix.gcc }}-${{ matrix.arch }}" deploy-cpu-bases-manifest: - if: inputs.cpu + if: ${{ github.event_name != 'workflow_dispatch' || inputs.cpu }} name: "cpu-base-manifest" runs-on: ubuntu-latest needs: deploy-cpu-bases @@ -123,7 +123,7 @@ jobs: ############## Intel OneApi CPU ####################### ####################################################### deploy-oneapi-bases: - if: inputs.intel + if: ${{ github.event_name != 'workflow_dispatch' || inputs.intel }} name: "oneapi-base" runs-on: ubuntu-latest env: @@ -181,7 +181,7 @@ jobs: ################### Nvidia nvhpc ###################### ####################################################### deploy-nvidia-bases: - if: inputs.nvidia + if: ${{ github.event_name != 'workflow_dispatch' || inputs.nvidia }} name: "nvidia-bases-${{ matrix.arch }}" runs-on: ${{ matrix.runner }} env: @@ -270,7 +270,7 @@ jobs: tags: "devitocodes/bases:cpu-nvc${{ matrix.extra_tag }}-${{ matrix.arch }}" deploy-nvidia-bases-manifest: - if: inputs.nvidia + if: ${{ github.event_name != 'workflow_dispatch' || inputs.nvidia }} name: "nvidia-base-manifest" runs-on: ubuntu-latest needs: deploy-nvidia-bases @@ -305,7 +305,7 @@ jobs: ##################### AMD ############################# ####################################################### deploy-amd-bases: - if: inputs.amd + if: ${{ github.event_name != 'workflow_dispatch' || inputs.amd }} name: "amd-base" runs-on: ["self-hosted", "amdgpu"] env: diff --git a/devito/ir/clusters/algorithms.py b/devito/ir/clusters/algorithms.py index 788c31060b..ab688e3c81 100644 --- a/devito/ir/clusters/algorithms.py +++ b/devito/ir/clusters/algorithms.py @@ -286,7 +286,8 @@ def guard(clusters): conditionals.pop(cd, None) exprs[i] = e.func(*e.args, conditionals=conditionals) - # Combination mode is And by default and Or if all conditions are + # Combination `mode` is And by default. + # If all conditions are Or then Or combination `mode` is used. guards = {d: mode(*v, evaluate=False) for d, v in guards.items()} # Construct a guarded Cluster diff --git a/devito/types/dimension.py b/devito/types/dimension.py index ea07f5cff1..4010472bef 100644 --- a/devito/types/dimension.py +++ b/devito/types/dimension.py @@ -861,7 +861,10 @@ class ConditionalDimension(DerivedDimension): index into arrays. A typical use case is when arrays are accessed indirectly via the ``condition`` expression. relation: Or/And, default=And - How this ConditionalDimension will be combined with other ones. + How this ConditionalDimension will be combined with other ones during + lowering for example combining Function's ConditionalDimension with + an Equation's implicit_dim. All Dimensions within an equation + must have `Or` relation for the final combined condition to be Or. Examples -------- diff --git a/docker/Dockerfile.intel b/docker/Dockerfile.intel index a84f20fed7..c7a7a0f071 100644 --- a/docker/Dockerfile.intel +++ b/docker/Dockerfile.intel @@ -26,7 +26,7 @@ FROM base AS oneapi # Download the key to system keyring # https://www.intel.com/content/www/us/en/develop/documentation/installation-guide-for-intel-oneapi-toolkits-linux/top/installation/install-using-package-managers/apt.html#apt -SHELL /bin/bash -o pipefail +SHELL ["/bin/bash", "-o", "pipefail", "-c"] RUN wget --progress=dot:giga -O- https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB | gpg --dearmor > /usr/share/keyrings/oneapi-archive-keyring.gpg RUN echo "deb [signed-by=/usr/share/keyrings/oneapi-archive-keyring.gpg] https://apt.repos.intel.com/oneapi all main" > /etc/apt/sources.list.d/oneAPI.list @@ -37,7 +37,7 @@ RUN apt-get update -y && \ # Drivers mandatory for intel gpu # https://dgpu-docs.intel.com/driver/installation.html#ubuntu-install-steps -SHELL /bin/bash -o pipefail +SHELL ["/bin/bash", "-o", "pipefail", "-c"] RUN wget -qO - https://repositories.intel.com/graphics/intel-graphics.key | gpg --dearmor > /usr/share/keyrings/intel-graphics.gpg RUN echo "deb [arch=amd64 signed-by=/usr/share/keyrings/intel-graphics.gpg] https://repositories.intel.com/graphics/ubuntu jammy unified" > /etc/apt/sources.list.d/intel-gpu-jammy.list diff --git a/docker/Dockerfile.nvidia b/docker/Dockerfile.nvidia index 00a86ca96a..bf14afa451 100644 --- a/docker/Dockerfile.nvidia +++ b/docker/Dockerfile.nvidia @@ -19,7 +19,7 @@ RUN apt-get update && \ dh-autoreconf python3-venv python3-dev python3-pip # nodesource: nvdashboard requires nodejs>=10 -SHELL /bin/bash -o pipefail +SHELL ["/bin/bash", "-o", "pipefail", "-c"] RUN curl https://developer.download.nvidia.com/hpc-sdk/ubuntu/DEB-GPG-KEY-NVIDIA-HPC-SDK | gpg --yes --dearmor -o /usr/share/keyrings/nvidia-hpcsdk-archive-keyring.gpg RUN arch="$(uname -m)" && \ case "$arch" in \ @@ -92,7 +92,7 @@ ENV UCX_TLS=cuda,cuda_copy,cuda_ipc,sm,shm,self #ENV UCX_TLS=cuda,cuda_copy,cuda_ipc,sm,shm,self,rc_x,gdr_copy # Make simlink for path setup since ENV doesn't accept shell commands. -SHELL /bin/bash -o pipefail +SHELL ["/bin/bash", "-o", "pipefail", "-c"] RUN arch="$(uname -m)" && \ case "$arch" in \ x86_64) linux=Linux_x86_64 ;; \ diff --git a/examples/seismic/tutorials/02_rtm.ipynb b/examples/seismic/tutorials/02_rtm.ipynb index cf90a02585..d4554fd285 100644 --- a/examples/seismic/tutorials/02_rtm.ipynb +++ b/examples/seismic/tutorials/02_rtm.ipynb @@ -87,7 +87,7 @@ " return demo_model('marmousi2d-isotropic', data_path='../../../../data/',\n", " grid=grid, nbl=20)\n", " filter_sigma = (6, 6)\n", - " nshots = 301 # Need good covergae in shots, one every two grid points\n", + " nshots = 301 # Need good coverage in shots, one every two grid points\n", " nreceivers = 601 # One receiver every grid point\n", " t0 = 0.\n", " tn = 3500. # Simulation last 3.5 second (3500 ms)\n", diff --git a/examples/userapi/05_conditional_dimension.ipynb b/examples/userapi/05_conditional_dimension.ipynb index e38fa989d8..91cb9cdda8 100644 --- a/examples/userapi/05_conditional_dimension.ipynb +++ b/examples/userapi/05_conditional_dimension.ipynb @@ -71,7 +71,6 @@ "name": "stderr", "output_type": "stream", "text": [ - "NUMA domain count autodetection failed, assuming 1\n", "Operator `Kernel` ran in 0.01 s\n" ] }, @@ -159,77 +158,80 @@ "output_type": "stream", "text": [ "\n", - "Symbol defining a non-convex iteration sub-space derived from a ``parent``\n", - "Dimension, implemented by the compiler generating conditional \"if-then\" code\n", - "within the parent Dimension's iteration space.\n", - "\n", - "Parameters\n", - "----------\n", - "name : str\n", - " Name of the dimension.\n", - "parent : Dimension\n", - " The parent Dimension.\n", - "factor : int, optional, default=None\n", - " The number of iterations between two executions of the if-branch. If None\n", - " (default), ``condition`` must be provided.\n", - "condition : expr-like, optional, default=None\n", - " An arbitrary SymPy expression, typically involving the ``parent``\n", - " Dimension. When it evaluates to True, the if-branch is executed. If None\n", - " (default), ``factor`` must be provided.\n", - "indirect : bool, optional, default=False\n", - " If True, use `self`, rather than the parent Dimension, to\n", - " index into arrays. A typical use case is when arrays are accessed\n", - " indirectly via the ``condition`` expression.\n", - "relation: Or/And, default=And\n", - " How this ConditionalDimension will be combined with other ones.\n", - "\n", - "Examples\n", - "--------\n", - "Among the other things, ConditionalDimensions are indicated to implement\n", - "Function subsampling. In the following example, an Operator evaluates the\n", - "Function ``g`` and saves its content into ``f`` every ``factor=4`` iterations.\n", - "\n", - ">>> from devito import Dimension, ConditionalDimension, Function, Eq, Operator\n", - ">>> size, factor = 16, 4\n", - ">>> i = Dimension(name='i')\n", - ">>> ci = ConditionalDimension(name='ci', parent=i, factor=factor)\n", - ">>> g = Function(name='g', shape=(size,), dimensions=(i,))\n", - ">>> f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))\n", - ">>> op = Operator([Eq(g, 1), Eq(f, g)])\n", - "\n", - "The Operator generates the following for-loop (pseudocode)\n", - "\n", - ".. code-block:: C\n", - "\n", - " for (int i = i_m; i <= i_M; i += 1) {\n", - " g[i] = 1;\n", - " if (i%4 == 0) {\n", - " f[i / 4] = g[i];\n", - " }\n", - " }\n", + " Symbol defining a non-convex iteration sub-space derived from a ``parent``\n", + " Dimension, implemented by the compiler generating conditional \"if-then\" code\n", + " within the parent Dimension's iteration space.\n", "\n", - "Another typical use case is when one needs to constrain the execution of\n", - "loop iterations so that certain conditions are honoured. The following\n", - "artificial example uses ConditionalDimension to guard against out-of-bounds\n", - "accesses in indirectly accessed arrays.\n", + " Parameters\n", + " ----------\n", + " name : str\n", + " Name of the dimension.\n", + " parent : Dimension\n", + " The parent Dimension.\n", + " factor : int, optional, default=None\n", + " The number of iterations between two executions of the if-branch. If None\n", + " (default), ``condition`` must be provided.\n", + " condition : expr-like, optional, default=None\n", + " An arbitrary SymPy expression, typically involving the ``parent``\n", + " Dimension. When it evaluates to True, the if-branch is executed. If None\n", + " (default), ``factor`` must be provided.\n", + " indirect : bool, optional, default=False\n", + " If True, use `self`, rather than the parent Dimension, to\n", + " index into arrays. A typical use case is when arrays are accessed\n", + " indirectly via the ``condition`` expression.\n", + " relation: Or/And, default=And\n", + " How this ConditionalDimension will be combined with other ones during\n", + " lowering for example combining Function's ConditionalDimension with\n", + " an Equation's implicit_dim. All Dimensions within an equation\n", + " must have `Or` relation for the final combined condition to be Or.\n", "\n", - ">>> from sympy import And\n", - ">>> ci = ConditionalDimension(name='ci', parent=i,\n", - "... condition=And(g[i] > 0, g[i] < 4, evaluate=False))\n", - ">>> f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))\n", - ">>> op = Operator(Eq(f[g[i]], f[g[i]] + 1))\n", + " Examples\n", + " --------\n", + " Among the other things, ConditionalDimensions are indicated to implement\n", + " Function subsampling. In the following example, an Operator evaluates the\n", + " Function ``g`` and saves its content into ``f`` every ``factor=4`` iterations.\n", "\n", - "The Operator generates the following for-loop (pseudocode)\n", + " >>> from devito import Dimension, ConditionalDimension, Function, Eq, Operator\n", + " >>> size, factor = 16, 4\n", + " >>> i = Dimension(name='i')\n", + " >>> ci = ConditionalDimension(name='ci', parent=i, factor=factor)\n", + " >>> g = Function(name='g', shape=(size,), dimensions=(i,))\n", + " >>> f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))\n", + " >>> op = Operator([Eq(g, 1), Eq(f, g)])\n", "\n", - ".. code-block:: C\n", + " The Operator generates the following for-loop (pseudocode)\n", "\n", - " for (int i = i_m; i <= i_M; i += 1) {\n", - " if (g[i] > 0 && g[i] < 4) {\n", - " f[g[i]] = f[g[i]] + 1;\n", - " }\n", - " }\n", + " .. code-block:: C\n", "\n", - "\n" + " for (int i = i_m; i <= i_M; i += 1) {\n", + " g[i] = 1;\n", + " if (i%4 == 0) {\n", + " f[i / 4] = g[i];\n", + " }\n", + " }\n", + "\n", + " Another typical use case is when one needs to constrain the execution of\n", + " loop iterations so that certain conditions are honoured. The following\n", + " artificial example uses ConditionalDimension to guard against out-of-bounds\n", + " accesses in indirectly accessed arrays.\n", + "\n", + " >>> from sympy import And\n", + " >>> ci = ConditionalDimension(name='ci', parent=i,\n", + " ... condition=And(g[i] > 0, g[i] < 4, evaluate=False))\n", + " >>> f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))\n", + " >>> op = Operator(Eq(f[g[i]], f[g[i]] + 1))\n", + "\n", + " The Operator generates the following for-loop (pseudocode)\n", + "\n", + " .. code-block:: C\n", + "\n", + " for (int i = i_m; i <= i_M; i += 1) {\n", + " if (g[i] > 0 && g[i] < 4) {\n", + " f[g[i]] = f[g[i]] + 1;\n", + " }\n", + " }\n", + "\n", + " \n" ] } ], @@ -327,7 +329,7 @@ " #pragma omp for schedule(static,1)\n", " for (int x = x_m; x <= x_M; x += 1)\n", " {\n", - " #pragma omp simd aligned(f:16)\n", + " #pragma omp simd aligned(f:64)\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", " if (f[x + 1][y + 1] > 0)\n", @@ -407,7 +409,7 @@ " #pragma omp for schedule(static,1)\n", " for (int x = x_m; x <= x_M; x += 1)\n", " {\n", - " #pragma omp simd aligned(f,g:16)\n", + " #pragma omp simd aligned(f,g:64)\n", " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", " if (y < 5 && g[x + 1][y + 1] != 0)\n", @@ -715,7 +717,7 @@ { "data": { "text/plain": [ - "np.int32(10)" + "10" ] }, "execution_count": 11, @@ -825,6 +827,280 @@ " assert all(g.data[:, i] == ids_np[i])" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Example C: Combining ConditionaDimension\n", + "\n", + "In some cases, a `ConditionaDimension` might be used in combination with an implicit_dim to handle specific cases. This combination can be made mutually exclusive (And) or inclusive (Or).\n", + "\n", + "As an example, let's consider the following case:\n", + "\n", + "- Set all even x indices to 1 using the standard subsampling `factor` ConditionaDimension\n", + "- Set the edges to 2" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "from sympy import Or\n", + "from devito import CondEq\n", + "\n", + "x, y = grid.dimensions\n", + "c_even = ConditionalDimension(name='ceven', parent=x, factor=2, relation=Or)\n", + "c_edge = ConditionalDimension(name='cedge', parent=grid.dimensions[0], condition=CondEq(x, x.symbolic_max), relation=Or)\n", + "c_edge_a = ConditionalDimension(name='cedge', parent=grid.dimensions[0], condition=CondEq(x, x.symbolic_max-1))\n", + "\n", + "f = Function(name=\"f\", grid=grid, dimensions=(c_even, y), shape=(5, 10), space_order=0)\n", + "\n", + "op = Operator(Eq(f, 1, implicit_dims=c_edge))" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/* Devito generated code for Operator `Kernel` */\n", + "\n", + "#define _POSIX_C_SOURCE 200809L\n", + "#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);\n", + "#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;\n", + "\n", + "#include \"stdlib.h\"\n", + "#include \"math.h\"\n", + "#include \"sys/time.h\"\n", + "#include \"xmmintrin.h\"\n", + "#include \"pmmintrin.h\"\n", + "#include \"omp.h\"\n", + "\n", + "struct dataobj\n", + "{\n", + " void *restrict data;\n", + " int * size;\n", + " unsigned long nbytes;\n", + " unsigned long * npsize;\n", + " unsigned long * dsize;\n", + " int * hsize;\n", + " int * hofs;\n", + " int * oofs;\n", + " void * dmap;\n", + "} ;\n", + "\n", + "struct profiler\n", + "{\n", + " double section0;\n", + "} ;\n", + "\n", + "\n", + "int Kernel(const int cevenf, struct dataobj *restrict f_vec, const int x_M, const int x_m, const int y_M, const int y_m, const int nthreads, struct profiler * timers)\n", + "{\n", + " float (*restrict f)[f_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[f_vec->size[1]]) f_vec->data;\n", + "\n", + " /* Flush denormal numbers to zero in hardware */\n", + " _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);\n", + " _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);\n", + "\n", + " START(section0)\n", + " #pragma omp parallel num_threads(nthreads)\n", + " {\n", + " #pragma omp for schedule(static,1)\n", + " for (int x = x_m; x <= x_M; x += 1)\n", + " {\n", + " if (x == x_M || (x)%(cevenf) == 0)\n", + " {\n", + " #pragma omp simd aligned(f:64)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", + " {\n", + " f[x / cevenf][y] = 1;\n", + " }\n", + " }\n", + " }\n", + " }\n", + " STOP(section0,timers)\n", + "\n", + " return 0;\n", + "}\n", + "\n" + ] + }, + { + "data": { + "text/plain": [ + "PerformanceSummary([(PerfKey(name='section0', rank=None),\n", + " PerfEntry(time=0.000124, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "print(op)\n", + "op()" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n", + " [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]\n" + ] + } + ], + "source": [ + "print(f.data)\n", + "assert np.all(f.data == 1)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/* Devito generated code for Operator `Kernel` */\n", + "\n", + "#define _POSIX_C_SOURCE 200809L\n", + "#define START(S) struct timeval start_ ## S , end_ ## S ; gettimeofday(&start_ ## S , NULL);\n", + "#define STOP(S,T) gettimeofday(&end_ ## S, NULL); T->S += (double)(end_ ## S .tv_sec-start_ ## S.tv_sec)+(double)(end_ ## S .tv_usec-start_ ## S .tv_usec)/1000000;\n", + "\n", + "#include \"stdlib.h\"\n", + "#include \"math.h\"\n", + "#include \"sys/time.h\"\n", + "#include \"xmmintrin.h\"\n", + "#include \"pmmintrin.h\"\n", + "#include \"omp.h\"\n", + "\n", + "struct dataobj\n", + "{\n", + " void *restrict data;\n", + " int * size;\n", + " unsigned long nbytes;\n", + " unsigned long * npsize;\n", + " unsigned long * dsize;\n", + " int * hsize;\n", + " int * hofs;\n", + " int * oofs;\n", + " void * dmap;\n", + "} ;\n", + "\n", + "struct profiler\n", + "{\n", + " double section0;\n", + "} ;\n", + "\n", + "\n", + "int Kernel(const int cevenf, struct dataobj *restrict f_vec, const int x_M, const int x_m, const int y_M, const int y_m, const int nthreads, struct profiler * timers)\n", + "{\n", + " float (*restrict f)[f_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[f_vec->size[1]]) f_vec->data;\n", + "\n", + " /* Flush denormal numbers to zero in hardware */\n", + " _MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON);\n", + " _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);\n", + "\n", + " START(section0)\n", + " #pragma omp parallel num_threads(nthreads)\n", + " {\n", + " #pragma omp for schedule(static,1)\n", + " for (int x = x_m; x <= x_M; x += 1)\n", + " {\n", + " if (x == x_M - 1 && (x)%(cevenf) == 0)\n", + " {\n", + " #pragma omp simd aligned(f:64)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", + " {\n", + " f[x / cevenf][y] = 1;\n", + " }\n", + " }\n", + " }\n", + " }\n", + " STOP(section0,timers)\n", + "\n", + " return 0;\n", + "}\n", + "\n" + ] + }, + { + "data": { + "text/plain": [ + "PerformanceSummary([(PerfKey(name='section0', rank=None),\n", + " PerfEntry(time=0.000129, gflopss=0.0, gpointss=0.0, oi=0.0, ops=0, itershapes=[]))])" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# NBVAL_IGNORE_OUTPUT\n", + "op = Operator(Eq(f, 1, implicit_dims=c_edge_a))\n", + "print(op)\n", + "f.data.fill(0)\n", + "op()" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]\n", + " [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]\n" + ] + } + ], + "source": [ + "print(f.data)\n", + "assert np.all(f.data[:-1, :] == 0)\n", + "assert np.all(f.data[-1, :] == 1)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -838,7 +1114,7 @@ "metadata": { "celltoolbar": "Tags", "kernelspec": { - "display_name": "Python 3", + "display_name": "devitoenv (3.11.0)", "language": "python", "name": "python3" }, @@ -852,7 +1128,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.11" + "version": "3.11.0rc1" } }, "nbformat": 4, diff --git a/tests/test_buffering.py b/tests/test_buffering.py index 5cc5c5b440..0a0037c223 100644 --- a/tests/test_buffering.py +++ b/tests/test_buffering.py @@ -780,7 +780,7 @@ def test_multi_cond(): eqs.append(Eq(f.forward, T+1, implicit_dims=ctend)) # run operator with buffering - op = Operator(eqs, opt=('streaming', 'buffering')) + op = Operator(eqs, opt='buffering') op.apply(time_m=0, time_M=ntmod-2) for i in range(nt):