Prerequisite Concepts for Python in High-Performance Computing (HPC)
A High-Level Introduction to HPC Topics with Python, including Memory Layout, BLAS, LAPACK, SIMD and Threading
Introduction
This article introduces Python and key high-performance computing (HPC) concepts and libraries. It focuses on the main concepts rather than programming.
It covers:
- Compiled and interpreted languages
- Python libraries
- Memory, CPU caching, and memory locality
- C-style and Fortran-style ordering
- BLAS and LAPACK
- Parallel computing
- SIMD
- Threads, processes, and the GIL
Statement on AI
This article was not generated using AI. AI was used only to improve the grammar and sentence structure of my writing and to create helpful diagrams.
Compiled and Interpreted Languages
Compiled languages such as C#, C++, and Java convert code into instructions that can be executed efficiently by the computer.
Interpreted languages such as Python and JavaScript execute code through an interpreter. They are flexible and allow code to be run and modified easily. Python is also dynamically typed, so variable types do not need to be declared explicitly.
Although Python can be slower than compiled languages, it is widely used in HPC because of its powerful libraries and tools.
The standard Python implementation is CPython, which is written primarily in C. Other implementations include PyPy, which uses JIT compilation; Jython for the Java platform; and IronPython for .NET.
Python Libraries
The following provide tools for scientific computing, data analysis, parallel computing, machine learning, and visualisation:
Python is widely used in machine learning because it makes it easy to load data and build computational pipelines.
NumPy provides core tools for numerical computing, especially arrays. SciPy builds on NumPy by adding scientific tools for tasks such as signal processing, numerical integration, and optimisation.
Memory and Addresses
Memory can be thought of as a consecutive collection of storage locations, each with an address.
For example, memory might contain locations with addresses 0, 1, and 2. If we ask the computer for the value stored at address 2, it accesses that location and retrieves the data stored there.
Byte address: 0 1 2 3 ...
+-----+-----+-----+-----+-----+
Memory: | | | | | ... |
+-----+-----+-----+-----+-----+An array can be thought of as a consecutive region of memory with a starting address, a size, and offsets used to access its elements.
Array element = start address + offsetFor example, suppose an array of five elements begins at memory address 10. The elements could be accessed conceptually as:
Start address
10
|
+--------------+---------------+
| | |
v v v
Offset 0 1 2 3 4
| | | | |
v v v v v
+----+ +----+ +----+ +----+ +----+
Address | 10 | | 11 | | 12 | | 13 | | 14 |
+----+ +----+ +----+ +----+ +----+
10+0 10+1 10+2 10+3 10+4Therefore, when we use array indexing, the computer uses the array’s starting location and the requested offset to determine where the required element is stored.
Many programming languages use zero-based indexing, meaning the first element of an array has index 0 rather than 1. This comes from the relationship between arrays and memory addresses.
An array index can be thought of as an offset from the array’s starting memory address:
- The 1st element has an offset of 0
- The 2nd has an offset of 1
- The 3rd has an offset of 2
Programming languages provide different ways of allocating memory. For example, in C++, memory can be allocated using mechanisms such as new or malloc. A program requests a certain amount of memory, and the system provides access to an appropriate region of memory. The underlying idea is similar across different languages:
- Request a region of memory, specifying the type and number of elements required
- Reserve it for the program
- Access the data stored within that region
CPU Caching and Memory Locality
CPUs are extremely fast, but accessing data from main memory is comparatively slow. If the CPU had to retrieve every individual value directly from memory before processing it, memory access could become a major performance bottleneck.
Modern CPUs address this problem with caches, including L1, L2, and L3. When the CPU accesses a memory location, it may also load nearby data into its faster cache because it's likely to be needed next. This is why storing related data close together in memory can significantly improve performance.
First access
|
v
+------+------+------+------+------+
Memory: | 10 | 11 | 12 | 13 | 14 |
+------+------+------+------+------+
\______ CPU cache loads ______/
nearby memory locationsPython Lists vs Arrays
Python does not have a built-in array type that behaves like the contiguous numerical arrays used in languages such as C++. A Python list is different because it can contain objects of different types.
For example, a list could contain an integer, a string, and a floating-point number. Internally, the list essentially stores references to objects that may be located in different parts of memory.
A numerical array, by contrast, normally contains elements of the same data type stored contiguously in memory. This improves spatial locality and allows the CPU cache to work more efficiently.
Multidimensional Arrays and Memory Layout
Although we think of matrices as two-dimensional arrays, their elements are ultimately stored in a one-dimensional sequence of memory locations. The multidimensional structure describes how we interpret that sequence.
There are two important memory layouts:
- C-style (row-major) ordering stores the elements of each row consecutively.
- Fortran-style (column-major) ordering stores the elements of each column consecutively.
The choice of layout affects performance because accessing consecutive memory locations makes better use of the CPU cache.
A multidimensional array can therefore be viewed as a one-dimensional block of memory together with information describing its dimensions. For example, 100 consecutive elements could be interpreted as a 10 × 10 matrix without necessarily copying the underlying data.
C-Style vs Fortran-Style Ordering
The most efficient way to traverse an array depends on its memory layout. With C-style ordering, traversing elements in row order generally provides good memory locality. With Fortran-style ordering, traversing elements in column order provides good locality.
Warning: Ignoring the memory layout can reduce performance because the program may repeatedly jump between distant memory locations instead of accessing consecutive locations. This reduces the effectiveness of CPU caching.
C-Style (Row-Major)
In C-Style Ordering, a matrix is stored in memory row by row.
where A[i,j] is the address of the element in row i, column j.
Get George Marklow’s stories in your inbox
Join Medium for free to get updates from this writer.
For example,
Matrix:
+---+---+
| 4 | 7 |
A = +---+---+
| 1 | 6 |
+---+---+
is stored in memory as follows:
Index: 0 1 2 3
+-------+-------+-------+-------+
Value: | 4 | 7 | 6 | 1 |
+-------+-------+-------+-------+
-----> memory order ----->Fortran-Style (Column-Major)
In Fortran-style Ordering, the matrix is stored in memory column by column.
For example,
Matrix:
+---+---+
| 4 | 7 |
A = +---+---+
| 1 | 6 |
+---+---+is stored in memory as follows:
Index: 0 1 2 3
+-------+-------+-------+-------+
Value: | 4 | 1 | 7 | 6 |
+-------+-------+-------+-------+
-----> memory order ----->BLAS and LAPACK
BLAS (Basic Linear Algebra Subprograms) defines routines for common linear algebra operations, providing a standard interface that different vendors can optimise for their hardware.
For example, implementations such as Intel MKL can take advantage of specialised CPU features to accelerate calculations. This means that software can use the same BLAS interface while benefiting from a highly optimised implementation for the available hardware.
BLAS routines are commonly divided into three levels:
A rough rule-of-thumb way to remember what O(n), O(n²), and O(n³) algorithms might look like
- O(n): There is one loop that iterates through the vector once — that’s it!
- O(n²): There are typically two loops, one nested inside the other.
- O(n³): There are typically three loops, one nested inside another nested loop.
LAPACK (Linear Algebra Package) builds on BLAS and provides higher-level linear algebra routines. In other words, LAPACK uses BLAS fundamental linear algebra operations to solve bigger linear algebra problems. These include
- solving systems of linear equations
- eigenvalue problems
- matrix factorisations, i.e. decomposing matrices
- singular value decompositions (SVDs)
On a personal note: the analogy I like to use for BLAS vs. LAPACK is to think of NumPy vs. SciPy, where BLAS is NumPy and LAPACK is SciPy. SciPy is built on top of NumPy but has more functionality and does a lot more routines that build on BLAS.
There are variations, such as:
Parallel Computing
Modern computing is very parallel. This really took off in the mid-2000s, greatly improving performance.
Modern CPUs support parallelism through multiple cores and instructions that perform the same operation on multiple pieces of data at once.
GPUs are designed for highly parallel workloads, using thousands of processing cores to perform many calculations at the same time.
Parallelism can also involve connecting multiple computers and sharing calculations between them. This approach is common in HPC clusters.
SIMD
At this point, it is useful to move up the hierarchy and look at the CPU itself.
CPU-Level
SIMD Acceleration, a form of parallelism, stands for Single Instruction Multiple Data.
Suppose a CPU has cached some data. When it accesses one memory location, it may also load nearby memory locations into the cache because they are likely to be accessed next.
The CPU may need to perform the same calculation on several values. SIMD allows it to process these values at the same time using vector registers and instructions.
For example, AVX (Advanced Vector Extensions) is a SIMD instruction set. A CPU must support AVX to use it. Other instruction-set extensions, such as SSE2, provide additional features for improving performance.
SIMD uses vector registers, which can hold several values at once. For example, AVX provides 256-bit registers that can hold eight 32-bit floats or four 64-bit doubles. The CPU can then perform the same calculation on all these values at once.
SIMD instructions often describe the operation and data type. For example, VADDPD adds packed double-precision values, while VSQRTPD calculates their square roots.
Consider adding two arrays, a and b. Without SIMD, the CPU will add each pair of values separately:
Array a Array b Result
+------+ +------+ +-------------+
| a[0] | + | b[0] | --> | a[0] + b[0] |
+------+ +------+ +-------------+
| a[1] | + | b[1] | --> | a[1] + b[1] |
+------+ +------+ +-------------+
| a[2] | + | b[2] | --> | a[2] + b[2] |
+------+ +------+ +-------------+
| a[3] | + | b[3] | --> | a[3] + b[3] |
+------+ +------+ +-------------+With SIMD, the VADDPD instruction can perform four additions at once. In this simplified example, instead of four separate instructions taking four cycles, SIMD operates with one instruction in one cycle:
Array a Array b Result
+------+ +------+ +-------------+
| a[0] | | b[0] | | a[0] + b[0] |
+------+ +------+ +-------------+
| a[1] | | b[1] | | a[1] + b[1] |
+------+ +------+ +-------------+
VADDPD( | a[2] | , | b[2] | ) -> | a[2] + b[2] |
+------+ +------+ +-------------+
| a[3] | | b[3] | | a[3] + b[3] |
+------+ +------+ +-------------+This gives four times the throughput, allowing the CPU to process more data in the same amount of time.
In summary, SIMD is a low-level CPU feature that is not normally used directly from Python. However, several Python libraries can take advantage of SIMD:
- NumPy: can benefit from SIMD through optimised numerical routines, including BLAS implementations.
- NumExpr: evaluates numerical expressions using optimised, vectorised operations.
- Numba: uses JIT compilation and can generate SIMD instructions for suitable Python code.
Data Parallelism
SIMD usually provides parallelism within a single CPU core. Multiple CPU cores can provide further parallelism by processing different parts of a calculation at the same time.
Using multiple CPU cores can significantly improve the performance of computationally intensive applications.
Thread and Multi-Process Execution
There are two main execution units: processes and threads.
A Process is an independent running program with its own memory. Processes are isolated from each other, so communication requires mechanisms provided by the operating system or other libraries (e.g. MPI). For example, running several Python programs creates several Python processes. Think of a process as having its own ‘neighbourhood’.
A Thread is an execution stream within a process. Multiple threads in the same process can share and access the same memory, making communication between them easier.
A Word on the GIL
You may have heard of the Global Interpreter Lock (GIL), a legacy feature that normally allows only one thread to execute Python bytecode at a time. For example, if a program has four threads, only one can execute Python code at any given moment.
The GIL simplifies parts of CPython, particularly memory management and garbage collection. A thread must acquire the GIL before executing Python bytecode, preventing other threads from doing so at the same time.
The main disadvantage is that CPU-bound programs may not benefit from multithreading because the threads must compete for the GIL.
Note: At the time of writing, PEP 703, which makes the GIL optional in CPython, has been implemented. I recommend reading the PEP for a more detailed explanation of these changes.
Thanks for reading! Let me know what you think in the comments section below, and don’t forget to subscribe 👍
