sklam · May 23, 2022 · Mar 11, 2022 · Mar 16, 2022 · Apr 1, 2022 · Apr 5, 2022
diff --git a/docs/source/cuda/examples.rst b/docs/source/cuda/examples.rst
diff --git a/docs/source/cuda/laplace_final.svg b/docs/source/cuda/laplace_final.svg
diff --git a/docs/source/cuda/laplace_initial.svg b/docs/source/cuda/laplace_initial.svg
@@ -0,0 +1,76 @@
+import unittest
+
+from numba.cuda.testing import CUDATestCase, skip_on_cudasim
+from numba.tests.support import captured_stdout
+import numpy as np
+
+
+@skip_on_cudasim("cudasim doesn't support cuda import at non-top-level")
+class TestCpuGpuCompat(CUDATestCase):
+    """
+    Test compatibility of CPU and GPU functions
+    """
+
+    def setUp(self):
+        # Prevent output from this test showing up when running the test suite
+        self._captured_stdout = captured_stdout()
+        self._captured_stdout.__enter__()
+        super().setUp()
+
+    def tearDown(self):
+        # No exception type, value, or traceback
+        self._captured_stdout.__exit__(None, None, None)
+        super().tearDown()
+
+    def test_ex_cpu_gpu_compat(self):
+        # ex_cpu_gpu_compat.import.begin
+        from math import pi
+
+        import numba
+        from numba import cuda
+        # ex_cpu_gpu_compat.import.end
+
+        # ex_cpu_gpu_compat.allocate.begin
+        X = cuda.to_device([1, 10, 234])
+        Y = cuda.to_device([2, 2, 4014])
+        Z = cuda.to_device([3, 14, 2211])
+        results = cuda.to_device([0.0, 0.0, 0.0])
+        # ex_cpu_gpu_compat.allocate.end
+
+        # ex_cpu_gpu_compat.define.begin
+        @numba.jit
+        def business_logic(x, y, z):
+            return 4 * z * (2 * x - (4 * y) / 2 * pi)
+        # ex_cpu_gpu_compat.define.end
+
+        # ex_cpu_gpu_compat.cpurun.begin
+        print(business_logic(1, 2, 3))  # -126.79644737231007
+        # ex_cpu_gpu_compat.cpurun.end
+
+        # ex_cpu_gpu_compat.usegpu.begin
+        @cuda.jit
+        def f(res, xarr, yarr, zarr):
+            tid = cuda.grid(1)
+            if tid < len(xarr):
+                # The function decorated with numba.jit may be directly reused
+                res[tid] = business_logic(xarr[tid], yarr[tid], zarr[tid])
+        # ex_cpu_gpu_compat.usegpu.end
+
+        # ex_cpu_gpu_compat.launch.begin
+        f.forall(len(X))(results, X, Y, Z)
+        print(results)
+        # [-126.79644737231007, 416.28324559588634, -218912930.2987788]
+        # ex_cpu_gpu_compat.launch.end
+
+        expect = [
+            business_logic(x, y, z) for x, y, z in zip(X, Y, Z)
+        ]
+
+        np.testing.assert_equal(
+            expect,
+            results.copy_to_host()
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()
@@ -0,0 +1,152 @@
+import unittest
+
+from numba.cuda.testing import (CUDATestCase, skip_if_cudadevrt_missing,
+                                skip_on_cudasim)
+from numba.tests.support import captured_stdout
+
+
+@skip_if_cudadevrt_missing
+@skip_on_cudasim("cudasim doesn't support cuda import at non-top-level")
+class TestLaplace(CUDATestCase):
+    """
+    Test simple vector addition
+    """
+
+    def setUp(self):
+        # Prevent output from this test showing up when running the test suite
+        self._captured_stdout = captured_stdout()
+        self._captured_stdout.__enter__()
+        super().setUp()
+
+    def tearDown(self):
+        # No exception type, value, or traceback
+        self._captured_stdout.__exit__(None, None, None)
+        super().tearDown()
+
+    def test_ex_laplace(self):
+
+        # set True to regenerate the figures that
+        # accompany this example
+        plot = False
+
+        # ex_laplace.import.begin
+        import numpy as np
+        from numba import cuda
+        # ex_laplace.import.end
+
+        # ex_laplace.allocate.begin
+        # Use an odd problem size.
+        # This is so there can be an element truly in the "middle" for symmetry.
+        size = 1001
+        data = np.zeros(size)
+
+        # Middle element is made very hot
+        data[500] = 10000
+        buf_0 = cuda.to_device(data)
+
+        # This extra array is used for synchronization purposes
+        buf_1 = cuda.device_array_like(buf_0)
+
+        niter = 10000
+        # ex_laplace.allocate.end
+
+        if plot:
+            import matplotlib.pyplot as plt
+            fig, ax = plt.subplots(figsize=(16 * 0.66, 9 * 0.66))
+            plt.plot(
+                np.arange(len(buf_0)),
+                buf_0.copy_to_host(),
+                lw=3,
+                marker="*",
+                color='black'
+            )
+
+            plt.title('Initial State', fontsize=24)
+            plt.xlabel('Position', fontsize=24)
+            plt.ylabel('Temperature', fontsize=24)
+
+            ax.set_xticks(ax.get_xticks(), fontsize=16)
+            ax.set_yticks(ax.get_yticks(), fontsize=16)
+            plt.xlim(0, len(data))
+            plt.ylim(0, 10001)
+            plt.savefig('laplace_initial.svg')
+
+        # ex_laplace.kernel.begin
+        @cuda.jit
+        def solve_heat_equation(buf_0, buf_1, timesteps, k):
+            i = cuda.grid(1)
+
+            # Don't continue if our index is outside the domain
+            if i >= len(buf_0):
+                return
+
+            # Prepare to do a grid-wide synchronization later
+            grid = cuda.cg.this_grid()
+
+            for step in range(timesteps):
+                # Select the buffer from the previous timestep
+                if (step % 2) == 0:
+                    data = buf_0
+                    next_data = buf_1
+                else:
+                    data = buf_1
+                    next_data = buf_0
+
+                # Get the current temperature associated with this point
+                curr_temp = data[i]
+
+                # Apply formula from finite difference equation
+                if i == 0:
+                    # Left wall is held at T = 0
+                    next_temp = curr_temp + k * (data[i + 1] - (2 * curr_temp))
+                elif i == len(data) - 1:
+                    # Right wall is held at T = 0
+                    next_temp = curr_temp + k * (data[i - 1] - (2 * curr_temp))
+                else:
+                    # Interior points are a weighted average of their neighbors
+                    next_temp = curr_temp + k * (
+                        data[i - 1] - (2 * curr_temp) + data[i + 1]
+                    )
+
+                # Write new value to the next buffer
+                next_data[i] = next_temp
+
+                # Wait for every thread to write before moving on
+                grid.sync()
+        # ex_laplace.kernel.end
+
+        # ex_laplace.launch.begin
+        solve_heat_equation.forall(len(data))(
+            buf_0, buf_1, niter, 0.25
+        )
+        # ex_laplace.launch.end
+
+        results = buf_1.copy_to_host()
+        if plot:
+            fig, ax = plt.subplots(figsize=(16 * 0.66, 9 * 0.66))
+            plt.plot(
+                np.arange(len(results)),
+                results, lw=3,
+                marker="*",
+                color='black'
+            )
+            plt.title(f"T = {niter}", fontsize=24)
+            plt.xlabel('Position', fontsize=24)
+            plt.ylabel('Temperature', fontsize=24)
+
+            ax.set_xticks(ax.get_xticks(), fontsize=16)
+            ax.set_yticks(ax.get_yticks(), fontsize=16)
+
+            plt.ylim(0, max(results))
+            plt.xlim(0, len(results))
+            plt.savefig('laplace_final.svg')
+
+        # Integral over the domain should be equal to its initial value.
+        # Note that this should match the initial value of data[500] above, but
+        # we don't assign it to a variable because that would make the example
+        # code look a bit oddly verbose.
+        np.testing.assert_allclose(results.sum(), 10000)
+
+
+if __name__ == "__main__":
+    unittest.main()
@@ -0,0 +1,109 @@
+import unittest
+
+from numba.cuda.testing import CUDATestCase, skip_on_cudasim
+from numba.tests.support import captured_stdout
+
+
+@skip_on_cudasim("cudasim doesn't support cuda import at non-top-level")
+class TestMonteCarlo(CUDATestCase):
+    """
+    Test monte-carlo integration
+    """
+
+    def setUp(self):
+        # Prevent output from this test showing up when running the test suite
+        self._captured_stdout = captured_stdout()
+        self._captured_stdout.__enter__()
+        super().setUp()
+
+    def tearDown(self):
+        # No exception type, value, or traceback
+        self._captured_stdout.__exit__(None, None, None)
+        super().tearDown()
+
+    def test_ex_montecarlo(self):
+        # ex_montecarlo.import.begin
+        import numba
+        import numpy as np
+        from numba import cuda
+        from numba.cuda.random import (
+            create_xoroshiro128p_states,
+            xoroshiro128p_uniform_float32,
+        )
+        # ex_montecarlo.import.end
+
+        # ex_montecarlo.define.begin
+        # number of samples, higher will lead to a more accurate answer
+        nsamps = 1000000
+        # ex_montecarlo.define.end
+
+        # ex_montecarlo.kernel.begin
+        @cuda.jit
+        def mc_integrator_kernel(out, rng_states, lower_lim, upper_lim):
+            """
+            kernel to draw random samples and evaluate the function to
+            be integrated at those sample values
+            """
+            size = len(out)
+
+            gid = cuda.grid(1)
+            if gid < size:
+                # draw a sample between 0 and 1 on this thread
+                samp = xoroshiro128p_uniform_float32(rng_states, gid)
+
+                # normalize this sample to the limit range
+                samp = samp * (upper_lim - lower_lim) + lower_lim
+
+                # evaluate the function to be
+                # integrated at the normalized
+                # value of the sample
+                y = func(samp)
+                out[gid] = y
+        # ex_montecarlo.kernel.end
+
+        # ex_montecarlo.callfunc.begin
+        @cuda.reduce
+        def sum_reduce(a, b):
+            return a + b
+
+        def mc_integrate(lower_lim, upper_lim, nsamps):
+            """
+            approximate the definite integral of `func` from
+            `lower_lim` to `upper_lim`
+            """
+            out = cuda.to_device(np.zeros(nsamps, dtype="float32"))
+            rng_states = create_xoroshiro128p_states(nsamps, seed=42)
+
+            # jit the function for use in CUDA kernels
+
+            mc_integrator_kernel.forall(nsamps)(
+                out, rng_states, lower_lim, upper_lim
+            )
+            # normalization factor to convert
+            # to the average: (b - a)/(N - 1)
+            factor = (upper_lim - lower_lim) / (nsamps - 1)
+
+            return sum_reduce(out) * factor
+        # ex_montecarlo.callfunc.end
+
+        # ex_montecarlo.launch.begin
+        # define a function to integrate
+        @numba.jit
+        def func(x):
+            return 1.0 / x
+
+        mc_integrate(1, 2, nsamps)  # array(0.6929643, dtype=float32)
+        mc_integrate(2, 3, nsamps)  # array(0.4054021, dtype=float32)
+        # ex_montecarlo.launch.end
+
+        # values computed independently using maple
+        np.testing.assert_allclose(
+            mc_integrate(1, 2, nsamps), 0.69315, atol=0.001
+        )
+        np.testing.assert_allclose(
+            mc_integrate(2, 3, nsamps), 0.4055, atol=0.001
+        )
+
+
+if __name__ == "__main__":
+    unittest.main()