inducer-pycuda

PyCUDA provides Pythonic access to NVIDIA's CUDA parallel computation API. It is designed to make CUDA programming more convenient and safer through several key features:

• RAII-style Object Cleanup: Resource cleanup is tied to the lifetime of Python objects. PyCUDA manages dependencies to ensure, for example, that a context is not detached before all memory allocated within it has been freed.
• High-level Abstractions: Provides convenient wrappers like pycuda.driver.SourceModule for managing CUDA source code and pycuda.gpuarray.GPUArray for managing data on the GPU.
• Full Driver API Access: Offers complete access to the CUDA driver API and includes interoperability support for OpenGL.
• Automatic Error Checking: CUDA errors are automatically caught and translated into standard Python exceptions.
• Performance: The core layer is implemented in C++ to ensure minimal overhead.

Metaprogramming in PyCUDA involves writing Python code that generates CUDA source code at runtime. This approach is used to solve several common CUDA optimization challenges:

• Automated Tuning: Instead of relying on unreliable heuristics for parameters like threads per block or shared memory size, you can benchmark different configurations at runtime and select the fastest one for the current hardware.
• Dynamic Data Types: Generate specialized kernels for specific types (e.g., float vs double) exactly when they are needed, avoiding the need to precompile multiple versions.
• Problem Specialization: Generate code tailored to the specific problem size or constraints, which allows for optimizations that generic code cannot achieve.
• Constant Optimization: Compiling problem sizes or specific values into the code as constants rather than variables. This reduces register pressure and improves performance (e.g., making multiplications more efficient).
• Manual Loop Unrolling: Since nvcc may not always unroll loops as expected with #pragma unroll, you can use Python to dynamically unroll loops to the required size during code generation.

PyCUDA allows binding memory to texture and surface units for specialized access patterns.

Texture References

Use TextureReference() to bind linear memory or an Array to a texture unit.

• set_array(array): Binds an Array object. The texture reference keeps the array alive.
• set_address(devptr, bytes, allow_offset=False): Binds a chunk of linear memory.
• set_address_2d(devptr, descr, pitch): Binds a 2D chunk of global memory.

Surface References

Use SurfaceReference() (constructed via Module.get_surfref) to bind an Array for read/write access.

• set_array(array, flags=0): Binds the surface to an Array.

Creating Arrays for Textures

To convert NumPy arrays or GPU arrays into Array objects suitable for textures/surfaces, use:

• np_to_array(nparray, order, allowSurfaceBind=False)
• gpuarray_to_array(gpuparray, order, allowSurfaceBind=False)
• matrix_to_array(matrix, order)
• make_multichannel_2d_array(matrix, order)

Managed memory (available in CUDA 6.0+) creates a virtual memory space visible to both CPU and GPU. The OS migrates pages between them automatically.

Warning: Accessing managed memory on the host while a GPU kernel is executing is strictly forbidden and will cause a segmentation fault that terminates the Python interpreter immediately.

Allocation Functions:

• managed_empty(shape, dtype, order="C", mem_flags=0)
• managed_zeros(shape, dtype, order="C", mem_flags=0)
• managed_empty_like(array, mem_flags=0)
• managed_zeros_like(array, mem_flags=0)

Usage Pattern: Managed arrays can be used on the host, passed to a kernel, and accessed on the host again without explicit memcpy calls. Always call context.synchronize() before accessing the array on the host after a kernel launch.

Visibility Control: Use ManagedAllocation.attach(mem_flags, stream=None) to change visibility using mem_attach_flags (e.g., GLOBAL or HOST).

from pycuda.autoinit import context
import pycuda.driver as cuda
import numpy as np
from pycuda.compiler import SourceModule

# 1. Allocate managed memory
a = cuda.managed_empty(shape=10, dtype=np.float32, mem_flags=cuda.mem_attach_flags.GLOBAL)
a[:] = np.linspace(0, 9, len(a)) # Fill on host

# 2. Use in kernel
mod = SourceModule("""
__global__ void doublify(float *a) { a[threadIdx.x] *= 2; }
""")
doublify = mod.get_function("doublify")
doublify(a, grid=(1,1), block=(len(a),1,1))

# 3. CRITICAL: Synchronize before host access
context.synchronize()

# 4. Access on host
median = np.median(a)

For high-quality pseudorandom or quasirandom numbers, use the CURAND-based classes in pycuda.curandom.

Important Resource Note: These generators run on the GPU, and each thread uses its own generator. Creating these objects is resource-intensive. On older Tesla devices, you may be limited to ~256 active generators; Fermi devices allow up to ~1024. If you encounter errors creating these objects, reduce the number of active threads/generators.

Pseudorandom Generators

• XORWOWRandomNumberGenerator: Provides pseudorandom numbers with a period of at least $2^{190}$. (CUDA 3.2+)
• MRG32k3aRandomNumberGenerator: Provides pseudorandom numbers. (CUDA 4.1+)

Quasirandom Generators

Quasirandom sequences are designed to fill n-dimensional space more evenly than pseudorandom sequences, though they are more expensive to generate.

• Sobol32RandomNumberGenerator: Period of $2^{32}$. (CUDA 3.2+)
• ScrambledSobol32RandomNumberGenerator: Scrambled version of Sobol32. (CUDA 4.0+)
• Sobol64RandomNumberGenerator: Period of $2^{64}$. (CUDA 4.0+)
• ScrambledSobol64RandomNumberGenerator: Scrambled version of Sobol64. (CUDA 4.0+)

Frequent calls to pycuda.driver.mem_alloc (e.g., when using gpuarray.GPUArray) can be slow. DeviceMemoryPool mitigates this by holding onto memory blocks instead of returning them to the system, allowing for faster reuse of similarly-sized blocks.

Key Components:

• DeviceMemoryPool: Manages a pool of linear device memory.
• allocate(size): Returns a PooledDeviceAllocation object.
• PooledDeviceAllocation: Represents the allocated memory. It can be cast to an int to get the starting device address. Calling .free() explicitly returns the memory to the pool.

Management Methods:

• free_held(): Frees all unused memory currently held by the pool.
• stop_holding(): Instructs the pool to immediately free memory returned to it instead of holding it. This is useful for cleanup.

Warning: Allocations made outside the pool may encounter out-of-memory errors if the pool has already claimed most of the available device memory.