Certainly, gfortran’s default implementation works well as a quick and easy solution for small vectors and matrices outside of hot paths. Also, I remember discussions about improving matmul(transpose(A),B) somewhat recently. I don’t remember for which version it was, but it’s the sort of things that is improved regularly.

I would love an option to align arrays in gfortran. It’s very important to take advantage of automatic vectorisation but I haven’t found a reliable way to do it.

I'm on gfortran 10.1.1, so I'd only be missing out on very recent improvements (i.e., to trunk).

Note that these benchmarks were dynamically sized. If you write a fortran program like

  real(kind=8), dimension(10,10) :: A, B, C
  ! fill A and B somehow
  C = matmul(A, B)

> I would love an option to align arrays in gfortran. It’s very important to take advantage of automatic vectorisation but I haven’t found a reliable way to do it.

I wish more compilers could also use masks to vectorize without padding. If multiplying 7x7 matrices with AVX512, the obvious solution is to just mask the loads/stores of columns, without the need for padding. Of course, padding would also ensure all your loads/stores are aligned, which can be nice.

LoopVectorization can't do many of these yet, so its performance will fall off a cliff shortly after the largest size on the plots (and at much smaller sizes for CPUs with a smaller L2 cache). I had to add code to perform packing/tiling in my actual matmul code on top of what it did. So that MLIR can generate that sort of code already looks promising.

Still, the work of telling it what to do isn't easy.

I'm not involved in any of those projects, so everything I say here is pure speculation. But I imagine pragmas and the like would be important whenever the compiler doesn't know the sizes at compile time. Otherwise, you probably don't want it to generate massive amounts of code through multiple extra blocking loops, massive unrolling in a main kernel, and multiple clean up kernels for every random loop nest.

Just writing three loops and letting the compiler optimize it was much faster for `A * B'`, so it must be a pretty naive implementation getting called.