Merge pull request #242 from IntelPython/feature/add_numba_dpex_n_implementations

Feature/add numba dpex n implementations

Author: adarshyoga

dpbench/benchmarks/npbench/azimint_hist/azimint_hist_numba_dpex_n.py

@@ -0,0 +1,61 @@
+# SPDX-FileCopyrightText: 2014 Jérôme Kieffer et al.
+# SPDX-FileCopyrightText: 2021 ETH Zurich and the NPBench authors
+# SPDX-FileCopyrightText: 2022 - 2023 Intel Corporation
+#
+# SPDX-License-Identifier: BSD-3-Clause
+
+"""
+Jérôme Kieffer and Giannis Ashiotis. Pyfai: a python library for
+high performance azimuthal integration on gpu, 2014. In Proceedings of the
+7th European Conference on Python in Science (EuroSciPy 2014).
+"""
+
+import dpnp as np
+from numba_dpex import dpjit
+
+
+@dpjit
+def get_bin_edges_parallel(a, bins):
+    bin_edges = np.zeros((bins + 1,), dtype=np.float64)
+    a_min = a.min()
+    a_max = a.max()
+    delta = (a_max - a_min) / bins
+    for i in range(bin_edges.shape[0]):
+        bin_edges[i] = a_min + i * delta
+
+    bin_edges[-1] = a_max  # Avoid roundoff error on last point
+    return bin_edges
+
+
+@dpjit
+def compute_bin(x, bin_edges):
+    # assuming uniform bins for now
+    n = bin_edges.shape[0] - 1
+    a_min = bin_edges[0]
+    a_max = bin_edges[-1]
+
+    # special case to mirror NumPy behavior for last bin
+    if x == a_max:
+        return n - 1  # a_max always in last bin
+
+    return int(n * (x - a_min) / (a_max - a_min))
+
+
+@dpjit
+def histogram_parallel(a, bins, weights):
+    hist = np.zeros((bins,), dtype=a.dtype)
+    bin_edges = get_bin_edges_parallel(a, bins)
+
+    for i in range(a.shape[0]):
+        bin = compute_bin(a[i], bin_edges)
+        hist[bin] += weights[i]
+
+    return hist, bin_edges
+
+
+@dpjit
+def azimint_hist(data, radius, npt):
+    histu = np.histogram(radius, npt)[0]
+    # histw = np.histogram(radius, npt, weights=data)[0]
+    histw = histogram_parallel(radius, npt, weights=data)[0]
+    return histw / histu

dpbench/benchmarks/npbench/azimint_naive/azimint_naive_numba_dpex_n.py

@@ -0,0 +1,27 @@
+# SPDX-FileCopyrightText: 2014 Jérôme Kieffer et al.
+# SPDX-FileCopyrightText: 2021 ETH Zurich and the NPBench authors
+# SPDX-FileCopyrightText: 2022 - 2023 Intel Corporation
+#
+# SPDX-License-Identifier: BSD-3-Clause
+
+"""
+Jérôme Kieffer and Giannis Ashiotis. Pyfai: a python library for
+high performance azimuthal integration on gpu, 2014. In Proceedings of the
+7th European Conference on Python in Science (EuroSciPy 2014).
+"""
+
+import dpnp as np
+from numba_dpex import dpjit
+
+
+@dpjit
+def azimint_naive(data, radius, npt):
+    rmax = radius.max()
+    res = np.zeros(npt, dtype=np.float64)
+    for i in range(npt):
+        r1 = rmax * i / npt
+        r2 = rmax * (i + 1) / npt
+        mask_r12 = np.logical_and((r1 <= radius), (radius < r2))
+        values_r12 = data[mask_r12]
+        res[i] = values_r12.mean()
+    return res

dpbench/benchmarks/npbench/channel_flow/channel_flow_numba_dpex_n.py

@@ -0,0 +1,228 @@
+# SPDX-FileCopyrightText: 2017 Lorena A. Barba, Gilbert F. Forsyth.
+# SPDX-FileCopyrightText: 2018 Barba, Lorena A., and Forsyth, Gilbert F.
+# SPDX-FileCopyrightText: 2021 ETH Zurich and the NPBench authors
+# SPDX-FileCopyrightText: 2022 - 2023 Intel Corporation
+#
+# SPDX-License-Identifier: BSD-3-Clause
+
+"""
+CFD Python: the 12 steps to Navier-Stokes equations.
+Journal of Open Source Education, 1(9), 21,
+https://doi.org/10.21105/jose.00021
+"""
+
+import dpnp as np
+from numba_dpex import dpjit
+
+
+@dpjit
+def build_up_b(rho, dt, dx, dy, u, v):
+    b = np.zeros_like(u)
+    b[1:-1, 1:-1] = rho * (
+        1
+        / dt
+        * (
+            (u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx)
+            + (v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)
+        )
+        - ((u[1:-1, 2:] - u[1:-1, 0:-2]) / (2 * dx)) ** 2
+        - 2
+        * (
+            (u[2:, 1:-1] - u[0:-2, 1:-1])
+            / (2 * dy)
+            * (v[1:-1, 2:] - v[1:-1, 0:-2])
+            / (2 * dx)
+        )
+        - ((v[2:, 1:-1] - v[0:-2, 1:-1]) / (2 * dy)) ** 2
+    )
+
+    # Periodic BC Pressure @ x = 2
+    b[1:-1, -1] = rho * (
+        1
+        / dt
+        * (
+            (u[1:-1, 0] - u[1:-1, -2]) / (2 * dx)
+            + (v[2:, -1] - v[0:-2, -1]) / (2 * dy)
+        )
+        - ((u[1:-1, 0] - u[1:-1, -2]) / (2 * dx)) ** 2
+        - 2
+        * (
+            (u[2:, -1] - u[0:-2, -1])
+            / (2 * dy)
+            * (v[1:-1, 0] - v[1:-1, -2])
+            / (2 * dx)
+        )
+        - ((v[2:, -1] - v[0:-2, -1]) / (2 * dy)) ** 2
+    )
+
+    # Periodic BC Pressure @ x = 0
+    b[1:-1, 0] = rho * (
+        1
+        / dt
+        * (
+            (u[1:-1, 1] - u[1:-1, -1]) / (2 * dx)
+            + (v[2:, 0] - v[0:-2, 0]) / (2 * dy)
+        )
+        - ((u[1:-1, 1] - u[1:-1, -1]) / (2 * dx)) ** 2
+        - 2
+        * (
+            (u[2:, 0] - u[0:-2, 0])
+            / (2 * dy)
+            * (v[1:-1, 1] - v[1:-1, -1])
+            / (2 * dx)
+        )
+        - ((v[2:, 0] - v[0:-2, 0]) / (2 * dy)) ** 2
+    )
+
+    return b
+
+
+@dpjit
+def pressure_poisson_periodic(nit, p, dx, dy, b):
+    pn = np.empty_like(p)
+
+    for q in range(nit):
+        pn = p.copy()
+        p[1:-1, 1:-1] = (
+            (pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dy**2
+            + (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dx**2
+        ) / (2 * (dx**2 + dy**2)) - dx**2 * dy**2 / (
+            2 * (dx**2 + dy**2)
+        ) * b[
+            1:-1, 1:-1
+        ]
+
+        # Periodic BC Pressure @ x = 2
+        p[1:-1, -1] = (
+            (pn[1:-1, 0] + pn[1:-1, -2]) * dy**2
+            + (pn[2:, -1] + pn[0:-2, -1]) * dx**2
+        ) / (2 * (dx**2 + dy**2)) - dx**2 * dy**2 / (
+            2 * (dx**2 + dy**2)
+        ) * b[
+            1:-1, -1
+        ]
+
+        # Periodic BC Pressure @ x = 0
+        p[1:-1, 0] = (
+            (pn[1:-1, 1] + pn[1:-1, -1]) * dy**2
+            + (pn[2:, 0] + pn[0:-2, 0]) * dx**2
+        ) / (2 * (dx**2 + dy**2)) - dx**2 * dy**2 / (
+            2 * (dx**2 + dy**2)
+        ) * b[
+            1:-1, 0
+        ]
+
+        # Wall boundary conditions, pressure
+        p[-1, :] = p[-2, :]  # dp/dy = 0 at y = 2
+        p[0, :] = p[1, :]  # dp/dy = 0 at y = 0
+
+
+@dpjit
+def channel_flow(nit, u, v, dt, dx, dy, p, rho, nu, F):
+    udiff = 1
+    stepcount = 0
+
+    while udiff > 0.001:
+        un = u.copy()
+        vn = v.copy()
+
+        b = build_up_b(rho, dt, dx, dy, u, v)
+        pressure_poisson_periodic(nit, p, dx, dy, b)
+
+        u[1:-1, 1:-1] = (
+            un[1:-1, 1:-1]
+            - un[1:-1, 1:-1] * dt / dx * (un[1:-1, 1:-1] - un[1:-1, 0:-2])
+            - vn[1:-1, 1:-1] * dt / dy * (un[1:-1, 1:-1] - un[0:-2, 1:-1])
+            - dt / (2 * rho * dx) * (p[1:-1, 2:] - p[1:-1, 0:-2])
+            + nu
+            * (
+                dt
+                / dx**2
+                * (un[1:-1, 2:] - 2 * un[1:-1, 1:-1] + un[1:-1, 0:-2])
+                + dt
+                / dy**2
+                * (un[2:, 1:-1] - 2 * un[1:-1, 1:-1] + un[0:-2, 1:-1])
+            )
+            + F * dt
+        )
+
+        v[1:-1, 1:-1] = (
+            vn[1:-1, 1:-1]
+            - un[1:-1, 1:-1] * dt / dx * (vn[1:-1, 1:-1] - vn[1:-1, 0:-2])
+            - vn[1:-1, 1:-1] * dt / dy * (vn[1:-1, 1:-1] - vn[0:-2, 1:-1])
+            - dt / (2 * rho * dy) * (p[2:, 1:-1] - p[0:-2, 1:-1])
+            + nu
+            * (
+                dt
+                / dx**2
+                * (vn[1:-1, 2:] - 2 * vn[1:-1, 1:-1] + vn[1:-1, 0:-2])
+                + dt
+                / dy**2
+                * (vn[2:, 1:-1] - 2 * vn[1:-1, 1:-1] + vn[0:-2, 1:-1])
+            )
+        )
+
+        # Periodic BC u @ x = 2
+        u[1:-1, -1] = (
+            un[1:-1, -1]
+            - un[1:-1, -1] * dt / dx * (un[1:-1, -1] - un[1:-1, -2])
+            - vn[1:-1, -1] * dt / dy * (un[1:-1, -1] - un[0:-2, -1])
+            - dt / (2 * rho * dx) * (p[1:-1, 0] - p[1:-1, -2])
+            + nu
+            * (
+                dt / dx**2 * (un[1:-1, 0] - 2 * un[1:-1, -1] + un[1:-1, -2])
+                + dt / dy**2 * (un[2:, -1] - 2 * un[1:-1, -1] + un[0:-2, -1])
+            )
+            + F * dt
+        )
+
+        # Periodic BC u @ x = 0
+        u[1:-1, 0] = (
+            un[1:-1, 0]
+            - un[1:-1, 0] * dt / dx * (un[1:-1, 0] - un[1:-1, -1])
+            - vn[1:-1, 0] * dt / dy * (un[1:-1, 0] - un[0:-2, 0])
+            - dt / (2 * rho * dx) * (p[1:-1, 1] - p[1:-1, -1])
+            + nu
+            * (
+                dt / dx**2 * (un[1:-1, 1] - 2 * un[1:-1, 0] + un[1:-1, -1])
+                + dt / dy**2 * (un[2:, 0] - 2 * un[1:-1, 0] + un[0:-2, 0])
+            )
+            + F * dt
+        )
+
+        # Periodic BC v @ x = 2
+        v[1:-1, -1] = (
+            vn[1:-1, -1]
+            - un[1:-1, -1] * dt / dx * (vn[1:-1, -1] - vn[1:-1, -2])
+            - vn[1:-1, -1] * dt / dy * (vn[1:-1, -1] - vn[0:-2, -1])
+            - dt / (2 * rho * dy) * (p[2:, -1] - p[0:-2, -1])
+            + nu
+            * (
+                dt / dx**2 * (vn[1:-1, 0] - 2 * vn[1:-1, -1] + vn[1:-1, -2])
+                + dt / dy**2 * (vn[2:, -1] - 2 * vn[1:-1, -1] + vn[0:-2, -1])
+            )
+        )
+
+        # Periodic BC v @ x = 0
+        v[1:-1, 0] = (
+            vn[1:-1, 0]
+            - un[1:-1, 0] * dt / dx * (vn[1:-1, 0] - vn[1:-1, -1])
+            - vn[1:-1, 0] * dt / dy * (vn[1:-1, 0] - vn[0:-2, 0])
+            - dt / (2 * rho * dy) * (p[2:, 0] - p[0:-2, 0])
+            + nu
+            * (
+                dt / dx**2 * (vn[1:-1, 1] - 2 * vn[1:-1, 0] + vn[1:-1, -1])
+                + dt / dy**2 * (vn[2:, 0] - 2 * vn[1:-1, 0] + vn[0:-2, 0])
+            )
+        )
+
+        # Wall BC: u,v = 0 @ y = 0,2
+        u[0, :] = 0
+        u[-1, :] = 0
+        v[0, :] = 0
+        v[-1, :] = 0
+
+        udiff = (np.sum(u) - np.sum(un)) / np.sum(u)
+        stepcount += 1
+
+    return stepcount

dpbench/benchmarks/npbench/compute/compute_numba_dpex_n.py

@@ -0,0 +1,17 @@
+# SPDX-FileCopyrightText: 2023 Stefan Behnel, Robert Bradshaw,
+#   Dag Sverre Seljebotn, Greg Ewing, William Stein, Gabriel Gellner, et al.
+# SPDX-FileCopyrightText: 2021 ETH Zurich and the NPBench authors
+# SPDX-FileCopyrightText: 2022 - 2023 Intel Corporation
+#
+# SPDX-License-Identifier: BSD-3-Clause
+
+# https://cython.readthedocs.io/en/latest/src/userguide/numpy_tutorial.html
+
+import dpnp as np
+from numba_dpex import dpjit
+
+
+@dpjit
+def compute(array_1, array_2, a, b, c):
+    # return np.clip(array_1, 2, 10) * a + array_2 * b + c
+    return np.minimum(np.maximum(array_1, 2), 10) * a + array_2 * b + c

dpbench/benchmarks/npbench/crc16/crc16_numba_dpex_n.py

@@ -0,0 +1,27 @@
+# SPDX-FileCopyrightText: 2018 Øystein Sture
+#
+# SPDX-License-Identifier: BSD-3-Clause
+
+import dpnp as np
+from numba_dpex import dpjit
+
+
+# Adapted from https://gist.github.com/oysstu/68072c44c02879a2abf94ef350d1c7c6
+@dpjit
+def crc16(data, poly=0x8408):
+    """
+    CRC-16-CCITT Algorithm
+    """
+    crc = 0xFFFF
+    for b in data:
+        cur_byte = 0xFF & b
+        for _ in range(0, 8):
+            if (crc & 0x0001) ^ (cur_byte & 0x0001):
+                crc = (crc >> 1) ^ poly
+            else:
+                crc >>= 1
+            cur_byte >>= 1
+    crc = ~crc & 0xFFFF
+    crc = (crc << 8) | ((crc >> 8) & 0xFF)
+
+    return crc & 0xFFFF