CUDA by Numba Examples Half 1 | by Carlos Costa | Medium

Observe this sequence to find out about CUDA programming from scratch with Python

GPUs (graphics processing items), because the title implies, have been initially developed for laptop graphics. Since then, they’ve turn out to be ubiquitous in nearly each space that requires excessive computational throughput. This progress has been enabled by the event of GPGPU (common objective GPU) interfaces, which permit us to program GPUs for general-purpose computing. The most typical of those interfaces is CUDA, adopted by OpenCL and most not too long ago, HIP.

CUDA was initially designed to be appropriate with C. Later variations prolonged it to C++ and Fortran. Within the Python ecosystem, one of many methods of utilizing CUDA is thru Numba, a Simply-In-Time (JIT) compiler for Python that may goal GPUs (it additionally targets CPUs, however that’s exterior of our scope). With Numba, one can write kernels instantly with (a subset of) Python, and Numba will compile the code on-the-fly and run it. Whereas it doesn’t implement the entire CUDA API, its supported options are sometimes sufficient to acquire spectacular speedups in comparison with CPUs (for all lacking options, see the Numba documentation).

Numba just isn’t the one choice, nonetheless. CuPy provides each excessive degree capabilities which depend on CUDA below the hood, low-level CUDA assist for integrating kernels written in C, and JIT-able Python capabilities (just like Numba). PyCUDA gives much more fine-grained management of the CUDA API. Extra not too long ago, Nvidia launched the official CUDA Python, which is able to certainly enrich the ecosystem. All of those tasks can go machine arrays to one another, you aren’t locked into utilizing just one.

The aim of this sequence is to supply a studying platform for frequent CUDA patterns by way of examples written in Numba CUDA. What this sequence just isn’t, is a complete information to both CUDA or Numba. The reader could check with their respective documentations for that. The construction of this tutorial is impressed by the ebook CUDA by Instance: An Introduction to Normal-Goal GPU Programming by Jason Sanders and Edward Kandrot. If you happen to ultimately develop out of Python and need to code in C, it is a wonderful useful resource.

We’ll discover ways to run our first Numba CUDA kernel. We may even discover ways to use CUDA effectively for embarrassingly parallel duties, that’s, duties that are utterly impartial from one another. Lastly, we’ll discover ways to time our kernel runtimes from the CPU.

Click on right here to seize the code in Google Colab.

The largest benefit GPUs have over CPUs is their skill to execute the identical directions in parallel. A single CPU core will run directions serially, one after the opposite. Parallelizing over a CPU requires using its a number of cores (bodily or digital) on the identical time. An ordinary trendy laptop has 4–8 cores. Then again, trendy GPUs have lots of if not hundreds of compute cores. See Determine 1 for a comparability between these two. GPU cores are usually slower and may solely execute easy directions, however their sheer quantity often makes up for these shortcomings manyfold. The caveat is that to ensure that GPUs to have an fringe of CPUs, the algorithms they run have to be parallelizable.

I consider there are 4 most important features to grokking GPU programming. The primary I already talked about: understanding methods to suppose and design algorithms which are parallel by nature. This may be exhausting each as a result of some algorithms are designed serially, but in addition as a result of there might be some ways of parallelizing the identical algorithm.

The second side is studying methods to map buildings which sit on the host akin to vectors and pictures, onto GPU constructs akin to threads and blocks. Recurring patterns and helper capabilities can help us on this, however on the finish of the day, experimentation shall be essential to get probably the most out of your GPU.

The third is comprehending the asynchronous execution mannequin that drives GPU programming. Not solely GPUs and CPUs execute directions independently from one another, GPUs have streams which permit a number of processing streams to run in the identical GPU. This asynchronicity is essential when designing optimum processing flows.

The fourth and ultimate side is the relation between summary ideas and concrete code: that is achieved by studying the API and its nuances.

As you learn this primary chapter, attempt to establish these ideas within the following examples!

We’ll begin by establishing the environment: a Numba model increased than 0.55 and a supported GPU.

The principle workhorse of Numba CUDA is the cuda.jit decorator. It’s used to outline capabilities which is able to run within the GPU.

We’ll begin by defining a easy operate, which takes two numbers and shops them on the primary component of the third argument. Our first lesson is that kernels (GPU capabilities that launch threads) can’t return values. We get round that by passing inputs and outputs. This can be a frequent sample in C, however not quite common in Python.

As you could have seen, earlier than we name the kernel, we have to allocate an array on the machine. As well as, if we need to show the returned worth, we have to copy it again to the CPU. You could be asking your self why we selected to allocate a float32 (single-precision float). It is because, whereas supported in most trendy GPUs, double precision arithmetic can take 4x or longer than single precision arithmetic. So it is higher to get used to utilizing np.float32 and np.complex64 as a substitute of float/np.float64 and complicated/np.complex128.

Whereas the kernel definition appears just like a CPU operate, the kernel name is a bit completely different. Particularly, it has sq. brackets earlier than the arguments:

add_scalars[1, 1](2.0, 7.0, dev_c)

These sq. brackets check with the variety of blocks in a grid, and the variety of threads in a block, respectively. Let’s speak a bit bit extra about what these imply as we be taught to parallelize with CUDA.

The anatomy of a CUDA grid

When a kernel is launched it has a grid related to it. A grid consists of blocks; a block consists of threads. Determine 2 reveals a one dimensional CUDA grid. The grid within the determine has 4 blocks. The variety of blocks in a grid is held in a particular variable which might be accessed contained in the kernel known as gridDim.x. The .x is refers back to the first dimensional of the grid (the one one on this case). Two dimensional grids even have .y and three dimensional grids, .z variables. As of 2022, there aren’t any four-dimensional grids or increased. Additionally contained in the kernel, you’ll find out which block is being executed by way of using blockIdx.x, which on this case will run from 0 to three.

Every block has a sure variety of threads, held within the variable blockDim.x. Thread indices are held within the variable threadIdx.x, which on this instance will run from 0 to 7.

Importantly, threads in numerous blocks are scheduled to run in another way, have entry to completely different reminiscence areas, and differ in another methods (see CUDA Refresher: The CUDA Programming Mannequin for a quick dialogue). For now, we’ll skip these particulars.

After we launched the kernel in our first instance with parameters [1, 1], we informed CUDA to run one block with one thread. Passing a number of blocks with a number of threads, will run the kernel many occasions. Manipulating threadIdx.x and blockIdx.x will enable us to uniquely establish every thread.

As a substitute of summing two numbers, let’s attempt to sum two arrays. Suppose the arrays every have 20 components. Like within the determine above, we will launch a kernel with 8 threads per block. If we wish every thread to deal with just one array component, we’ll then want at the very least 4 blocks. Launching 4 blocks, 8 threads every, our grid will then launch 32 threads.

Now we have to determine methods to map the thread indices to the array indices. threadIdx.x runs from 0 to 7, so on their very own they aren’t in a position to index our array. As well as, completely different blocks have the identical threadIdx.x. Then again, they’ve completely different blockIdx.x. To acquire a singular index for every thread, we will mix these variables:

i = threadIdx.x + blockDim.x * blockIdx.x

For the primary block, blockIdx.x = 0 and i will run from 0 to 7. For the second block, blockIdx.x = 1. Since blockDim.x = 8, i will run from 8 to fifteen. Equally, for blockIdx.x = 2, i will run from 16 to 23. Within the fourth and ultimate block, i will run from 24 to 31. See Desk 1 beneath.

We solved one downside: methods to map every thread to every component within the array… however now we now have a difficulty the place some threads would overflow the array, because the array has 20 components and i goes as much as 32-1. The answer is straightforward: for these threads, do not do something!

Let’s see the code.

In newer variations of Numba, we get a warning noting that we known as the kernel with host arrays. Ideally, we need to keep away from transferring information round from host to machine, as that is very gradual. We needs to be calling the kernel with machine arrays in all arguments. We will try this by transferring the array from host to machine beforehand:

dev_a = cuda.to_device(a)dev_b = cuda.to_device(b)

Furthermore, the calculation of distinctive indices per thread can get previous rapidly. Fortunately Numba gives the quite simple wrapper cuda.grid which is known as with the grid dimension as the one argument. The brand new kernel will appear to be this:

What occurs once we change the dimensions of the array? One simple means out is to easily change the grid parameters (variety of blocks and threads per block) with a view to launch at the very least as many threads as there are components within the array.

There may be some science and a few artwork to setting these parameters. For the “science”, we’ll say that (a) they need to be a a number of of two, sometimes between 32 and 1024, and (b) they need to be chosen in order to maximise occupancy (what number of threads are lively on the identical time). Nvidia gives a spreadsheet that may assist calculating these. For the “artwork”, nothing can predict the conduct of your kernels, so if you happen to actually need to optimize these parameters, you could profile your code with typical inputs. In apply, a “cheap” variety of threads for contemporary GPUs is 256.

Earlier than we transfer on from summing vectors, we have to speak about {hardware} limits. GPUs can’t run an arbitrary variety of threads and blocks. Sometimes every block can’t have greater than 1024 threads, and a grid can’t have greater than 2¹⁶ − 1 = 65535 blocks. This isn’t to say which you could launch 1024 × 65535 threads… there are limits to the variety of threads that may be launched based mostly on how a lot reminiscence their registers occupy, amongst different concerns. Furthermore, one have to be cautious of attempting to course of giant arrays which don’t match within the GPU RAM all of sudden. In these instances, one could profit from processing the arrays piecewise, both utilizing a single GPU or a number of GPUs.

INFO: In Python, {hardware} limits might be obtained by way of Nvidia’s cuda-python library by way of the operate cuDeviceGetAttribute of their documentation. See the Appendix on the finish of this part for an instance.

Grid-stride loops

In instances the place the variety of blocks per grid exceeds the {hardware} restrict however the array suits in reminiscence, as a substitute of utilizing one thread per array component, we will use one thread to course of a number of components. We’ll accomplish that through the use of a way known as grid-stride loops. In addition to overcoming {hardware} limitations, grid-stride loop kernels profit from reusing threads, by minimizing thread creation/destruction overhead. Mark Harris’ weblog submit CUDA Professional Tip: Write Versatile Kernels with Grid-Stride Loops goes into element about among the advantages of grid-strided loops.

The concept behind this system is so as to add a loop within the CUDA kernel to course of a number of enter components. The stride of this loop, because the title implies, is the same as the variety of threads in a grid. This fashion, if the full variety of threads within the grid (threads_per_grid = blockDim.x * gridDim.x) is smaller than the variety of components of the array, as quickly because the kernel is finished processing the index cuda.grid(1) it’ll course of the index cuda.grid(1) + threads_per_grid and so forth till all array components have been processed. With out additional ado, let’s take a look at the code.

This code is similar to the above, with the completely different that we’re beginning at cuda.grid(1), however executing extra samples, one each threads_per_grid till we hit the tip of the array.

Now, which one in every of these kernels is quicker?

GPU programming is all about pace. Subsequently it is very important measure code execution precisely.

CUDA kernels are machine capabilities which are launched by the host (CPU), however in fact they’re executed on the GPU. The GPU and the CPU don’t talk except we inform them to. So when the GPU kernel is launched, the CPU will merely proceed operating directions, be they launching extra kernels or executing different CPU capabilities. If we place a time.time() name earlier than and after kernel launch, we shall be timing solely how lengthy it takes for the kernel to launch, to not run.

One operate we will use to make sure that the GPU has “caught up” is the cuda.synchronize(). Calling this operate will cease the host from executing every other code till the GPU finishes execution of each kernel that has been launched in it.

To time a kernel execution, we will then merely time how lengthy it takes for the kernel to run after which synchronize. There are two caveats to this. First, we have to use time.perf_counter() or time.perf_counter_ns()and never time.time(). time.time() doesn’t depend the time that the host is sleeping ready for the GPU to complete execution. The second caveat is that timing code from the host just isn’t ultimate as there are overheads associated to this. Later, we’ll clarify how one can use CUDA occasions to time kernels from the machine. Mark Harris has one other glorious weblog submit about this subject entitled Tips on how to Implement Efficiency Metrics in CUDA C/C++.

When utilizing Numba, we now have one element we should take note of. Numba is a Simply-In-Time compiler, that means that the capabilities are solely compiled when they’re known as. Subsequently timing the primary name of the operate may even time the compilation step which is normally a lot slower. We should bear in mind to all the time compile the code first by launching the kernel after which synchronizing it to make sure that nothing is left to run within the GPU. This ensures that the following kernel runs instantly with out compilation. Additionally observe that the dtype of the array needs to be the identical, as Numba compiles a singular operate for every mixture of argument dtypes.

For easy kernels, we will additionally measure the all through of the algorithm, which equals the variety of floating level operations per second. It’s often measured in GFLOP/s (giga-FLOP per second). Our including operation comprises just one FLOP: addition. As such, the throughput is given by:

To finish this tutorial, let’s craft a 2D kernel to use logarithmic correction to a picture.

Given a picture I(x, y) with values between 0 and 1, the log-corrected picture is given by

Iᵪ(x, y) = γ log₂ (1 + I(x, y))

First let’s seize some information!

As you’ll be able to see, the information is actually saturated on the decrease finish. There are nearly no values above 0.6.

Let’s write the kernel.

Allow us to make observe of the 2 for loops. Discover that that the primary for loop begins at iy and the second, innermost loop begins at ix. We might have simply chosen i0 to start out at ix and i1 to start out at iy as a substitute, which might even really feel extra pure. So why did we select this order? It seems that the reminiscence entry sample for the primary selection is extra environment friendly. Because the first grid index is the quickest one, we need to make it match our quickest dimension: the final one.

If you happen to don’t need to take my phrase for it (and also you shouldn’t!) you might have now discovered methods to time kernel executions, and you may check out the 2 variations. For small arrays such because the one utilizing right here, the distinction is negligible, however for bigger arrays (say 10,000 by 10,000), I’ve measured a speedup of about 10%. Not tremendous spectacular, but when I might offer you a ten% enchancment with a single swapping of variables, who wouldn’t take it?

And that’s it! We will now see extra particulars within the corrected picture.

As an train, strive timing completely different launches with completely different grids to seek out the optimum grid dimension in your machine.

On this tutorial you discovered the fundamentals of Numba CUDA. You discovered methods to create easy CUDA kernels, and transfer reminiscence to GPU to make use of them. You additionally discovered methods to iterate over 1D and 2D arrays utilizing a way known as grid-stride loops.

For fine-grained management over the precise attributes of your GPU, you’ll be able to depend on the lower-level, official CUDA Python package deal offered by Nvidia.