Examples added (#197)

* Added Cython examples

1. Cython/sycl_direct_linkage

Example of native extension "cdef import"-ing sycl C++ classes directly
from CL/sycl.hpp

Queue is created within the function, adding measurable overhead.

Extension uses GEMV to compute column-wise total of a C-contiguous matrix,
and illustrates linking to oneMKL.

2. Cython/sycl_bufer

Example of native extension building on the above, but illustrating getting
the queue from dpctl.

3. Cython/usm_memory

Example of native extension allocating USM shared memory via dpctl, and using
it as a buffer underlying NumPy array.

Cython functions dispatches to a SYCL code that works with USM pointer.

One function populates USM memory underneath NumPy array with random numbers
using ``oneapi::mkl::rng::device`` function used in SYCL kernel, with random
number being parameters of European vanilla options.

The second function uses SYCL to price these options using Black-Scholes
formula.

* updated create_sycl_queues.py to run on current dpctl

* Extended sycl_buffer example to implement column-wise summation without MKL

* few examples illustrating MemoryUSM* objects

* have it black's way

* Adjusted examples section of global README.md

Author: oleksandr-pavlyk

README.md

@@ -76,11 +76,17 @@ Examples
 ========
 See examples in folder `examples`.
 
-Run examples:
+Run python examples:
 ```bash
-python examples/create_sycl_queues.py
+for script in `ls examples/python/`; do echo "executing ${script}"; python examples/python/${script}; done
 ```
 
+Examples of building Cython extensions with DPC++ compiler, that interoperate with dpCtl can be found in 
+folder `cython`.
+
+Each example in `cython` folder can be built using `CC=clang CXX=dpcpp python setup.py build_ext --inplace`.
+Please refer to `run.py` script in respective folders to execute extensions.
+
 Tests
 =====
 See tests in folder `dpctl/tests`.

examples/cython/sycl_buffer/README.md

@@ -0,0 +1,80 @@
+#1 Example of SYCL extension working NumPy array input via SYCL buffers
+
+
+#2 Decription
+
+Cython function expecting a 2D array in C-contiguous layout that
+computes column-wise total by using SYCL oneMKL (as GEMV call with
+an all units vector).
+
+Example illustrates compiling SYCL extension, linking to oneMKL.
+
+
+#2 Compiling
+
+```
+# make sure oneAPI is activated, $ONEAPI_ROOT must be set
+CC=clang CXX=dpcpp python setup.py build_ext --inplace
+```
+
+
+#2 Running
+
+```
+# SYCL_BE=PI_OPENCL sets SYCL backend to OpenCL to avoid a
+# transient issue with MKL's using the default Level-0 backend
+(idp) [08:16:12 ansatnuc04 simple]$ SYCL_BE=PI_OPENCL ipython
+Python 3.7.7 (default, Jul 14 2020, 22:02:37)
+Type 'copyright', 'credits' or 'license' for more information
+IPython 7.17.0 -- An enhanced Interactive Python. Type '?' for help.
+
+In [1]: import syclbuffer as sb, numpy as np, dpctl
+
+In [2]: x = np.random.randn(10**4, 2500)
+
+In [3]: %time m1 = np.sum(x, axis=0)
+CPU times: user 22.3 ms, sys: 160 µs, total: 22.5 ms
+Wall time: 21.2 ms
+
+In [4]: %time m = sb.columnwise_total(x)  # first time is slower, due to JIT overhead
+CPU times: user 207 ms, sys: 36.1 ms, total: 243 ms
+Wall time: 248 ms
+
+In [5]: %time m = sb.columnwise_total(x)
+CPU times: user 8.89 ms, sys: 4.12 ms, total: 13 ms
+Wall time: 12.4 ms
+
+In [6]: %time m = sb.columnwise_total(x)
+CPU times: user 4.82 ms, sys: 8.06 ms, total: 12.9 ms
+Wall time: 12.3 ms
+```
+
+Running bench.py:
+
+```
+========== Executing warm-up ==========
+NumPy result:  [1. 1. 1. ... 1. 1. 1.]
+SYCL(Intel(R) Core(TM) i7-10710U CPU @ 1.10GHz) result: [1. 1. 1. ... 1. 1. 1.]
+SYCL(Intel(R) Gen9 HD Graphics NEO) result: [1. 1. 1. ... 1. 1. 1.]
+Times for 'opencl:cpu:0'
+[2.864787499012891, 2.690436460019555, 2.5902308400254697, 2.5802528870408423, 2.538990616973024]
+Times for 'opencl:gpu:0'
+[1.9769684099592268, 2.3491444009705447, 2.293720397981815, 2.391633405990433, 1.9465659779962152]
+Times for NumPy
+[3.4011058019823395, 3.07286038500024, 3.0390414349967614, 3.0305576199898496, 3.002687797998078]
+```
+
+Running run.py:
+
+```
+(idp) [09:14:53 ansatnuc04 sycl_buffer]$ SYCL_BE=PI_OPENCL python run.py
+Result computed by NumPy
+[  0.27170187 -23.36798583   7.31326489  -1.95121928]
+Result computed by SYCL extension
+[  0.27170187 -23.36798583   7.31326489  -1.95121928]
+
+Running on:  Intel(R) Gen9 HD Graphics NEO
+[  0.27170187 -23.36798583   7.31326489  -1.95121928]
+Running on:  Intel(R) Core(TM) i7-10710U CPU @ 1.10GHz
+[  0.27170187 -23.36798583   7.31326489  -1.95121928]
+```

examples/cython/sycl_buffer/_buffer_example.pyx

@@ -0,0 +1,28 @@
+cimport numpy as cnp
+import numpy as np
+
+cimport dpctl as c_dpctl
+import dpctl
+
+cdef extern from "use_sycl_buffer.h":
+    int c_columnwise_total(c_dpctl.DPPLSyclQueueRef q, size_t n, size_t m, double *m, double *ct) nogil
+    int c_columnwise_total_no_mkl(c_dpctl.DPPLSyclQueueRef q, size_t n, size_t m, double *m, double *ct) nogil
+
+def columnwise_total(double[:, ::1] v, method='mkl'):
+    cdef cnp.ndarray res_array = np.empty((v.shape[1],), dtype='d')
+    cdef double[::1] res_memslice = res_array
+    cdef int ret_status
+    cdef c_dpctl.SyclQueue q
+    cdef c_dpctl.DPPLSyclQueueRef q_ref
+
+    q = c_dpctl.get_current_queue()
+    q_ref = q.get_queue_ref()
+
+    if method == 'mkl':
+        with nogil:
+            ret_status = c_columnwise_total(q_ref, v.shape[0], v.shape[1], &v[0,0], &res_memslice[0])
+    else:
+        with nogil:
+            ret_status = c_columnwise_total_no_mkl(q_ref, v.shape[0], v.shape[1], &v[0,0], &res_memslice[0])        
+
+    return res_array

examples/cython/sycl_buffer/bench.py

@@ -0,0 +1,51 @@
+import dpctl
+import syclbuffer as sb
+import numpy as np
+
+X = np.full((10 ** 4, 4098), 1e-4, dtype="d")
+
+# warm-up
+print("=" * 10 + " Executing warm-up " + "=" * 10)
+print("NumPy result: ", X.sum(axis=0))
+
+dpctl.set_default_queue("opencl", "cpu", 0)
+print(
+    "SYCL({}) result: {}".format(
+        dpctl.get_current_queue().get_sycl_device().get_device_name(),
+        sb.columnwise_total(X),
+    )
+)
+
+dpctl.set_default_queue("opencl", "gpu", 0)
+print(
+    "SYCL({}) result: {}".format(
+        dpctl.get_current_queue().get_sycl_device().get_device_name(),
+        sb.columnwise_total(X),
+    )
+)
+
+import timeit
+
+print("Times for 'opencl:cpu:0'")
+print(
+    timeit.repeat(
+        stmt="sb.columnwise_total(X)",
+        setup='dpctl.set_default_queue("opencl", "cpu", 0); '
+        "sb.columnwise_total(X)",  # ensure JIT compilation is not counted
+        number=100,
+        globals=globals(),
+    )
+)
+
+print("Times for 'opencl:gpu:0'")
+print(
+    timeit.repeat(
+        stmt="sb.columnwise_total(X)",
+        setup='dpctl.set_default_queue("opencl", "gpu", 0); sb.columnwise_total(X)',
+        number=100,
+        globals=globals(),
+    )
+)
+
+print("Times for NumPy")
+print(timeit.repeat(stmt="X.sum(axis=0)", number=100, globals=globals()))

examples/cython/sycl_buffer/run.py

@@ -0,0 +1,22 @@
+import syclbuffer as sb
+import numpy as np
+
+X = np.random.randn(100, 4)
+
+print("Result computed by NumPy")
+print(X.sum(axis=0))
+print("Result computed by SYCL extension")
+print(sb.columnwise_total(X))
+
+
+print("")
+# controlling where to offload
+import dpctl
+
+with dpctl.device_context("opencl:gpu"):
+    print("Running on: ", dpctl.get_current_queue().get_sycl_device().get_device_name())
+    print(sb.columnwise_total(X))
+
+with dpctl.device_context("opencl:cpu"):
+    print("Running on: ", dpctl.get_current_queue().get_sycl_device().get_device_name())
+    print(sb.columnwise_total(X))

examples/cython/sycl_buffer/setup.py

@@ -0,0 +1,67 @@
+import sys
+from os.path import join, exists, abspath, dirname
+from os import getcwd
+from os import environ
+from Cython.Build import cythonize
+
+
+def configuration(parent_package="", top_path=None):
+    from numpy.distutils.misc_util import Configuration
+    from numpy.distutils.system_info import get_info
+    import numpy as np
+    import dpctl
+
+    config = Configuration("", parent_package, top_path)
+
+    oneapi_root = environ.get("ONEAPI_ROOT", None)
+    if not oneapi_root:
+        raise ValueError("ONEAPI_ROOT must be set, typical value is /opt/intel/oneapi")
+
+    mkl_info = {
+        "include_dirs": [join(oneapi_root, "mkl", "include")],
+        "library_dirs": [
+            join(oneapi_root, "mkl", "lib"),
+            join(oneapi_root, "mkl", "lib", "intel64"),
+        ],
+        "libraries": [
+            "mkl_sycl",
+            "mkl_intel_ilp64",
+            "mkl_tbb_thread",
+            "mkl_core",
+            "tbb",
+            "iomp5",
+        ],
+    }
+
+    mkl_include_dirs = mkl_info.get("include_dirs")
+    mkl_library_dirs = mkl_info.get("library_dirs")
+    mkl_libraries = mkl_info.get("libraries")
+
+    pdir = dirname(__file__)
+    wdir = join(pdir)
+
+    eca = ["-Wall", "-Wextra", "-fsycl", "-fsycl-unnamed-lambda"]
+
+    config.add_extension(
+        name="syclbuffer",
+        sources=[
+            join(pdir, "_buffer_example.pyx"),
+            join(wdir, "use_sycl_buffer.cpp"),
+            join(wdir, "use_sycl_buffer.h"),
+        ],
+        include_dirs=[wdir, np.get_include(), dpctl.get_include()] + mkl_include_dirs,
+        libraries=["sycl"] + mkl_libraries,
+        runtime_library_dirs=mkl_library_dirs,
+        extra_compile_args=eca,  # + ['-O0', '-g', '-ggdb'],
+        extra_link_args=["-fPIC"],
+        language="c++",
+    )
+
+    config.ext_modules = cythonize(config.ext_modules, include_path=[pdir, wdir])
+    return config
+
+
+if __name__ == "__main__":
+    from numpy.distutils.core import setup
+
+    setup(configuration=configuration)

examples/cython/sycl_buffer/use_sycl_buffer.cpp

@@ -0,0 +1,109 @@
+#include <CL/sycl.hpp>
+#include "use_sycl_buffer.h"
+#include <oneapi/mkl.hpp>
+#include "dppl_sycl_types.h"
+
+int
+c_columnwise_total(DPPLSyclQueueRef q_ref, size_t n, size_t m, double *mat, double *ct) {
+
+    sycl::queue q = *(reinterpret_cast<sycl::queue *>(q_ref));
+
+    sycl::buffer<double, 1> mat_buffer = sycl::buffer(mat, sycl::range<1>(n * m));
+    sycl::buffer<double, 1> ct_buffer = sycl::buffer(ct, sycl::range<1>(m));
+
+    double *ones = reinterpret_cast<double *>(malloc(n * sizeof(double)));
+    {
+	sycl::buffer<double, 1> ones_buffer = sycl::buffer(ones, sycl::range<1>(n));
+
+	try {
+	    auto ev = q.submit([&](sycl::handler &cgh) {
+				   auto ones_acc = ones_buffer.get_access<sycl::access::mode::read_write>(cgh);
+				   cgh.fill(ones_acc, double(1.0));
+			       });
+	    
+	    ev.wait_and_throw();
+	}
+	catch (sycl::exception const& e) {
+	    std::cout << "\t\tCaught synchronous SYCL exception during fill:\n"
+		      << e.what() << std::endl << "OpenCL status: " << e.get_cl_code() << std::endl;
+	    goto cleanup;
+	}
+
+	try {
+	    oneapi::mkl::blas::row_major::gemv(
+		q,
+		oneapi::mkl::transpose::trans,
+		n, m, double(1.0), mat_buffer, m,
+		ones_buffer, 1,
+		double(0.0), ct_buffer, 1);
+	    q.wait();
+	}
+	catch (sycl::exception const &e) {
+	    std::cout << "\t\tCaught synchronous SYCL exception during GEMV:\n"
+		      << e.what() << std::endl << "OpenCL status: " << e.get_cl_code() << std::endl;
+	    goto cleanup;
+	}
+    }
+    
+    free(ones);
+    return 0;
+
+  cleanup:
+    free(ones);
+    return -1;
+}
+
+inline size_t upper_multiple(size_t n, size_t wg) { return wg * ((n + wg - 1)/wg); }
+
+int
+c_columnwise_total_no_mkl(DPPLSyclQueueRef q_ref, size_t n, size_t m, double *mat, double *ct) {
+
+    sycl::queue q = *(reinterpret_cast<sycl::queue *>(q_ref));
+
+    sycl::buffer<double, 2> mat_buffer = sycl::buffer(mat, sycl::range<2>(n, m));
+    sycl::buffer<double, 1> ct_buffer = sycl::buffer(ct, sycl::range<1>(m));
+
+    auto e = q.submit(
+	[&](sycl::handler &h) {
+	    sycl::accessor ct_acc {ct_buffer, h, sycl::write_only};
+	    h.parallel_for(
+		sycl::range<1>(m),
+		[=](sycl::id<1> i){
+		    ct_acc[i] = 0.0;
+		});
+	});
+
+    constexpr size_t wg = 256;
+    auto e2 = q.submit(
+	[&](sycl::handler &h) {
+
+	    sycl::accessor mat_acc {mat_buffer, h, sycl::read_only};
+	    sycl::accessor ct_acc {ct_buffer, h};
+	    h.depends_on(e);
+
+	    sycl::range<2> global {upper_multiple(n, wg), m};
+	    sycl::range<2> local {wg, 1};
+
+	    h.parallel_for(
+		sycl::nd_range<2>(global, local),
+		[=](sycl::nd_item<2> it) {
+		    size_t i = it.get_global_id(0);
+		    size_t j = it.get_global_id(1);
+		    double group_sum = sycl::ONEAPI::reduce(
+			it.get_group(),
+			(i < n) ? mat_acc[it.get_global_id()] : 0.0, 
+			std::plus<double>()
+			);
+		    if (it.get_local_id(0) == 0) {
+			sycl::ONEAPI::atomic_ref<
+			    double,
+			    sycl::ONEAPI::memory_order::relaxed,
+			    sycl::ONEAPI::memory_scope::system,
+			    sycl::access::address_space::global_space>(ct_acc[j]) += group_sum;
+		    }
+		});
+	});
+
+    e2.wait_and_throw();
+    return 0;
+}

examples/cython/sycl_buffer/use_sycl_buffer.h

@@ -0,0 +1,7 @@
+#include <CL/sycl.hpp>
+#include "dppl_sycl_types.h"
+
+extern int c_columnwise_total(
+    DPPLSyclQueueRef q, size_t n, size_t m, double *mat, double *ct);
+extern int c_columnwise_total_no_mkl(
+    DPPLSyclQueueRef q, size_t n, size_t m, double *mat, double *ct);