4 ms·

I think Eigen's template strategy is probably hard to beat from the perspective of combining operations. How well Fortran does is probably mostly up to the comp

by celrod 6y ago

I think Eigen's template strategy is probably hard to beat from the perspective of combining operations. How well Fortran does is probably mostly up to the compiler implementation. In some of my benchmarks, gfortran sometimes seems to fuse the operations and end up much slower than if it has performed them separately in succession, but I wouldn't be surprised if compiling with `-fexternal-blas` (and linking MKL) would've solved that. If you want to try external BLAS/LAPACK with Eigen, I'd look at: https://eigen.tuxfamily.org/dox/TopicUsingBlasLapack.html https://eigen.tuxfamily.org/dox/TopicUsingBlasLapack.html I have `A * B`, `A * B'`, `A' * B` and `A' * B'` small-single-threaded-matmul benchmarks here: https://chriselrod.github.io/LoopVectorization.jl/latest/examples/matrix_multiplication/ https://chriselrod.github.io/LoopVectorization.jl/latest/exa... I compared triple nested loops with Clang, icc, ifort, gfortran, Julia, and LoopVectorization.jl with matmul routines from ifort, gfortran, OpenBLAS, MKL, and Eigen. While gfortran's builtin hit over 40 GFLOPS with `A * B` and `A' * B'`, it failed to get half that if only one argument was transposed. I'm supposed awkward fusing at the start of this post because if it had done one after the other, it should have still hit >40 GFLOPS when only 1 matrix was transposed.

kergonath 6y ago
These benchmarks are very interesting! From my experience gfortran’s MATMUL is very good for small matrices, but OpenBLAS gets better from ~10x10 elements. Not quite sure that’s the same as your benchmark; its colour code is a bit confusing. Would you mind making the data available? 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.
- celrod 6y ago
  My color coding is bad (and there's an open issue about it), but I haven't come up with a great solution. One idea I should add was to make everything relative to LoopVectorization, but that wouldn't help for those interested in other comparisons like Eigen vs gfortran. 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) the compiler will take advantage of the known dimensions and inline the call, making it much faster. > 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.
  - kergonath 6y ago
    I think what the plots need is an option to show or hide each curve. It’s a rather large increase in complexity compared to just generating images, though.
certik 6y ago
Thanks for the benchmarks. Do you see some intrinsic reason why a Fortran compiler couldn't do these optimizations? I would think it would know all the information so it would be the ideal place to optimize it.
- celrod 6y ago
  There is no reason it could not. Those optimizations just have to be implemented. Flang (the one merged into LLVM) is using MLIR, which has all the required code-gen abilities. That just leaves the cost modeling / deciding which optimizations/transformations to apply.
  - certik 6y ago
    Yes, once MLIR is more mature, we plan to have a backend for LFortran to also use MLIR.
    - celrod 6y ago
      For BLAS in particular, this paper can give you an idea of some of MLIR's capabilities: https://arxiv.org/pdf/2003.00532.pdf https://arxiv.org/pdf/2003.00532.pdf (But maybe you already know them better than I do.) 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.
waltpad 6y ago
If you have a look at the different fortran implementations of BLAS gemm (matrix multiplication), you'll see that the transposed matrix cases are treated specifically. In fact, IIRC, the gemm function has flags to indicate for each matrix if it is transposed.
- celrod 6y ago
  That's what I did when calling OpenBLAS and MKL, but I confess I don't know the internal details of a non-inlined `matmul` call in gfortran when you don't use `-fexternal-blas`. 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.

by celrod 6y ago

kergonath 6y ago

These benchmarks are very interesting! From my experience gfortran’s MATMUL is very good for small matrices, but OpenBLAS gets better from ~10x10 elements. Not quite sure that’s the same as your benchmark; its colour code is a bit confusing. Would you mind making the data available? 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.

celrod 6y ago
My color coding is bad (and there's an open issue about it), but I haven't come up with a great solution. One idea I should add was to make everything relative to LoopVectorization, but that wouldn't help for those interested in other comparisons like Eigen vs gfortran. 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) the compiler will take advantage of the known dimensions and inline the call, making it much faster. > 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.
- kergonath 6y ago
  I think what the plots need is an option to show or hide each curve. It’s a rather large increase in complexity compared to just generating images, though.

certik 6y ago

Thanks for the benchmarks. Do you see some intrinsic reason why a Fortran compiler couldn't do these optimizations? I would think it would know all the information so it would be the ideal place to optimize it.

celrod 6y ago
There is no reason it could not. Those optimizations just have to be implemented. Flang (the one merged into LLVM) is using MLIR, which has all the required code-gen abilities. That just leaves the cost modeling / deciding which optimizations/transformations to apply.
- certik 6y ago
  Yes, once MLIR is more mature, we plan to have a backend for LFortran to also use MLIR.
  - celrod 6y ago
    For BLAS in particular, this paper can give you an idea of some of MLIR's capabilities: https://arxiv.org/pdf/2003.00532.pdf https://arxiv.org/pdf/2003.00532.pdf (But maybe you already know them better than I do.) 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.

waltpad 6y ago

If you have a look at the different fortran implementations of BLAS gemm (matrix multiplication), you'll see that the transposed matrix cases are treated specifically. In fact, IIRC, the gemm function has flags to indicate for each matrix if it is transposed.

celrod 6y ago
That's what I did when calling OpenBLAS and MKL, but I confess I don't know the internal details of a non-inlined `matmul` call in gfortran when you don't use `-fexternal-blas`. 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.