What it takes to transpose a matrix

Classical CPU architecture is a poor choice for performing matrix-oriented computations. Developers have to spend much effort to create efficient algorithms even for the problems appearing trivial on the surface. Matrix transpose is probably the most striking representative of such problems. It exposes numerous performance-related issues, such as high memory latency, shortcomings of cache organization, and failure of the compiler to automatically vectorize the code. Performance of more mathematically complex matrix operations is usually bounded by slow math instructions (such as div or sqrt), which mask the majority of other issues. Matrix transpose, on the other hand, being stripped of any slow instructions, is open to all kinds of inefficiencies. In this article we are going to gradually build a sequence of progressively more efficient implementations of matrix transpose, with the most sophisticated implementation being up to x25 times faster than the naive one. During each step we will locate the bottleneck, figure out what has caused it, and think of a solution to overcome it. This article is intended to serve as an introduction to optimizing matrix algorithms for x86_64, presented from the perspective of a real-world problem.

We will work with the following formal problem statement: given NxN src matrix consisting of 1-byte elements, and dst matrix of the same shape (non-overlapping with src), copy each src[r,c] element into dst[c,r]. Below is the straightforward C++ prototype:

Here you can see that I imposed some restrictions: matrices are of square shape, they are filled with one-byte elements, and transpose happens out-of-place. The most important restriction is that element size is one byte long. The reason behind it is that it will make things more challenging than with larger element sizes. Increasing element size to, let’s say, int64_t will create more pressure on memory throughput, making it primary bottleneck and therefore hiding other issues. In general, all imposed restrictions can be lifted without loss in performance. However, to create a truly "universal" transposer would require too much boilerplate code not fitted for the purpose of this article.

Numbers in square brackets after section names will indicate performance of the corresponding algorithm, expressed as average number of CPU cycles per single matrix element. They were measured for a matrix of size 2112x2112, selected to represent medium size. For larger matrices the difference in algorithm performance will be even more prominent. The hardware used for testing is Skylake 7700HQ CPU with Turbo Boost off.

Naive implementation is trivial and can be coded in less than a minute. Just scan src row by row, column by column, and write elements to the corresponding location in dst matrix.

Below animation demonstrates the progress of the above algorithm being applied to a small 22x22 matrix. Color denotes status of the element: already processed (green), current (blue), or not yet processed (red). Fading effect encodes which elements are currently in cache and their age since they were last accessed. In this toy simulation cache size was set to 10 entries with last-recently-used (LRU) eviction policy, where each cache line is 8 bytes long. Real x86_64 systems have caches of thousands of 64-byte long cache lines.

The implementation we currently have is the simplest possible and, as it usually happens, also the slowest one. We will use it as a baseline to compare more advanced implementations against it. But before moving forward, we would like to figure out the source of its inefficiency. To do it, let’s trace mentally the function step by step and look at what happens at hardware level. We are primarily interested in the memory subsystem. Since there are two logical streams of data — read from src and write to dst — we will analyze them separately.

While our expectations in general work good, for some N we can observe high spikes in running time. They occur usually around some "good" numbers, such as powers of two: 256, 2048, 4096. By looking into PMU counters we can detect that these spikes coincide with suspiciously bad cache behaviour — something we didn’t anticipate. What exactly could the be the cause?

The issue stems from cache organization. Usually we model cache as a black box, which is able to store K cache lines with last recently used eviction policy. In another words, if we access K different cache lines, then they all will happily settle down in cache, and only after accessing K+1’th distinct cache line some older cache line will have to be evicted. Such behaviour requires the cache to be fully associative, i.e. every slot in the cache must be able to store any cache line. In reality fully associative caches are rarely used because of performance issues — imagine, to search for a cache line, hardware would have to compare target address with all K (tens or even hundreds of thousands) cache slots simultaneously. Instead, archetypical cache is organized as a collection of multi-way sets. Logically, each set acts as a small, independent, fully-associative cache with capacity equal to number of ways (typically 4-16). Given index of the cache line, it can be stored only inside single set, which is computed using modulo operator. Example below shows the computations that are performed when you dereference a pointer with value 57690₁₀. First, cache line index is computed by dividing address by cache line size: ⌊57690₁₀/64⌋=901. And next, taking the remainder of dividing index by number of sets gives the set index: 901 % 64 = 5. Set number 5 is the only set where cache line backing address 57690₁₀ can be cached. Note that all other cache lines sharing the same remainder have to compete for the same set number 5. Such organization is very efficient. Since cache line size and number of sets in practice are both powers of two, the result of division and modulo operators is just subsequence of bits inside binary representation of the address, so in fact nothing needs to be "computed". And once the set index is known, only 8 ways must be checked simultaneously to find whether cache line at given address is cached or not.

For the majority of programs multi-set design is no different from fully-associative cache since sequence of accessed addresses in a random program is irregular. For a typical program, the majority of memory accesses are caused by walking through complex hierarchies of objects by dereferencing pointers into the heap. Hence all the cache sets are loaded approximately evenly. However, algorithms operating on multi-dimensional data are a notable exception from this rule. Majority of such algorithms walk along columns or planes of the data, thus inducing strided access pattern. Naive transposer is exactly the algorithm of this kind. Its inner loop fills in single column, thus accessing elements standing apart by N: ColumnBegin+i*N, i=0,1,2,…​ Strided access pattern is notorious for its interference with multi-set cache organization. Suppose that ColumnBegin=12345678₁₀, N=256, and number of sets is 64. The elements of this column are mapped to the following sequence of set indices:

You can see that the pattern repeats and that only 16 sets are involved instead of all 64. The remaining 48 sets have zero load. This effectively reduces cache size from 512 entries down to just 128. Since the latter is smaller than N=256, data overflows into L2, ruining performance for this specific N. For each level of cache there exist some values of N, for which one column of data should fit into cache but it doesn’t because of multi-set design. These bad N values may differ from one level of cache to another one because caches have different number of sets and ways.

Now it becomes clear that it is not enough to compare N with number of cache entries to predict performance. To ensure perfect caching we need to compute load of each of the sets by simulating first N iterations. If load for all sets is less or equal to number of ways, then Nx64 block will fit into the desired cache. If load exceeds number of ways for at least one set, then there will be overflow, leading to suboptimal performance. Of course, simulation is too cumbersome to do in practice. Fortunately, it is possible to avoid cache aliasing (how this issue is called) without knowing the exact geometry of caches. Let us state sufficient condition on N that ensures perfect caching: N must equal to cache line size multiplied by some odd number. All of N=64, N=192, N=6464, N=6592 will guarantee perfect caching for all cache levels. To prove the stated condition it is enough to show that for every two iterations i1, i2, such that 0 < |i1 - i2| < 64, accessed addresses are mapped to different cache sets.

The last transition follows from the fact that since NumSets is always power of two, then NumSets and OddNumber are co-prime. This is where we depend on the oddity. Final statement is true since 0 < |i1 - i2| < NumSets = 64. Net result is that if N equals to some odd number times cache line size, then any N consecutive iterations will access cache lines mapped to different sets.

How can the knowledge of which sizes are good and which are bad be of any help to us? We cannot require user to work only with matrices of good sizes — we have to support arbitrary-sized matrices. But what we can do is to require user to pad matrix up to an odd number of cache lines during memory allocation. If user wants matrix of, let’s say, logical size 4000x4000, then they should allocate matrix of physical size 4032x4032 instead. Top left block of size 4000x4000 would carry the actual data, while the remaining fringe is just a placeholder to ensure that the stride value is proportional to an odd number of cache lines. In terms of memory overhead requirement for padding is costly only for small matrices. With N growing to infinity, the overhead fades to zero very soon. For N=4000 it is less than 2%.

In addition to padding, we will also require that the beginning of matrices is aligned by cache line size (64). Together with requirement for padding this means that each row will be covered by integer number of cache lines. The regularity of such layout grants us a number of benefits: it is much easier to analyze; cache lines backing matrices become "matrix-private" in the sense that they carry only the data of the matrix and never the data outside of it; finally, but of no least importance, is that such layout is aesthetically pleasing.

After imposing new requirements and rerunning our previous examples, all timings become as expected. There are no more unpleasant spikes. For the sake of brevity, from now on, we will assume that matrix size N itself is proportional to an odd number of cache lines and that its beginning is aligned with the beginning of cache line.

In the naive algorithm, we have perfectly good read performance but very poor write performance. The latter is caused primarily by the fact that single Nx64 rectangle of dst is not able to fit into small L1d cache and eventually has to settle down in caches of deeper levels or RAM, from where it has to be reloaded 64 times. But also poor performance to a lesser extent is the result of the mechanics of write operation. Unlike read operations, which can discard cache line after the load is complete, write operations have to load cache line, modify it, and store it back into cache. They cannot throw it away. In this sense, writes are more expensive that reads in our algorithm. This leads us to a straightforward idea of reversing the order of scans. Let’s scan destination matrix in row-major order, while source matrix in column-major order.

We will call new algorithm the "reverse" one. Its code is exactly the same as for "naive" algorithm except that src and dst are swapped. Note that now it is dst matrix that is processed in a natural way: row by row, column by column.

By changing the order of scans we optimize for writes while sacrifiying reads. Now writes are perfect instead of reads. Since destination matrix is scanned in row-major order, a group of 64 consecutive dst elements is mapped to a single cache line. This allows for the cache line to be fully assembled from scratch in a store buffer and then flushed into cache, all of this without waiting for the load of the cache line to complete. Reads become bad, of course, for the same reasons as previously writes were. They are served from the fastest cache where Nx64 rectangle fits.

To see the difference between naive and reverse algorithms, let’s apply them to a matrix of size N=2112 and look at PMU counters related to caching. Table below lists such counters for transposing 1925 different N=2112 matrices, covering in total 8GiB memory region. Two right columns contain counter values normalized by total number of elements (2112x2112x1925). Normalized values close to integers are of primary interest to us.

First let’s compare writes. As was already mentioned before, l2_rqsts.rfo_{miss|hit} count L2 misses and hits with respect to write operation. Naive algorithm demonstrates one hit per each element and no misses, proving that for each processed element cache line has to be fetched from L2. This is exactly where we expect to find it, since single Nx64 block cannot fit into L1 but fits into L2: 32KiB < 2112x64=132KiB ≤ 256KiB. On the other hand, both counters are zero for the reverse algorithm, indicating that there was no need to search L2 or higher-level caches at all. Another way to see that reverse algorithm significantly improves write performance is resource_stalls.sb counter. Roughly speaking, it measures number of cycles when execution couldn’t advance because store buffer had no free entry. Its value is approximately 40% of all cycles for the naive algorithm, but is nearly zero for the reverse algorithm.

Comparing read performance is slightly more complicated. Presence of three L1 hits per each element and no L2 misses for the naive algorithm proves our reasoning that all data comes from L1. But why three loads and not one as we expected? To answer this question, let’s look at the generated machine code. Below is the inner loop of the naive transpose algorithm generated by GCC-10.

It turns out that in lines 9 and 11 there are extra loads. They reload addresses stored in src._data and dst._data into registers during each iteration of the algorithm. Since these values are constant, they obviously are always served from L1 and do not have any noticeable negative impact on performance. However, they increase counter values in the same way as truly heavy data loads do. That’s why we observe two extra loads from L1d per each processed element than we expected. This happens for both algorithms since generated code is identical with the exception that references to src and dst are swapped. Of course, different compilers and even different options of the same compiler may result in different generated code, with bigger or fewer number of loads per each element.

Now with extra loads demystified, we can clearly see that performance of reads mirrors that of writes. All data loads of the naive algorithm are served from L1 with no accesses to deeper level caches. But for the reverse algorithm data loads are served from L2. The latter can be deduced from l1_miss≈l2_hit≈1.0 and no accesses to L3 (l2_miss≈0.0).

To conclude the analysis, let’s look into the most important metric — number of cycles per element. Here we can see that reverse algorithm is faster by 1.4 times. Not bad taking into account that the only thing we had to do is to swap src and dst references. Although read and write performance in terms of events per element are entirely symmetric for two algorithms, write "events" turn out to be more expensive. This explains the superiority of reverse algorithm compared to the naive one.

Both naive and reverse algorithms are inefficient because the way they access memory causes continual eviction of L1d cache lines into deeper-level caches (or into RAM) and further reload of these cache lines back into L1d. In other words, working set is too huge compared to L1d cache. For both algorithms it equals to 64×N+ε. As soon as N exceeds size(L1d)/64≃512, continual eviction and reload of cache lines takes place. During this time program flow is blocked by large number of cycles during each iteration, leading to poor performance. The larger N, the deeper memory subsystem layers are accessed and the more prominent degradation becomes. If we want to make any progress on optimization, we need to somehow make working set fit into L1d no matter how large N is.

Luckily, the transpose operation preserves 2D data locality. If a group of elements are located close in source matrix, then they will be close in destination matrix. They are not scattered all over the dst matrix, but closely reside in a submatrix of the same shape and size. Not all operations are so good. For example, when squaring a matrix, each output element depends on one full row and one full column of the input matrix. In another example — computing partial sums — each output element is the sum of all i⋅j preceding elements. But the majority of matrix operations preserve 2D data locality, since matrix as a data structure is designed to represent spatially related data. Transpose also falls into category of such operations.

Data locality is a very powerful property because it allows to split matrix into arbitrary non-overlapping blocks and process one block at a time. The result of a single block transpose will be also a block of the same shape and size, therefore allowing caching to do its best. If block size is small enough, then it can be processed entirely within L1d boundaries without overflowing into deeper-level caches — exactly what we need. Image below demonstrates how we can split 12x12 matrix into 4x4 blocks. Note an obvious fact that locations of the blocks themselves are also transposed. For example, when we are working on top-right AC block, we are transposing its elements and writing them into bottom-left block of the destination matrix. More theoretically it can be said that what we are doing is a composition of 1+(12/4)(12/4)=10 transpositions. One of them is the outer transposition of the 3x3 block matrix itself, where each element is a single 4x4 block. All other 9 transpositions are byte-wise transpositions of the respective blocks. We can perform transpositions in any order we like among all 10! variants and will get identical results in the end. In practice we will process blocks one at a time and write elements to the final positions in dst matrix, thereby implicitly performing outer transposition.

Guided by the idea, let’s treat matrix as consisting of cache aligned blocks of size 64x64, so that every row of a block maps to exactly one cache line. With such setup, no matter how large original NxN matrix is, the working set of any transposer is bounded from above by 2×64×64+ε = 8KiB+ε < size(L1d). This bound consists of one 64x64 block of the source matrix, one 64x64 block of the destination matrix, and some small, constant-sized volume of memory for the algorithm itself.

Source code for block transpose requires four nested loops. Two outer loops iterate over the blocks, while two inner loops iterate over the elements of a single block. The order of scans (row-major or column-major) is of no particular importance thanks to small working set, with only marginal difference in performance between all four possible options. The piece of code below scans blocks in row-major order, while the elements of a single block in a column-major order, exactly as in the animation.

The efficient usage of cache should be evident from the animation. Now reading can cause expensive stalls only during the first 64 iterations of transposing each block, when the block is being loaded into cache, while the data for the remaining 63*64 iterations is cheaply served from L1d. Situation with writes is similar. Rows of dst block are fully assembled in store buffer and flushed at some point later asynchronously to code execution.

As before, per-element PMU counters can be used to confirm our findings. If we apply all three algorithms to a matrix of N=2112, we can observe that block algorithm is the only one not affected by any cache misses. All L2-related counters are zero, meaning that there was no need to send requests to L2 or deeper caches. This is true for both reads (mem_load_retired.l2_{hit,miss}) and writes (l2_rqsts.rfo_{hit,miss}). In this sense block algorithm inherits all the benefits of naive and reverse algorithms without inheriting any of their drawbacks.

Now let’s rerun the tests for matrices larger than N=2112. We can now clearly see what makes block algorithm so good. Its performance is constant for both small and big matrices — a property of the algorithms so much welcomed in practice.

It is hard to overestimate the importance of block decomposition for the matrix algorithms. Like all such algorithms, block transposer provides constantly good performance, and the code listing is still fewer than 20 lines of non-architectural code. Block transpose is probably the method of choice for the majority of applications. It is efficient, simple, and doesn’t depend on specifics of particular CPU architecture. Upcoming methods are even faster, but they lack these properties.

In our next improvement towards reducing memory latency we will stop relying on automatic prefetching and replace it with manually issued hints to prefetch data. In optimization parlance the former is usually called "hardware" prefetching to emphasize that prefetch requests are initiated by hardwired algorithms, while the latter is called "software" prefetching.

The benefits of hardware prefetching are many. Its killer feature is that it works fully automatically without any input from the coder as we have already seen. Of course, you still have to know basic dos and don’ts. For example, unnecessary using linked lists instead of arrays will destroy prefetching altogether. But for an average program without blatant flaws hardware prefetching works well. Another good thing about hardware prefetching is that, since it is implemented in hardware, it can take into account (at least theoretically) all the architectural details of specific CPU/RAM setup, such as cache sizes, parallelism, memory timings, etc. Moreover, behaviour of hardware prefetching can change dynamically along with program execution. Imagine, for example, that prefetcher mispredicts that a sequence of addresses a1,a2,a3,…​ will be accessed in the future and starts aggressively to preload them. With the help of counters implemented in the cache, it can figure out its mistake after some time, i.e. the fact that preloaded data is not accessed by the program. In this case prefetcher can stop prefetching this sequence of addresses to avoid cache pollution and to give opportunity to other streams of requests competing for memory access. Such dynamic algorithm can be implemented with not-so-complicated counters in hardware, but implementing it in software would be cumbersome to say the least.

And yet hardware prefetching is far from perfect. Its scope of prefetching is usually limited to 4KiB memory page, meaning that it won’t preload data that is beyond already accessed page. Another downside is that number of independently tracked streams is not very big. In our block transpose algorithm we need 64 of them to fully preload next block ahead of time. But probably the largest downside is that hardware prefetcher, working on a very low level, cannot grasp the high-level idea of the algorithm. For example, in block transpose algorithm it is quite improbable that it would predict correctly and preload block in position (1,0), while CPU is crunching top right block (0,nr_blocks-1). In general, in matrix algorithms, where operations are quite simple but are applied to large volume of data, it makes sense to try to guide prefetcher by issuing prefetch commands manually. Sometimes it improves performance, sometimes it doesn’t.

Prefetch instructions are architecture-dependent and are usually called through intrinsics — C/C++ library functions provided by the compiler. In case of x86_64 we will call __mm_prefetch(char* addr, _MM_HINT_NTA), which is translated by all popular C/C++ compilers into one of the corresponding x86_64 PREFETCHh instructions. This function takes two arguments: address and locality hint. Address can be any, not necessarily aligned with cache line boundary, and it can be even an invalid address, in which case all access violation errors are silently suppressed. The latter feature is useful since it allows to avoid conditional logics when working near boundary of vectors/matrices. Second argument is a locality hint that indicates our intention towards the data. _MM_HINT_NTA ("non-temporal") indicates that we are going to work with data only briefly, thereby helping CPU to take appropriate steps to evict data soon and avoid cache pollution. Note that issuing __mm_prefetch() is merely a hint and does not guarantee any solid outcome. CPU is free to entirely ignore it. Explanation of different locality hints (T0, T1, T2, NTA) also should be taken with a grain of salt due to vast number of subtleties in modern CPUs. Treat locality hints more like an opaque enumeration with no clear semantics attached. In general, the most reliable approach is to try different combinations of __mm_prefetch() placements and locality hints and to see which one improves performance the most.

Now let’s apply software prefetching to block transpose algorithm. While we are working on a single block, we would like to initiate preload of the next block in the background. Our changes are the following. First, we include immintrin.h header to get access to x86_64 prefetch functions. Second, each time we start processing new block, we also compute prf_origin — pointer to the next block of the src matrix. Finally, we issue __mm_prefetch per every row of the next block, 64 in total. It is quite doubtful that CPU would be able to process them efficiently if all of them were issued at once. So, instead we insert prefetch commands before processing every row of the current block. Such placement scatters prefetch commands evenly along the instruction sequence.

Table below compares performance of four previously described algorithms. Impact of prefetching is easier to demonstrate using huge matrix. Therefore testing was performed by transposing three different N=46400 matrices filled with random values. In addition to cycles and cycles per element (cpe) metrics I also included cycle_activity.stalls_total. The latter is PMU metric that counts number of cycles when execution was stalled for any reason. Since we know a priori that we have big issues with memory access, we can expect that most stalls in the compared algorithms are caused by the data being requested from memory subsystem but not delivered yet, and therefore CPU has to idle.

According to the table the largest gain in performance was achieved by switching to block algorithm. This is understandable because block algorithm makes caching perfect. Impact of software prefetcher is not so dramatic. However, given that it costs virtually nothing to implement, cpe reduction by 14% is also quite pleasant.

So far we were dealing with high latency of memory subsystem. Exploiting block property helped to solve this problem quite well. Since now on we will stick with the idea of processing one 64×64 cache aligned block at a time and will focus on improving the algorithm of transposing such blocks.

Block transpose algorithm we designed so far makes only one load and one store per every element. This is the lowest as it can be provided that we limit ourselves with processing only one element at a time. Now recall that elements of our matrix are bytes, and, therefore, loads and stores we perform act on single bytes. This doesn’t fully utilize the power of x86_64 since its registers and almost all instructions are 64-bit. Using narrower registers and instructions (8/16/32 bits) typically costs the same as if using native 64-bits instructions. So, let’s try to design an algorithm that will process multiple elements simultaneously.

Since single 64-bit register can hold 8 elements at once, it is tempting to apply block decomposition again and design an algorithm that is capable of transposing 8x8 matrices. Such an algorithm would start from loading rows of the input 8x8 matrix into eight 64-bit variables, then it would perform some black-box computations on these variables, finally producing eight 64-bit variables holding rows of the transposed matrix. These rows would be directly written into their proper positions of dst matrix. Now, given such black-box algorithm, if we decompose src matrix into 64x64 blocks first (to optimize caching), and then decompose every 64x64 block into 8x8 blocks, and then apply 8x8 transposer to each of these blocks, we will get properly transposed dst matrix.

Designing efficient 8x8 transposer with 64-bit variables provides some challenges, though. Every row of the destination matrix depends on all rows of the source matrix. It would be very inefficient to use brute force transposer that would compute every output row independently: to compute single output row, every input row must be masked, shifted properly and ORed with the result. Overall complexity of such transposer would be proportional to 8x8, thereby providing no improvement in performance over the naive algorithm. A better solution is to use block decomposition again; now it will be the third time when we make use of it.

The idea is to reinterpret 8x8 matrix recursively as 2x2 block matrices. This requires log₂(8)=3 levels of block decomposition. At first we reinterpret 8x8 matrix as 2x2 block matrix with blocks of size 4x4 as its elements. Next we reinterpret each of these blocks as 2x2 matrix with 2x2 blocks as its elements. Finally each 2x2 block is a 2x2 block matrix itself. Transposing all 2x2 block matrices at all three levels (1+4+16=21 matrices in total) is equivalent to transposing original 8x8 matrix thanks to block property. We are free to select any order in which to process block matrices: small matrices first, large matrices first, or any mix of thereof. Image below demonstrates the progress of the algorithm "bottom-to-top", starting from the smallest matrices.

Note an important fact that at every level we have only 2x2 block matrices, although with blocks of different sizes. This is beneficial for us because every output row of 2x2 block transpose depends only on two input rows, and this holds true for every block size. Look, for example, at row b0. It can be easily computed by applying bitwise operations to rows a0 and a1:

Similarly, row d0 depends only on rows c0 and c4 but requires different masks and shift values. In this manner we can do all computations row by row, thereby working on multiple matrices in parallel. Total complexity of such approach is α*8*log₂(8), which is significantly faster than β*8*8 complexity of the naive transposer, even if we take into account that α is slightly larger than β because of a more complicated instruction sequence per each step. We will call our new transposer Vec64 to highlight that we process multiple elements in parallel. Obviously, the benefits grow as word width increases. For 32-bit architecture (decomposition into 2 levels) the benefits over naive transposer would diminish, while for hypothetical 128-bit architecture (4 levels) they would be even more prominent.

Another challenge of our new algorithm is where to place prefetch calls. Vec64 is not very homogenic: it is hard to split it equally into 64 equal parts without complicating code too much. One of the solutions is to issue prefetch calls in batches of eight, hoping that such size is not very big to cause any performance issues.

Full listing of our new transposer is presented below. There are two versions in fact: the faster 64-bit and the slower 32-bit. It would be natural to code such regular algorithm with a bunch of macros and templates to avoid trivial mistakes. However, I kept the code uncondensed for ease of comprehension. While reading the listing, keep in mind that x86_64 is a little endian architecture. Left-most elements in the image of 8x8 transposer are loaded into least significant bytes of the 64-bit registers. That’s why all masks and directions of shifts are also reversed.

The outcome is x1.8 boost in performance compared to the previous block prefetch algorithm. Now it takes less than one cycle to process single element.

As the final remark I would like to emphasize that in this algorithm we used block decomposition three times, and every time we did it for different reasons:

It should be obvious by now that dividing matrices into blocks and processing them independently is the key to designing any efficient matrix algorithm, not only transposer.

As was mentioned earlier, block algorithm becomes more efficient as word size increases. General purpose registers are only 64 bit long, so there is nothing more we can do about them. On the other hand, x86_64 carries dedicated SIMD registers and instructions with widths of 128, 256 and 512 bits depending on CPU generation. For example, with 256 bit registers we could create efficient 32×32 transposer by using log₂(32)=5 levels of 2×2 block decomposition. Alas, working with SIMD extensions is much more challenging than issuing + or >>. A multitude of high-level instructions allows to do the same thing in different ways and with different performance results. Atop of that, availability of SIMD extensions heavily depends on CPU architecture and generation. Here we will focus on relatively modern generations of x86_64 CPUs.

First let’s do short review of what SIMD extensions x64_64 CPUs offer. The very first widely used SIMD extension for x86 CPUs was MMX ("multimedia extensions"), now already forgotten. It introduced a set of new registers which could be treated as vectors of 8/16/32-bit integers and a dedicated instruction set that allowed to perform operations on the elements of these registers in parallel. As time passed, more and more extensions were added. Newly released extensions expanded register width, increased number of registers or introduced brand new instructions. Extensions below are grouped by target vector width and are listed in historical order:

I will assume that we have access to AVX2 but not the newer extensions. AVX2 gives us access to 256 bit registers and integer instructions. This means that we will be able to create efficient 32x32 transposer by using log₂(32)=5 levels of block decomposition.

As before, we will work level by level, row by row. Every row will be represented by a variable of __m256i type, which we will treat as a vector of type uint8_t[32]. This type is defined in immintrin.h header, shipped with all major compilers (intel, gcc, clang) which support AVX2. Variable of this type can be naturally mapped to a 256-bit AVX register by the compiler in the same way as uint64_t variable can be mapped to a general purpose 64-bit register. The same header also provides the intrinsics to work with these registers. Each intrinsic is mapped to respective CPU instruction, thus making it unnecessary to manually write assembler code. In our code we will rely on only three instructions: _mm256_shuffle_epi8, _mm256_blendv_epi8, _mm256_permute2x128_si256. In addition, two more instructions will be used implicitly: _mm256_load_si256 and _mm256_store_si256. They move data between memory and 256-bit registers. We do not need to call them manually since compilers can insert such instructions automatically on encountering _mm256i r = *ptr; and *ptr = r; respectively. Semantics of the other three instructions is more complicated and is provided below.

Shuffle. _mm256_shuffle_epi8 instruction accepts a source 256-bit register and a 256-bit control register, and, treating input as a vector of 8-bit values, shuffles them according to the control register, returning the result. Control register specifies for every index (0..31) of the output vector either the index of the source vector where to copy an element from, or carries a special value indicating that target location must be zeroed. Such semantics makes shuffle a very powerful instruction. If this instruction didn’t have any limitations, it would be possible to do any of the following with just single short-running instruction:

In contrast, coding the same operations with general purpose instructions would take a long sequence of mask, shift and OR instructions. Alas, _mm256_shuffle_epi8 has one limitation: elements may be shuffled only within every left and right 128-bit lane but not across the lanes. It is not possible, for example, to copy an element from src[30] (left lane) to dst[12] (right lane). So, technically speaking, _mm256_shuffle_epi8 is a pair of independent general shuffle operations applied to left and right 128-bit lanes.

In our algorithm, we will use shuffle instruction in the first four stages (out of five) to reposition elements of the registers holding rows of the matrix from the previous stage. In the first level, we swap adjacent elements, in the second level — adjacent groups of two elements each, and so on. Diagram below demonstrates the application of the shuffle during the third stage, when we swap adjacent groups of four elements each.

Permute. In the fifth stage of the algorithm we are transposing 2x2 block matrices where each block is 16x16, effectively meaning that we have to swap lanes. There is not way in AVX2 to perform general shuffle that requires to move data between lanes with a single instruction. But since our case is very special one — just swap the lanes — we can do this with another instruction — _mm256_permute2x128_si256. Given two input data registers a and b and a control register, this instruction allows to construct a register where each of the two 128-bit lanes takes any of the four options: a[0:127], a[128:255], b[0:127], b[128:255]. Setting control register to 0x01 makes permute instruction to return a with swapped lanes; b can be anything in our case.

Blend. Blend is another highly useful instruction. It accepts two source registers and a mask, and constructs result where every i-th element is selected from either src1[i] or src2[i] based on the mask.

In total, a combination of shuffle/permute and blend allows us to build a value consisting of almost arbitrary selection of elements (remember cross-lane restriction) from two input registers, a and b, in the following way:

Using this rule we are ready to build 32x32 transposer. There will be log2(32)=5 levels. Recall that at each level we build and intermediate matrix where each row is the result of computations applied to two rows of the matrix from the previous level. Every such computation can be expressed by the same sequence of operations. The difference will be in which rows we take as input, which masks we apply, and on the final level instead of shuffle we have to use permute. Image below demonstrates how to compute the first two rows on the very first stage, when we transpose each 2x2 matrix inside 32x32 matrix. What originally was , now becomes . And what was , now becomes . At the second stage the same sequence of operations applies, but now we work on pairs of elements. At the fifth stage we are working on groups of 16 elements, and instead of shuffle we use permute.

On the global scale full sequence of computations can be expressed with the following network. Each vertex corresponds to a 256-bit value that holds entire row of 32 elements of some matrix, additionally associated with the instruction that was used to compute this value. Edges encode dependencies which must be computed before we can compute target value itself. In the first column of vertices we load the rows of the src matrix. Next five pairs of columns correspond to five stages of the algorithm. First column in each pair is the result of shuffle (or permute), second column — blendv. Thus, the values computed in the penultimate column contain the rows of the transposed matrix, and the final column is dedicated to saving them into the memory that backs dst. Top vertices standing apart from the remaining network correspond to the loading of 256-bit masks, which we need to pass to shuffle and blendv. These masks should be computed beforehand.

Writing code for such complex transposer is not a simple task. Manual coding in entirely unrolled format (as Vec64) is out of the question because it would require to write more than 400 lines of highly repetitive code. Using nested loops (level by level, block by block, row by row) is definitely possible and the transposer would be quite compact, but its performance would be at the mercy of the compiler. For example, its choice which loops to unroll and to what extent will have impact on how well execution units will be loaded. Third alternative is to use code generation with the help of some scripting language, and probably it is the cleanest way to create medium-sized kernels such as ours.

To make use of code generation, we need to explicitly instantiate the above network in terms of vertices and edges. With each vertex we associate a piece of code that consists of executing single vector instruction and writing the result into uniquely-named variable. Arguments of the instructions are the variables associated with the left-adjacent vertices.

One of the benefits of such approach is that it is possible to play with different orderings of the vertices. Any valid topological sorting will produce the same result but performance may vary. Computations may be performed column by column, in a n-ary tree manner, according to a randomly generated topological sorting, or in any other clever way you can think of. A priory it is hard to say which sorting will be better in terms of performance due to microarchitectural details. Still, some general rules can be applied. On the one hand, we would like to generate code that at any step keeps number of already computed and still required in the future values as few as possible. This is because compiler maps each variable to some register, and when it runs out of registers (for AVX2, there are only 16 of them), it has to perform a spill: a piece of code that saves value into stack and restores its back later when it is needed. Obviously, the fewer spills — the better. Another thing we need to consider is that invocations of the same instruction should be spread. Issuing 32 identical instructions one by one is not the best idea, since for each instruction type, there is only a few of the execution units which can handle it. If we issue many identical instructions one by one, other execution units will starve for work. Much better idea is to interleave instructions of different types, raising the chances that they will be executed in parallel.

Long story short, the solution which appears to work well empirically is to perform all computations level-by-level, similar to scalar vectorization. At every step two rows of the matrix from the previous level serve as input and four instructions are applied to them: two shuffles and two blends, thus producing two rows of the next-level matrix. This is repeated level by level, pair of rows by pair of rows. The order in which values are computed according to these rules is denoted by the numbers near the vertices in the network above.

Code-generated 32x32 transposer must be plugged into the overall algorithm, again with carefully chosen strategy for preloading. Since we employ code generation, the natural choice is to spread all 16 preload instructions at equally-long chunks of the generated code. Why 16? Recall that we need 4 invocations of the 32x32 transposer to finish single 64x64 block. Hence during each invocation we need to preload 16 rows of the next block. There are 393 instructions in total, so we insert preload each 393/16=24 instructions into generated code. In the code sample, three origin variables correspond to the offsets of the 32x32 submatrices we are working on: source matrix, destination, and prefetch respectively. stride values indicate the number of bytes that must be added to a corresponding pointer to make it point to the same column of the next row. For our current transposer both strides must be set to N.

Switching from general purpose instruction set to AVX turned out to be quite complex. But it comes with a benefit of further reduction in running time by 15-30% compared to 64-bit version:

This trend should continue with increase in vector size, provided that performance is not stuck in memory access and that latency of longer-vector instructions is the same as of their shorter counterparts. Unfortunately, vector sizes can’t be increased indefinitely. In the majority of horizontal vector instructions, such as shuffle, every element of the output register depends on every element of the input register. Such logic requires O(log(n)) stages of digital circuitry, which increases latency and hence defeats the purpose of increasing vector size. We already saw that shuffle in fact is not a purely 256-bit horizontal SIMD instruction, but a pair of simultaneously executing 128-bit shuffles. For CPUs, it is quite unrealistic that vector instructions will be increased from 512 (provided by AVX-512) to longer values.

Although performance has improved, thorough observation of PMU counters reveals that memory subsystem is again the bottleneck and not the SIMD code. According to resource_stalls.sb counter, there are almost 16 stall cycles per each 64 bytes of input data. Values in the table were collected by transposing three matrices of N=46400.

Memory and the execution are the most common culprits of suboptimal performance. One of them is always the bottleneck. Fixing memory issues makes execution the bottleneck, and optimizing execution uncovers more issues with memory. This process can be repeated many times up to the point when both are working close to their theoretical limits.

The fact that we get nearly sixteen SB stalls per each cache line is unexpected. Each two invocations of 32x32 transposer fully assemble 32 cache lines, which should definitely fit into SB (Skylake has 56 entries). And since they are assembled in full, we expect that they will be flushed out without sending any load requests to RAM first. And yet there are SB stalls. To figure out what is going on, we can look to offcore_requests.demand_rfo counter. Its value, normalized by number of processed cache lines, reveals that there is one RFO message per each cache line. RFO messages are sent to other cores to agree on the ownership of cache lines. There can be only one exclusive owner of given cache line at any given time, or otherwise different cores may have different versions of the same cache line, ruining the correct behaviour of the programs which concurrently modify the data. Once the responses from all other cores are received, the requesting core becomes the sole owner of the cache line and has the right to modify it locally. In case of our transposer it is sure that the region of dst where we are writing to is not yet in L1d or L2, leading us to the conclusion that current core is not yet the owner of the respective cache lines, and that’s why it has to send the RFOs before the modification can be applied. Apparently, waiting for the RFO responses is the bottleneck of our SIMD transposer.

Such type of issues can be circumvented if we write the result of 64x64 transpose not directly into dst matrix, but into temporary L1d-local buffer. After transpose of the 64x64 block is complete, we can move data from the buffer into dst using non-temporal hint. The idea that we have to write data twice is counterintuitive, but it reduces waiting times for inter-core communication. During the first stage we write data into temporary buffer, which is reused for each 64x64 block and that’s why it settles down in L1d, with all of its cache lines being privately owned by the current core. The latter means that when some cache line of this buffer needs to be updated, there is no need to send RFO and wait for responses from other cores: cache line is marked in L1d as privately owned by this core with no possibility of other cores having a copy of it. In the second stage we have to copy data from buffer into dst. Using ordinary move instructions would generate RFO message for each cache line and wait while all other cores respond that they have relinquished ownership if they had one. Waiting for RFO responses is what is currently limiting performance of the transposer. Using non-temporal move instructions, on the other hand, relaxes the ownership rules. No RFOs are sent and, since cache lines are modified in full, they can be moved directly to memory.

Of course, non-temporal moves do not come without limitations, or otherwise they would be used ubiquitously. First, to create any useful effect, they have to be issued in burst covering the whole cache line. Second, you have to be sure that no other threads are accessing the data concurrently with non-temporal move. The latter is a common contract in matrix-oriented code.

The changes we need to make to our program are small: preallocate 64x64 buffer, instruct 32x32 SIMD transposer to write result into this buffer, and finally use a burst of _mm256_stream_si256 instructions to copy data from the buffer into dst. 32x32 transpose kernel remains the same.