toby

In August last year I gave a talk at the 2017 HPC Advisory Council conference in Perth, where I showed some of work I had been doing on optimising a seismic code for large scale seismic wave propagation simulations at The University of Western Australia. It was the first time I had organised some volunteers to help me demonstrate a point, with some hilarious results!

The performance enhancements I contributed amounted to a factor of 2x increase in the performance of the code. This ended up saving our research group about 4 million core hours of compute time. In this talk I also do a direct comparison between the performance of AWS cloud computing nodes and the Pawsey Supercomputer Magnus.

PDF version of article

It can be confusing as to how CUDA hardware threads (CUDA cores) are mapped to the available
work. This article shows how to disentangle the CUDA hardware thread mapping.

Blocks

When a CUDA kernel executes, it executes hardware threads in batches of 32 threads called warps.
Threads are referenced according to 3D rectangular cuboid called a block, with axes (x,y,z). At present blocks are limited to a maximum 1024 threads. The maximum length of a block in the z direction is 64 threads and any of the axes for the block may have a “length” of 1 thread. A common choice for the dimensions of a block is (16,16,1) so that a block has 256 threads which is a multiple of the 32-thread warp size.

Grid

Blocks are referenced according to a 3D rectangular cuboid structure called a grid, also with axes (x,y,z). The maximum number of blocks along a dimension of the grid is 65536, and the minimum number is 1. Parameters that you specify to the kernel include both the dimensions of the block, and the dimensions of the grid.

Structures available during kernel execution

During a CUDA kernel execution, the following structures are available to the kernel code:

• gridDim.x, gridDim.y, gridDim.z; the number of blocks along the x, y, and z dimensions of the grid
• blockIdx.x, blockIdx.y, blockIdx.z; the zero-based index of the block along the
  x,y, and z dimensions of the grid
• blockDim.x, blockDim.y, blockDim.z; the number of threads along the x, y, and z
  dimensions of the block.
• threadIdx.x, threadIdx.y, threadIdx.z; the zero-based index of the threadalong the x, y, and z dimensions of the block.

Mapping order

The order of threads in blocks and the order of blocks in grids are both arranged according to a
FORTRAN-style column-major format. That is, they are arranged along x axis first, then the y axis, and finally the z axis. Care must be taken when CUDA kernels are applied to multidimensional arrays whose ordering is in row-based (C-style) format, as the thread ordering is opposite to the layout of the
array. Note: The order in which the blocks and threads are executed in is undefined and cannot be
relied upon in the execution of the kernel.

Example mapping

In the diagram below is a mapping where the grid has size (5,4,1), that is, it has 5 blocks in the x
direction, 4 blocks in the y direction, and 1 block in the z direction. Each block is of size (5,5,1) with 5 threads along the x and y directions, and 1 thread along the z direction. At the grid level (on the left), the tuple for each block is the 3D index, e.g. (0,0,0), and below the 3D index is the 1D index of each block, e.g. [0]. At the block level (on the right), a similar indexing scheme applies, where the tuple is the 3D index of the thread within the block and the number in the square bracket is the 1D index of the thread within the block. During execution, the CUDA threads are mapped to the problem in an undefined manner.

In the diagram below we show randomly completed threads and blocks as green to highlight the fact that the order of execution for threads is undefined.

3D index of a block relative to its grid

This is just (blockIdx.x, blockIdx.y, blockIdx.z)

3D index of a thread relative to its block

This is just (threadIdx.x, threadIdx.y, threadIdx.z)

Obtaining a thread’s 3D index relative to the grid

In order to get a thread’s 3D index relative to the grid, one has to multiply the index of the block by the number of threads along that dimension of the block, and add the component of 3D thread index. In a CUDA kernel it looks like this:

int x_index = blockIdx.x*blockDim.x + threadIdx.x;
int y_index = blockIdx.y*blockDim.y + threadIdx.y;
int z_index = blockIdx.z*blockDim.z + threadIdx.z;

1D index of a block relative to the grid

The 1D index of a block relative to the grid is given by the dot product of the block indices with the vector containing the number of blocks to stride over in each dimension.
blockId_grid = [blockIdx.x, blockIdx.y, blockIdx.z]·[1, gridDim.x,gridDim.x*gridDim.y]

1D index of a thread relative to its block

Similarly, the 1D index of a thread within a block is obtained as the dot product between the vector of indices and the vector containing the number threads to stride over in each dimension.
threadId_block = [threadIdx.x, threadIdx.y, threadIdx.z]·[1,
blockDim.x, blockDim.x*blockDim.y]

Obtaining a thread’s 1D index relative to the grid.

We obtain the 1D index of a thread relative to the grid by first obtaining the 1D index of the block,
multiplying it by the number of threads per block, and adding the 1D index of the thread within the block.
threadId_grid = blockId_grid *(blockDim.x*blockDim.y*blockDim.z)
+threadId_block

In C code within a kernel this could take the form:

size_t blockId_grid = blockIdx.x + blockIdx.y*gridDim.x + blockIdx.z*gridDim.x*gridDim.y;
size_t threads_per_block = blockDim.x*blockDim.y*blockDim.z;
size_t threadId_grid= blockId_grid * threads_per_block +threadIdx.x + threadIdx.y*blockDim.x +threadIdx.z*blockDim.x*blockDim.y;

A slight performance improvement might be obtained by using int instead of size_t, however you do risk integer overflow if the number of elements exceeds 2147483647 or approximately 8.6GB on the GPU for 32-bit floating point arrays.

If any of the dimensions of the grid have size 1 (e.g. gridDim.z==1), then the corresponding
blockIdx term for that dimension (e.g. blockIdx.z*gridDim.x*gridDim.y), is always
zero and you might consider removing that term from the calculation of blockId_grid.
Similarly, if any of the dimensions of a block have size 1 (e.g. blockDim.y==1), then the
corresponding threadIdx term (e.g threadIdx.y*blockDim.x), is always zero and may be
removed from the calculation of threadId_grid.

Reference:

Information in this sheet was sourced from the CUDA C programming guide, located
here.

This work is licensed under a Creative Commons Attribution-ShareAlike 3.0 Australia License.

 

In this video I explore how Python and its associated packages offer a holistic solution for your business data processing needs.

According to a PriceWaterhouseCoopers report, 74% of top Australian CEO’s “say availability of key skills – including technology skills – is a threat” to their business. One way to address this shortage of skills is for companies to train their staff in STEM (Science Technology Engineering and Mathematics) skills. The same report estimates that moving just 1 percent of the workforce to STEM roles would add $57.4 billion to Australia’s gross domestic product.

So how can we train staff in STEM skills? A practical way is to equip them with powerful IT tools that lead to innovation.  Most of us are aware of Python as a computer language, but are we aware of the awesome power it has for STEM applications?

Here are 10 powerful STEM abilities that Python unlocks

1. Write code that is both powerful and easy to read

Python is designed to be a general-purpose programming language that is easy to write and understand. In Python the same algorithm can be expressed in far fewer lines of code than other languages such as C++. Fewer lines of code means less time is spent in software development and fewer chances for errors to creep in. During execution of a Python program, the Python interpreter is careful at error and sanity checking, and when programming errors do occur, Python is very helpful in pointing out where the errors are.

Careful sanity checking can slow Python performance when compared to other languages like C. However, when developer time is more important than execution time – such as in exploratory data analysis, prototyping ideas, and testing hypotheses – the productivity benefits of writing simple and powerful code far outweigh potential drawbacks.

2. Write code that is platform-independent

When a Python script runs it is compiled to a platform-independent language called bytecode. Bytecode instructions are then translated into machine instructions by a program called a Python Virtual Machine. With this design the same Python code can run on multiple platforms such as Windows, MacOS and Linux, and across multiple architectures such as ARM, x86 or PowerPC. You can, for example, develop a Python program on Windows and run it  on a Linux machine. Platform-independent code is also great for cloud computing where compute resources and software environments are frequently subject to change.

3. Join complex systems

Python speaks a wide range of protocols and data formats. You can for example, use Python to access databases using SQL, and send and receive data over the network using libraries like socket, MPI4PY, and ZeroMQ. Python can read and write binary and ASCII files and has support for data formats like HDF5. Python can run programs both locally and over remote connections using SSH. The flexibility that Python provides makes it ideal for joining complex systems.

4. Do sophisticated math with powerful science libraries

Thanks to the Numpy and Scipy packages, Python has first rate multidimensional arrays with a comprehensive set of linear algebra routines. These linear algebra routines call optimized low-level math routines for performance. The low-level routines are interchangeable, this means you can use highly optimised commercial libraries like Intel’s MKL.

In addition to having great array support, Scipy and Numpy support a range of capabilities like statistics and random number generation; interpolation and numerical integration; N-dimensional Fast Fourier Transforms; optimisation, image processing, and special functions like Bessel functions. You can even use Einstein notation for math with tensors!

5. Export publication-quality plots

With the Matplotlib package, Python supports almost limitless versatility in making and arranging beautiful 2D and 3D figures and exporting them to images. The plotting backends support a wide variety of image formats, including popular ones such as SVG, PDF, EPS, TIFF, PNG.  The majority of plots in my supernova research paper were made with Python. Plotting packages like Bokeh and Plotly provide 2D and 3D plots that are rendered by a web browser. Both Bokeh and Plotly give you the tools to make beautiful dashboards.

6. Provision cloud resources

The big three cloud providers Amazon, Google, and Microsoft all have software development kits that allow you to provision and maintain cloud computing resources from within Python programs. This is great because you can easily integrate cloud provisioning into Python workflows.

7. Explore machine learning and artificial intelligence algorithms

Python provides an easy on-ramp to machine learning and classification algorithms with the well-documented Tensor Flow, Theano, and Scikit-learn libraries. This is useful, for example, in image and speech recognition, clustering data, creating artificial neural networks, and data mining.

8. Get more from your data with analytics and big data capabilities

The Pandas package provides an interface for reading and writing Microsoft Excel spreadsheets and powerful tools for working with labelled and time-series data. Data analytics and business intelligence benefits enormously from Pandas, as it provides easy lookup and transformation of data. The big data processing framework Apache Spark has a Python interface so you can run distributed algorithms over large scale clusters.

9. Conduct reproducible research

Jupyter notebooks have the ability to store HTML content, equations, and code in one file  for easy sharing with colleagues and collaborators. This is great because readers can execute the code in the notebook and reproduce the results of the author, thus saving time. Python works great with Jupyter; 2D and 3D visualisations can be incorporated directly into notebooks.

10. Get it for free

The best reason is that Python is free and open source. The intellection property rights are held by the Python Software Foundation, which releases the software under a GPL-compatible open source license. This means you are free to distribute and use Python to an unlimited number of machines.

Conclusion

Python and its associated libraries bring a wide variety of capabilities that make it a powerful tool for STEM. In 2008 I started using Python for all my data processing needs and it has benefitted my research enormously. I am sure it can do the same for you.