# How do BLAS libraries implement support for transposed matrices?

I'm trying to understand how BLAS libraries implement fast GEMM with support for transposed matrices. Say, I'm only operating on square matrixes (with dimensions n by n) and I update a submatrix of c in the following manner c[x:x+6][y:y+16] += a[x:x+6][slice_start:slice_end] * b[slice_start:slice_end][y:y+16] with the following kernel

void kernel(float * a, float * b, float * c, int x, int y,
int slice_start, int slice_end, int n,
bool trans_a, bool trans_b) {

float result[6][16] = {0};

for (int k = slice_start; k < slice_end; k++) {
for (int i = 0; i < 6; i++) {
float aa = a[(x + i) * n + k];
for (int j = 0; j < 16; j++) {
result[i][j] += aa * b[k * n + y + j];
}
}
}

for (int i = 0; i < 6; i++) {
for (int j = 0; j < 16; j++) {
c[(x + i) * n + y + j] += result[i][j];
}
}
}


How can I add support for transposed matrices? If a is transposed, I can simply switch the indexes for aa retrieval from a[(x + i) * n + k] to a[k * n + (x + i)], which wouldn't result in worse performance, since elements of a are accessed much less frequently. However, for the same reason, a similar approach wouldn't work for b. Not to mention that this approach for a accesses elements outside of the kernel window, hindering parallelisation across many threads.

So I was wondering if someone could suggest an approach here? Is there a simple trick? Or is a separate kernel needed for each combination of transpositions necessary? Or is support for transposed matrixes handled outside of the micro-kernel somehow? Or something else entirely?

thanks :)

Blas kernels typically don't deal with transposed matrices. Instead, they will repack matrices (transposed or not) outside of the kernel. See https://www.cs.utexas.edu/users/flame/pubs/blis3_ipdps14.pdf for details of what the loop nest generally looks like.

Specifically, from the paper (highlighted for emphasis):

Our description starts with the outer-most loop, indexed by jc. This loop partitions C and B into (wide) column panels. Next, A and the current column panel of B are partitioned into column panels and row panels, respectively, so that the current column panel of C (of width nc) is updated as a sequence of rank-k updates (with k = kc), indexed by pc. At this point, the GotoBLAS approach packs the current row panel of B into a contiguous buffer, ˜B. If there is an L3 cache, the computation is arranged to try to keep ˜B in the L3 cache. The primary reason for the outer-most loop, indexed by jc, is to limit the amount of workspace required for ˜B, with a secondary reason to allow ˜B to remain in the L3 cache. Now, the current panel of A is partitioned into blocks, indexed by ic, that are packed into a contiguous buffer, ˜A. The block is sized to occupy a substantial part of the L2 cache, leaving enough space to ensure that other data does not evict the block. The GotoBLAS approach then implements the “block-panel” multiplication of ˜A ˜B as its inner kernel, making this the basic unit of computation

It's not explicitly spelled out, but the pace where you pack panels and blocks is generally where the parts are transposed if necessary.

• i've actually read parts of this paper independently of this, but never crossed that part i guess... could you point me to a specific page in this document where this is described?
– ilya
Feb 21 at 4:05
• updated with selected quotes Feb 21 at 16:16
• so essentially (assuming we're in row major order), if a matrix isn't transposed, you partition it into horizontal panels (which means your micro-kernel has linear memory access), and if a matrix is transposed, we partition it into vertical panels? but even so the data isn't contiguous for transposed matrices? or are individual elements supposed to be moved around (not just blocks or panels) in order to ensure contiguous memory access?
– ilya
Feb 21 at 20:07
• the latter. the packing crates panels such that the microkernel will have contiguous memory by moving things around as needed. Feb 21 at 21:06
• how come there's virtually no difference in performance when performing GEMM on non-transposed vs transposed matrices? i just did 20000 iterations on two 1024x1024 matrices, and transposed operations actually came out slightly faster shouldn't the inverse be true given that you have to move individual floats around (which is pretty slow), not to mention that requiring additional temporary storage for transposing panels violates the BLAS spec of additional storage?
– ilya
Feb 21 at 21:54