What it takes to transpose a matrix

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 log2(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 log2(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.