<jupyter_start><jupyter_text>GPU programming with PyOpenCLLayout:1. Comparison of modern CPU and GPU2. Introduction to OpenCL 3. Comparison with other ways of programming GPU4. My very first kernel: Element-wise addition5. Parallel programming design patterns6. Metaprogramming in PyOpenCL7. Conclusions  Comparison of modern CPU and GPU    They are not that different: both are many-core processors.* CPUs are still faster (MHz) and have lower latency to memory* GPUs are have more compute power and larger bandwidth to memory  Introduction to OpenCL* Vendor & platform neutral parallel programming language to address CPU, GPU, FPGA and other types of accelerators. * Initated by Apple in 2008 and managed by the Khronos group, currently at version 3.0. * There is a strict separation between host code and device code (called kernels).* Kernels are written in a subset of C99 (without pointers, nor recursion) and compiled at _runtime_.* Bindings for host code, initially in C are available for all programming languages.* Sycl (Intel) is the descendent of OpenCL where C++ host code can be mixed with kernels (à la CUDA).* There is a direct mapping between CUDA kernels and OpenCL (version 1.2).  Definition:In parallel programming, a `kernel` is the core part of the computational code, without the outer loop. It runs in parallel on all cores of the device.<jupyter_code>import numpy as np

size = 1000
a = np.random.random(size)
b = np.random.random(size)
res = np.empty_like(a)

for idx in range(size):
    res[idx] = a[idx] + b[idx]

# ⚠ This is not vectorized numpy code ⚠<jupyter_output><empty_output><jupyter_text>The kernel is `res[i] = a[i] + b[i]`.The size of the problem, `size`, is precised at runtime, usually in the host code.The position of the working element `idx` is obtained via specific API calls.  GPU programming languages ...      PyOpenCL* Expose the full OpenCL API* Explicit kernels in C99 * Explicit memory management* Explicit device and context management* Asynchronous execution with queues and events* Runtime compilation (since the driver contains the compiler)* Many metaprogramming features to expose advanced parallel algorithms* Integration in Jupyter notebooks<jupyter_code>#Integration with Jupyter notebooks:
%load_ext pyopencl.ipython_ext

import pyopencl as cl
import pyopencl.array as cla
#Creation of the OpenCL context which holds the memory on the device
ctx = cl.create_some_context(interactive=True)

# The command queue is used to interact with the device asynchronously
queue = cl.CommandQueue(ctx, 
                         properties=cl.command_queue_properties.PROFILING_ENABLE)<jupyter_output>Choose platform:
[0] <pyopencl.Platform 'NVIDIA CUDA' at 0x5653b7405a40>
[1] <pyopencl.Platform 'Intel(R) OpenCL' at 0x5653b73fff00>
Choice [0]:0
Choose device(s):
[0] <pyopencl.Device 'GeForce GTX TITAN' on 'NVIDIA CUDA' at 0x5653b73dc750>
[1] <pyopencl.Device 'Quadro M2000' on 'NVIDIA CUDA' at 0x5653b73dc170>
Choice, comma-separated [0]:0
Set the environment variable PYOPENCL_CTX='0:0' to avoid being asked again.<jupyter_text>Preparation of some test data:* GPUs need to work on large chunks of data to be efficient, here we work on 64 MB of data* Allocation of the memory on the device and transfer of the data* Convention: the `_d` suffix indicates the object lives on the device<jupyter_code># Array creation
shape = (2**24,)
print(f"size of the probem {shape}")
a = np.random.random(shape).astype(np.float32)
b = np.random.random(shape).astype(np.float32)
print(f"Allocated memory for one buffer: {a.nbytes}")

#Reference result, calculated in float32
ref = a+b 

# Send data to the devicem using the `queue`
a_d = cla.to_device(queue, a)
b_d = cla.to_device(queue, b)
#prepare output buffer
res_d = cla.empty_like(a_d)

print(type(res_d), res_d)<jupyter_output>size of the probem (16777216,)
Allocated memory for one buffer: 67108864
<class 'pyopencl.array.Array'> [0. 0. 0. ... 0. 0. 0.]<jupyter_text>Structure of an OpenCL `kernel`:* Kernels are `void` functions: they do not return anything.* Kernels are declared with the `kernel` keyword* Arrays are declared `global` which indicates the memory of the device* The position of the worker is obtained with `get_global_id()` function<jupyter_code>%%cl_kernel
// OpenCL kernel written in Jupyter:

kernel void add(global float* a,
                global float* b,
                global float* res){
    int idx = get_global_id(0);
    res[idx] = a[idx] + b[idx];
}<jupyter_output><empty_output><jupyter_text>Integration with Jupyter:* The `%%cl_kernel` on the first line indicates the cell contains OpenCL code* Comments need to be indicated like in C* The code is compiled when the cell is executed* Kernels are exposed in the Python namespace automatically!<jupyter_code>%%cl_kernel
// Parallel programming hello-world

kernel void add(global float* a,
                global float* b,
                global float* res){
    int idx = get_global_id(0);
    res[idx] = a[idx] + b[idx];
}
# Execution of a kernel:

evt = add(queue, shape, None,             # Indicate the queue, the size of the problem and the chunking (None=auto)
          a_d.data, b_d.data, res_d.data) # This needs to match the signature of the kernel
print(evt)                                # Since execution is asynchronous, the event is the handle to the processing 
print("The results are correct: ",
np.allclose(ref, res_d.get()))
print(f"Execution time on GPU: {1e-6*(evt.profile.end-evt.profile.start):.3f} ms\nExecution time on CPU:")
%timeit a+b<jupyter_output>Execution time on GPU: 0.906 ms
Execution time on CPU:
18.7 ms ± 126 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)<jupyter_text>Parallel programming design patterns* Map: Element-wise operations like the `add` kernel.* Gather: Stencil like operation, for example convolutions.* Scatter: write at variable position, requires atomic operations.* Reduction: perform the sum for all elements in an array.* Scan: Perform the `cumsum`, sum of all previous elements, used in compactions.* Sort: Bitonic sort is one example of parallel sort.* ... Beside `Map`, all kernel require dozens to hundreeds of lines of code to implement the algorithm!PyOpenCL provides templates to all those algorithms making programmer's life simpler.  Generating `Map` kernels<jupyter_code>from pyopencl.elementwise import ElementwiseKernel

t_add = ElementwiseKernel(ctx, 
                          arguments="float* a, float* b, float* res", 
                          operation="res[i] = a[i] + b[i]")

res_d.fill(0.0)                  # Reset the destination array
t_add(a_d, b_d, res_d)           # Launch the kernel

print("Results are correct:",
np.allclose(ref, res_d.get()))
# Even more trivial:
c_d = a_d + b_d

%time a_d + b_d

type(c_d)<jupyter_output>CPU times: user 242 µs, sys: 100 µs, total: 342 µs
Wall time: 348 µs<jupyter_text>Reduction kernel, like the sum of all elements in an array:![Reduction kernel](kernel-code-sum-reduction.png)One can use `reductions` to perform the scalar product ... `a.b = sum(a[i]*b[i])` where `a[i]*b[i]` is the **map** operation and the sum is the **reduction**.Reduction are not limited to addition, any commutative operation is usable like `*`, `min`, `max`, ... **⚠** addition is not strictly commutative for floating point numbers, but that's another story.<jupyter_code># Dot product implemented as reduction kernel
from pyopencl.reduction import ReductionKernel
dot = ReductionKernel(ctx, 
                      dtype_out=np.float32, 
                      neutral="0",
                      reduce_expr="a + b", 
                      map_expr="a[i] * b[i]",
                      arguments="__global float* a, __global float* b")

print("Validation of the result in 32-bits:",
np.isclose(dot(a_d, b_d).get(), np.dot(a,b)))

# 32bits precision is not enough ! 
a64 = a.astype(np.float64)
b64 = b.astype(np.float64)
print("Validation of the result in 64-bits:",
np.isclose(dot(a_d, b_d).get(), np.dot(a64,b64)))
print("Performance evaluation on the GPU:")
%timeit dot(a_d, b_d).get()

print("Performance evaluation on CPU (32 bits):")
%timeit np.dot(a,b)
print("Performance evaluation on CPU (64 bits):")
%timeit np.dot(a64,b64)<jupyter_output>Performance evaluation on the GPU:
847 µs ± 3.74 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
Performance evaluation on CPU (32 bits):
9.25 ms ± 16.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Performance evaluation on CPU (64 bits):
7.09 ms ± 685 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)<jupyter_text>Visualisation of the code generated by PyOpenCLReduction can have multiple stages, here is it implemented with 2 kernels working with workgroups of 256 threads...* Stage 2 is reducing 256 values coming from stage 1* Stage 1 is reducing 256² values, i.e. 64k threads are used *simulateneously** Stage 0 is a simple loop over elements, each thread is reading 1k elements in a `for` loop<jupyter_code>print(dot.stage_2_inf.source) 
print(dot.stage_1_inf.source)
source = dot.stage_1_inf.source + "\n" + dot.stage_2_inf.source
nblines  = len(source.split('\n'))
print(f"Length of the source code: {nblines} lines")<jupyter_output>Length of the source code: 313 lines<jupyter_text>Scan kernel, sum of all previous elements from the array:![Scan kernel](300px-Prefix_sum_16.svg.png)Their typical application is *in memory* compaction or in compression algorithms. Once again, this is not limited to the `+` operation, it is valid for any commutative operation.In `numpy` this is implemented in the `cumsum` function:<jupyter_code>from pyopencl.scan import GenericScanKernel
cumsum = GenericScanKernel(ctx, 
                           np.float32,
                           arguments="__global float* ary, __global float* out",
                           input_expr="ary[i]",
                           scan_expr="a+b", 
                           neutral="0",
                           output_statement="out[i] = item;")
cumsum(a_d, res_d)
print("Results are corrects:", np.allclose(res_d.get(),np.cumsum(a64)))
print("Performances of `cumsum` implemented in pyopencl:")
%timeit cumsum(a_d, res_d).wait()
print("Performances of `cumsum` implemented in numpy (32 bits):")
%timeit np.cumsum(a)
print("Performances of `cumsum` implemented in numpy (64 bits):")
%timeit np.cumsum(a64)<jupyter_output>Performances of `cumsum` implemented in pyopencl:
8.47 ms ± 27 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
Performances of `cumsum` implemented in numpy (32 bits):
44.2 ms ± 109 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
Performances of `cumsum` implemented in numpy (64 bits):
50.2 ms ± 340 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
{'fraction_non_alphanumeric': 0.09533996837503488, 'fraction_numerical': 0.025858059715375314, 'mean_word_length': 4.057384760112888, 'pattern_counts': {'":': 0, '<': 35, '<?xml version=': 0, '>': 35, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '80849', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3525, 'n_tokens_neox': 3240, 'n_words': 1315, 'path': 'opencl/OpenCL.ipynb', 'repo_name': 'silx-kit/silx-training', 'size': '10751'}