diff --git a/examples/userapi/05_conditional_dimension.ipynb b/examples/userapi/05_conditional_dimension.ipynb index eb23270a09..42edf13c85 100644 --- a/examples/userapi/05_conditional_dimension.ipynb +++ b/examples/userapi/05_conditional_dimension.ipynb @@ -71,7 +71,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "Operator `Kernel` run in 0.01 s\n" + "Operator `Kernel` ran in 0.01 s\n" ] }, { @@ -239,7 +239,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Use case A (honour a condition)\n", + "# Use case A - Honour a condition\n", "\n", "Now it's time to show a more descriptive example. We define a conditional that increments by one all values of a Function that are larger than zero. We name our `ConditionalDimension` `ci`.\n", "In this example we want to update again by one all the elements in `f` that are larger than zero. Before updating the elements we reinitialize our data to the eye function." @@ -256,7 +256,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "Operator `Kernel` run in 0.01 s\n" + "Operator `Kernel` ran in 0.01 s\n" ] }, { @@ -313,63 +313,28 @@ "name": "stderr", "output_type": "stream", "text": [ - "Operator `Kernel` run in 0.01 s\n" + "Operator `Kernel` ran in 0.01 s\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "#define _POSIX_C_SOURCE 200809L\n", - "#include \"stdlib.h\"\n", - "#include \"math.h\"\n", - "#include \"sys/time.h\"\n", - "#include \"xmmintrin.h\"\n", - "#include \"pmmintrin.h\"\n", - "\n", - "struct dataobj\n", - "{\n", - " void *restrict data;\n", - " int * size;\n", - " int * npsize;\n", - " int * dsize;\n", - " int * hsize;\n", - " int * hofs;\n", - " int * oofs;\n", - "} ;\n", - "\n", - "struct profiler\n", - "{\n", - " double section0;\n", - "} ;\n", - "\n", - "\n", - "int Kernel(struct dataobj *restrict f_vec, struct profiler * timers, const int x_M, const int x_m, const int y_M, const int y_m)\n", + "/* Begin section0 */\n", + "START_TIMER(section0)\n", + "for (int x = x_m; x <= x_M; x += 1)\n", "{\n", - " float (*restrict f)[f_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[f_vec->size[1]]) f_vec->data;\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", - " struct timeval start_section0, end_section0;\n", - " gettimeofday(&start_section0, NULL);\n", - " /* Begin section0 */\n", - " for (int x = x_m; x <= x_M; x += 1)\n", + " #pragma omp simd aligned(f:32)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " #pragma omp simd aligned(f:32)\n", - " for (int y = y_m; y <= y_M; y += 1)\n", + " if (f[x + 1][y + 1] > 0)\n", " {\n", - " if (f[x + 1][y + 1] > 0)\n", - " {\n", - " f[x + 1][y + 1] = f[x + 1][y + 1] + 1;\n", - " }\n", + " f[x + 1][y + 1] = f[x + 1][y + 1] + 1;\n", " }\n", " }\n", - " /* End section0 */\n", - " gettimeofday(&end_section0, NULL);\n", - " timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;\n", - " return 0;\n", "}\n", - "\n" + "STOP_TIMER(section0,timers)\n", + "/* End section0 */\n" ] }, { @@ -398,7 +363,7 @@ "f.data[:] = np.eye(10)\n", "\n", "op2 = Operator(Eq(f, f + 1, implicit_dims=ci))\n", - "print(op2.ccode)\n", + "print(op2.body.body[-1])\n", "op2.apply()\n", "\n", "assert (np.count_nonzero(f.data - np.diag(np.diagonal(f.data)))==0)\n", @@ -426,64 +391,28 @@ "name": "stderr", "output_type": "stream", "text": [ - "Operator `Kernel` run in 0.01 s\n" + "Operator `Kernel` ran in 0.01 s\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "#define _POSIX_C_SOURCE 200809L\n", - "#include \"stdlib.h\"\n", - "#include \"math.h\"\n", - "#include \"sys/time.h\"\n", - "#include \"xmmintrin.h\"\n", - "#include \"pmmintrin.h\"\n", - "\n", - "struct dataobj\n", + "/* Begin section0 */\n", + "START_TIMER(section0)\n", + "for (int x = x_m; x <= x_M; x += 1)\n", "{\n", - " void *restrict data;\n", - " int * size;\n", - " int * npsize;\n", - " int * dsize;\n", - " int * hsize;\n", - " int * hofs;\n", - " int * oofs;\n", - "} ;\n", - "\n", - "struct profiler\n", - "{\n", - " double section0;\n", - "} ;\n", - "\n", - "\n", - "int Kernel(struct dataobj *restrict f_vec, struct dataobj *restrict g_vec, struct profiler * timers, const int x_M, const int x_m, const int y_M, const int y_m)\n", - "{\n", - " float (*restrict f)[f_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[f_vec->size[1]]) f_vec->data;\n", - " float (*restrict g)[g_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[g_vec->size[1]]) g_vec->data;\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", - " struct timeval start_section0, end_section0;\n", - " gettimeofday(&start_section0, NULL);\n", - " /* Begin section0 */\n", - " for (int x = x_m; x <= x_M; x += 1)\n", + " #pragma omp simd aligned(f,g:32)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " #pragma omp simd aligned(f,g:32)\n", - " for (int y = y_m; y <= y_M; y += 1)\n", + " if (y < 5 && g[x + 1][y + 1] != 0)\n", " {\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", + " f[x + 1][y + 1] = f[x + 1][y + 1] + g[x + 1][y + 1];\n", " }\n", " }\n", - " /* End section0 */\n", - " gettimeofday(&end_section0, NULL);\n", - " timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;\n", - " return 0;\n", "}\n", - "\n" + "STOP_TIMER(section0,timers)\n", + "/* End section0 */\n" ] }, { @@ -523,7 +452,7 @@ "op3 = Operator(Eq(f, f + g, implicit_dims=ci))\n", "op3.apply()\n", "\n", - "print(op3.ccode)\n", + "print(op3.body.body[-1])\n", "\n", "assert (np.count_nonzero(f.data - np.diag(np.diagonal(f.data)))==0)\n", "assert (np.count_nonzero(f.data) == 10)\n", @@ -556,63 +485,27 @@ "name": "stderr", "output_type": "stream", "text": [ - "Operator `Kernel` run in 0.01 s\n" + "Operator `Kernel` ran in 0.01 s\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "#define _POSIX_C_SOURCE 200809L\n", - "#include \"stdlib.h\"\n", - "#include \"math.h\"\n", - "#include \"sys/time.h\"\n", - "#include \"xmmintrin.h\"\n", - "#include \"pmmintrin.h\"\n", - "\n", - "struct dataobj\n", - "{\n", - " void *restrict data;\n", - " int * size;\n", - " int * npsize;\n", - " int * dsize;\n", - " int * hsize;\n", - " int * hofs;\n", - " int * oofs;\n", - "} ;\n", - "\n", - "struct profiler\n", + "/* Begin section0 */\n", + "START_TIMER(section0)\n", + "for (int x = x_m; x <= x_M; x += 1)\n", "{\n", - " double section0;\n", - "} ;\n", - "\n", - "\n", - "int Kernel(struct dataobj *restrict g_vec, struct dataobj *restrict h_vec, struct profiler * timers, const int x_M, const int x_m, const int y_M, const int y_m)\n", - "{\n", - " float (*restrict g)[g_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[g_vec->size[1]]) g_vec->data;\n", - " float (*restrict h)[h_vec->size[1]] __attribute__ ((aligned (64))) = (float (*)[h_vec->size[1]]) h_vec->data;\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", - " struct timeval start_section0, end_section0;\n", - " gettimeofday(&start_section0, NULL);\n", - " /* Begin section0 */\n", - " for (int x = x_m; x <= x_M; x += 1)\n", + " for (int y = y_m; y <= y_M; y += 1)\n", " {\n", - " for (int y = y_m; y <= y_M; y += 1)\n", + " if (y < 5 && g[x + 1][y + 1] != 0)\n", " {\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", + " h[x + 1][y + 1] = g[x + 1][y + 1] + h[x + 1][y + 1];\n", " }\n", " }\n", - " /* End section0 */\n", - " gettimeofday(&end_section0, NULL);\n", - " timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;\n", - " return 0;\n", "}\n", - "\n" + "STOP_TIMER(section0,timers)\n", + "/* End section0 */\n" ] }, { @@ -643,7 +536,7 @@ "op4 = Operator(Eq(h, h + g))\n", "op4.apply()\n", "\n", - "print(op4.ccode)\n", + "print(op4.body.body[-1])\n", "\n", "assert (np.count_nonzero(h.data) == 5)\n", "\n", @@ -654,7 +547,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Use case B (Function subsampling)\n", + "# Use case B - Function subsampling\n", "\n", "`ConditionalDimension`s are additionally indicated to implement `Function` subsampling. In the following example, an `Operator` subsamples `Function` `g` and saves its content into `f` every `factor=4` iterations." ] @@ -668,60 +561,24 @@ "name": "stderr", "output_type": "stream", "text": [ - "Operator `Kernel` run in 0.01 s\n" + "Operator `Kernel` ran in 0.01 s\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ - "#define _POSIX_C_SOURCE 200809L\n", - "#include \"stdlib.h\"\n", - "#include \"math.h\"\n", - "#include \"sys/time.h\"\n", - "#include \"xmmintrin.h\"\n", - "#include \"pmmintrin.h\"\n", - "\n", - "struct dataobj\n", + "/* Begin section0 */\n", + "START_TIMER(section0)\n", + "for (int i = i_m; i <= i_M; i += 1)\n", "{\n", - " void *restrict data;\n", - " int * size;\n", - " int * npsize;\n", - " int * dsize;\n", - " int * hsize;\n", - " int * hofs;\n", - " int * oofs;\n", - "} ;\n", - "\n", - "struct profiler\n", - "{\n", - " double section0;\n", - "} ;\n", - "\n", - "\n", - "int Kernel(struct dataobj *restrict f_vec, struct dataobj *restrict g_vec, const int i_M, const int i_m, struct profiler * timers)\n", - "{\n", - " float (*restrict f) __attribute__ ((aligned (64))) = (float (*)) f_vec->data;\n", - " float (*restrict g) __attribute__ ((aligned (64))) = (float (*)) g_vec->data;\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", - " struct timeval start_section0, end_section0;\n", - " gettimeofday(&start_section0, NULL);\n", - " /* Begin section0 */\n", - " for (int i = i_m; i <= i_M; i += 1)\n", + " if ((i)%(4) == 0)\n", " {\n", - " if ((i)%(4) == 0)\n", - " {\n", - " f[i / 4] = g[i];\n", - " }\n", + " f[i / 4] = g[i];\n", " }\n", - " /* End section0 */\n", - " gettimeofday(&end_section0, NULL);\n", - " timers->section0 += (double)(end_section0.tv_sec-start_section0.tv_sec)+(double)(end_section0.tv_usec-start_section0.tv_usec)/1000000;\n", - " return 0;\n", "}\n", - "\n", + "STOP_TIMER(section0,timers)\n", + "/* End section0 */\n", "\n", " Data in g \n", " [ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15.]\n", @@ -743,7 +600,7 @@ "f = Function(name='f', shape=(int(size/factor),), dimensions=(ci,))\n", "\n", "op5 = Operator([Eq(f, g)])\n", - "print(op5.ccode)\n", + "print(op5.body.body[-1])\n", "\n", "op5.apply()\n", "\n", @@ -767,6 +624,198 @@ "`ConditionalDimension` can be used as an alternative to save snaps to disk as illustrated in [08_snapshotting](https://github.com/devitocodes/devito/blob/master/examples/seismic/tutorials/08_snapshotting.ipynb \"08_snapshotting\"). The example subsamples the `TimeDimension`, capturing snaps of the our wavefield at equally spaced snaps. " ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Example A - Count nonzero elements\n", + "\n", + "Here we present an example of using `ConditionalDimension` by counting the nonzero elements of an array." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Data([[1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],\n", + " [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.]], dtype=float32)" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from devito.types import Scalar\n", + "from devito import Inc\n", + "\n", + "grid = Grid(shape=shape)\n", + "# Define function 𝑓. We will initialize f's data with ones on its diagonal.\n", + "f = Function(name='f', shape=shape, dimensions=grid.dimensions, space_order=0)\n", + "f.data[:] = np.eye(shape[0])\n", + "f.data[:]" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/* Begin section0 */\n", + "START_TIMER(section0)\n", + "for (int x = x_m; x <= x_M; x += 1)\n", + "{\n", + " for (int y = y_m; y <= y_M; y += 1)\n", + " {\n", + " if (f[x][y] != 0)\n", + " {\n", + " int r0 = 1;\n", + " g[1] += r0;\n", + " }\n", + " }\n", + "}\n", + "STOP_TIMER(section0,timers)\n", + "/* End section0 */\n" + ] + }, + { + "data": { + "text/plain": [ + "10" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#NBVAL_IGNORE_OUTPUT\n", + "# ConditionalDimension with Ne(f, 0) condition\n", + "ci = ConditionalDimension(name='ci', parent=f.dimensions[-1],\n", + " condition=Ne(f, 0))\n", + "\n", + "# g is our counter. g needs to be a Function to write in\n", + "g = Function(name='g', shape=(1,), dimensions=(ci,), dtype=np.int32)\n", + "\n", + "# Extra Scalar used to avoid reduction in gcc-5\n", + "eq0 = Inc(g[0], 1, implicit_dims=(f.dimensions + (ci,)))\n", + "\n", + "op0 = Operator([eq0])\n", + "op0.apply()\n", + "\n", + "assert (np.count_nonzero(f.data) == g.data[0])\n", + "\n", + "print(op0.body.body[-1])\n", + "g.data[0]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Example B - Find nonzero positions \n", + "\n", + "Here we present another example of using `ConditionalDimension`. An `Operator` saves the nonzero positions of an array. We know that we have `g.data[0]` nonzero elements from Example A." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "Operator `Kernel` ran in 0.01 s\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "int k = -1;\n", + "\n", + "/* Begin section0 */\n", + "START_TIMER(section0)\n", + "for (int x = x_m; x <= x_M; x += 1)\n", + "{\n", + " for (int y = y_m; y <= y_M; y += 1)\n", + " {\n", + " if (f[x][y] != 0)\n", + " {\n", + " k += 1;\n", + " g[k][0] = x;\n", + " g[k][1] = y;\n", + " }\n", + " }\n", + "}\n", + "STOP_TIMER(section0,timers)\n", + "/* End section0 */\n" + ] + } + ], + "source": [ + "#NBVAL_IGNORE_OUTPUT\n", + "count = g.data[0] # Number of nonzeros\n", + "\n", + "# Dimension used only to nest different size of Functions under the same dim\n", + "id_dim = Dimension(name='id_dim')\n", + "\n", + "# Conditional for nonzero element\n", + "ci = ConditionalDimension(name='ci', parent=f.dimensions[-1],\n", + " condition=Ne(f, 0))\n", + "g = Function(name='g', shape=(count, len(f.dimensions)),\n", + " dimensions=(ci, id_dim), dtype=np.int32, space_order=0)\n", + "\n", + "eqs = []\n", + "\n", + "# Extra Scalar used to avoid reduction in gcc-5\n", + "k = Scalar(name='k', dtype=np.int32)\n", + "eqi = Eq(k, -1)\n", + "eqs.append(eqi)\n", + "eqii = Inc(k, 1, implicit_dims=(f.dimensions + (ci,)))\n", + "eqs.append(eqii)\n", + "\n", + "for n, i in enumerate(f.dimensions):\n", + " eqs.append(Eq(g[k, n], f.dimensions[n], implicit_dims=(f.dimensions + (ci,))))\n", + "\n", + "op0 = Operator(eqs)\n", + "op0.apply()\n", + "print(op0.body.body[-1])\n", + "\n", + "g.data[:]\n", + "\n", + "ids_np = np.nonzero(f.data[:])\n", + "for i in range(len(shape)-1):\n", + " assert all(g.data[:, i] == ids_np[i])" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -794,9 +843,9 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.5" + "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 2 -} \ No newline at end of file +}