<jupyter_start><jupyter_text>PyCUDA: querying devices and simple kernels This notebook is based on the book *Hands on GPU programming with Python and CUDA* by <NAME> (ISBN 978-1-78899-391-3).  Querying device information Import the CUDA driver.<jupyter_code>import pycuda.driver as drv<jupyter_output><empty_output><jupyter_text>Initialize the driver.<jupyter_code>drv.init()<jupyter_output><empty_output><jupyter_text>Query the number of available devices.<jupyter_code>drv.Device.count()<jupyter_output><empty_output><jupyter_text>For convenience, assign the device to a variable.<jupyter_code>dev = drv.Device(0)<jupyter_output><empty_output><jupyter_text>Get the name of the device.<jupyter_code>dev.name()<jupyter_output><empty_output><jupyter_text>Verify the CUDA compute capability of this device.<jupyter_code>dev.compute_capability()<jupyter_output><empty_output><jupyter_text>Check the clock frequency of the GPU device.<jupyter_code>print(f'{dev.get_attribute(drv.device_attribute.CLOCK_RATE)/1000} MHz')<jupyter_output>1480.5 MHz<jupyter_text>Check the size of a warp.<jupyter_code>dev.get_attribute(drv.device_attribute.WARP_SIZE)<jupyter_output><empty_output><jupyter_text>Check the number of SMs on the GPU.<jupyter_code>dev.get_attribute(drv.device_attribute.MULTIPROCESSOR_COUNT)<jupyter_output><empty_output><jupyter_text>Check the maximum number of resident threads per SM.<jupyter_code>dev.get_attribute(drv.device_attribute.MAX_THREADS_PER_MULTIPROCESSOR)<jupyter_output><empty_output><jupyter_text>Blocks & grids Verify the maximum number of threads in a CUDA block.<jupyter_code>dev.get_attribute(drv.device_attribute.MAX_THREADS_PER_BLOCK)<jupyter_output><empty_output><jupyter_text>The maximum number of threads per dimension in a block varies.<jupyter_code>for max_dim in [drv.device_attribute.MAX_BLOCK_DIM_X, drv.device_attribute.MAX_BLOCK_DIM_Y, drv.device_attribute.MAX_BLOCK_DIM_Z]:
    print(f'{max_dim}: {dev.get_attribute(max_dim)}')<jupyter_output>MAX_BLOCK_DIM_X: 1024
MAX_BLOCK_DIM_Y: 1024
MAX_BLOCK_DIM_Z: 64<jupyter_text>Similarly, the maximum number of blocks per dimension in a grid varies.<jupyter_code>for max_dim in [drv.device_attribute.MAX_GRID_DIM_X, drv.device_attribute.MAX_GRID_DIM_Y, drv.device_attribute.MAX_GRID_DIM_Z]:
    print(f'{max_dim}: {dev.get_attribute(max_dim)}')<jupyter_output>MAX_GRID_DIM_X: 2147483647
MAX_GRID_DIM_Y: 65535
MAX_GRID_DIM_Z: 65535<jupyter_text>Memory Check the number of 4-byte registers per block.<jupyter_code>dev.get_attribute(drv.device_attribute.MAX_REGISTERS_PER_BLOCK)<jupyter_output><empty_output><jupyter_text>Check the number of 4-byte regeisters per SM.<jupyter_code>dev.get_attribute(drv.device_attribute.MAX_REGISTERS_PER_MULTIPROCESSOR)<jupyter_output><empty_output><jupyter_text>Check the maximum shared memory available for a block.<jupyter_code>dev.get_attribute(drv.device_attribute.MAX_SHARED_MEMORY_PER_BLOCK)<jupyter_output><empty_output><jupyter_text>Check the maximum shared memory available on an SM.<jupyter_code>dev.get_attribute(drv.device_attribute.MAX_SHARED_MEMORY_PER_MULTIPROCESSOR)<jupyter_output><empty_output><jupyter_text>Check the size of L2 cache.<jupyter_code>dev.get_attribute(drv.device_attribute.L2_CACHE_SIZE)<jupyter_output><empty_output><jupyter_text>Check the total memory of the device.<jupyter_code>dev.total_memory()/1024**3<jupyter_output><empty_output><jupyter_text>Simple computations Import and initialize the driver if not already done. Import `gpuarray` to get data into the device memory.<jupyter_code>import numpy as np
import pycuda.autoinit
from pycuda import gpuarray<jupyter_output><empty_output><jupyter_text>Create a numpy array in the host memory.<jupyter_code>host_data = np.array([1.1, 2.2, 3.3, 4.4], dtype=np.float32)<jupyter_output><empty_output><jupyter_text>Transfer the array to device memory and perform a simple computation on it.<jupyter_code>gpu_data = gpuarray.to_gpu(host_data)
gpu_data_new = 3.14*gpu_data<jupyter_output><empty_output><jupyter_text>Transfer the array from device memory to host memory.<jupyter_code>host_data_new = gpu_data_new.get()
host_data_new<jupyter_output><empty_output><jupyter_text>Let's measure performance. Create a host array of 10 million elements.<jupyter_code>host_data = np.random.uniform(size=10_000_000).astype(np.float32)<jupyter_output><empty_output><jupyter_text>Time a simple computation, including transferring data to and from device memory.<jupyter_code>%%timeit -n 1 -r 1
gpu_data = gpuarray.to_gpu(host_data)
gpu_data_new = gpu_data*np.float32(2.0)
host_data_new = gpu_data_new.get()<jupyter_output>25.6 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)<jupyter_text>Time the duration of the actual computation, so without data transfer overhead.<jupyter_code>gpu_data = gpuarray.to_gpu(host_data)
%%timeit -n 1 -r 1
gpu_data_new = gpu_data*np.float32(2.0)
%timeit np.float32(2.0)*host_data<jupyter_output>2.36 ms ± 323 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)<jupyter_text>Although it would in theory be possible to perform computation using this style of programming, this would be very inefficient.  Computations on a GPGPU are typically coded using kernel functions.  Elementwise kernels Required imports.<jupyter_code>import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import pycuda.autoinit
from pycuda import gpuarray
from pycuda.elementwise import ElementwiseKernel<jupyter_output><empty_output><jupyter_text>We will compute the Mandelbrot set.  First, initialize the initial values for $z$.<jupyter_code>min_val, max_val, nr_points = -1.8, 1.8, 2**14
lattice = (np.linspace(min_val, max_val, nr_points, dtype=np.complex64).reshape((1, -1)) +
           np.complex64(1.0j)*np.linspace(min_val, max_val, nr_points, dtype=np.complex64).reshape((-1, 1)))<jupyter_output><empty_output><jupyter_text>Transfer the data to the GPU.<jupyter_code>lattice_gpu = gpuarray.to_gpu(lattice)<jupyter_output><empty_output><jupyter_text>Create an empty array of the same shape on the GPU.<jupyter_code>graph_gpu = gpuarray.empty(lattice.shape, dtype=np.float32)<jupyter_output><empty_output><jupyter_text>Define the CUDA elementwise kernel.  * The first argument represents the parameters of the kernel function.  Note that 2D numpy arrays will be passed as 1D C arrays to the kernel function.  * The second argument is the source code for the kernel function.  Note that an index `i` for the arrays passed as arguments is automatically generated.  * The third argument is the name of the kernel function.  We can use the C types, exept for complex numbers, which should be typed as `pycuda::complex` for single precision floating point numbers.`int` on the device is 32-bit, so the Python 3 integer should be converted into a numpy 32-bit integer.<jupyter_code>kernel = ElementwiseKernel(
    'pycuda::complex<float> *lattice, float *graph, int max_iters, float upper_bound',
    '''
    graph[i] = 1.0f;
    pycuda::complex<float> c = lattice[i];
    pycuda::complex<float> z(0.0f, 0.0f);
    for (int j = 0; j < max_iters; ++j) {
        z = z*z + c;
        if (abs(z) > upper_bound) {
            graph[i] = 0.0f;
            break;
        }
    }
    ''',
    'madelbrot_kernel'
)<jupyter_output><empty_output><jupyter_text>Launch the kernel function, retrieve and display the result.<jupyter_code>kernel(lattice_gpu, graph_gpu, np.int32(255), np.float32(2.0))
graph = graph_gpu.get()
plt.imshow(graph);<jupyter_output><empty_output><jupyter_text>Scan kernels Required imports.<jupyter_code>import numpy as np
import pycuda.autoinit
from pycuda import gpuarray
from pycuda.scan  import InclusiveScanKernel<jupyter_output><empty_output><jupyter_text>We perform a scan on an array, computing the running sum of its elements.<jupyter_code>host_data = np.arange(1, 30, 3, dtype=np.int32)
host_data
gpu_data = gpuarray.to_gpu(host_data)<jupyter_output><empty_output><jupyter_text>Define and run the kernel.<jupyter_code>kernel = InclusiveScanKernel(np.int32, "a + b")
kernel(gpu_data)<jupyter_output><empty_output><jupyter_text>Note that the computation has been done in place, the result is just returned forconvenience.<jupyter_code>gpu_data.get()<jupyter_output><empty_output><jupyter_text>Reduction kernels Required imports.<jupyter_code>import numpy as np
import pycuda.autoinit
from pycuda import gpuarray
from pycuda.reduction import ReductionKernel
host_data = np.random.normal(size=100_000_000).astype(np.float32)
gpu_data = gpuarray.to_gpu(host_data)<jupyter_output><empty_output><jupyter_text>We define a kernel that computes the maximum of the array.  This can be done by an elementwise reduction with the first array element as the neutral element.<jupyter_code>max_kernel = ReductionKernel(np.float32,
                             neutral=str(host_data[0]),
                             reduce_expr='a > b ? a : b',
                             arguments='float *in')
max_kernel(gpu_data)<jupyter_output><empty_output><jupyter_text>The result is identical when computed on the CPU using numpy.<jupyter_code>host_data.max()<jupyter_output><empty_output><jupyter_text>A similar kernel can be defined for the minimum of an array.<jupyter_code>min_kernel = ReductionKernel(np.float32,
                             neutral='0.0',
                             reduce_expr='a > b ? b : a',
                             arguments='float *in')
min_kernel(gpu_data)
host_data.min()<jupyter_output><empty_output><jupyter_text>Note that the result as returned as a numpy array with a single element. The `ReductionKernel` is in fact more powerful, it can be used to implement a map-reduce operation.  We can implement the scalar product of two vectors easily.<jupyter_code>dot_kernel = ReductionKernel(np.float32,
                             neutral='0.0',
                             reduce_expr='a + b',
                             map_expr='v1[i]*v2[i]',
                             arguments='float *v1, float *v2')<jupyter_output><empty_output><jupyter_text>Create two vectors on the host system, and transfer to the device memory.<jupyter_code>host_v1 = np.random.uniform(size=10_000_000).astype(np.float32)
host_v2 = np.random.uniform(size=10_000_000).astype(np.float32)
gpu_v1 = gpuarray.to_gpu(host_v1)
gpu_v2 = gpuarray.to_gpu(host_v2)
dot_kernel(gpu_v1, gpu_v2)<jupyter_output><empty_output><jupyter_text>The result is nearly the same on the CPU, except for the least signicant digit due to rounding differences consistent with single precision arithmetic.<jupyter_code>host_v1@host_v2<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0935224586288416, 'fraction_numerical': 0.0216548463356974, 'mean_word_length': 6.0225763612217795, 'pattern_counts': {'":': 0, '<': 173, '<?xml version=': 0, '>': 175, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '43', 'dir': 'jupyter-structured-clean-dedup', 'id': '647116', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 3774, 'n_tokens_neox': 3591, 'n_words': 828, 'path': 'Python/PyCuda/simple_pycuda.ipynb', 'repo_name': 'Gjacquenot/training-material', 'size': '10583'}