Hacker News Reader

Comment by waynecochran | Hacker News Reader

I use C++ Eigen extensively. I wonder if it feasible to write a C++ programs that hands off blocks of linear algebra code to Fortran subroutines. Eigen does a wonderful job of late binding and lazy evaluation so that something like the following is fairly efficient:

     (L.transpose() * P).diagonal().array().square().mean();

which computes the squared average of column by column dot products between matrices L and P. I assume Fortran, since matrices are native types, does this kind of thing even better.

Aside: Fortran 77 was the first programming language I learned in 1985 and I hated programming for awhile after that. We called it the "F" language. Now that I am a bit older and do a lot of linear algebra coding... I wonder if I would like it now.

View on HN

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

I have `A * B`, `A * B'`, `A' * B` and `A' * B'` small-single-threaded-matmul benchmarks here: 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.

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.

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.

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.

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.

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.

Yes, once MLIR is more mature, we plan to have a backend for LFortran to also use MLIR.

For BLAS in particular, this paper can give you an idea of some of MLIR's capabilities: 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.

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.

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.

You can define EIGEN_USE_BLAS and then Eigen will delegate matrix operations to BLAS routines (such as gemm) instead of its own engine [0]. Then, you could link your program to a BLAS implementation written in Fortran.

However, when I last researched this for my previous work, OpenBLAS was the fastest open source implementation and is written in a mix of C and assembly, not Fortran. The repository contains quite a lot of Fortran, but it is mostly for the LAPACK implementation, which is not part of BLAS.

In our testing we got a nice speedup of our signal processing pipeline when enabling EIGEN_USE_BLAS, but this was on ARM64, so your experience may differ.

[0] https://eigen.tuxfamily.org/dox/TopicUsingBlasLapack.html

You can do this with Julia. Julia has an insanely good interop with essentially any language you can name, including cpp and Fortran. You will go through julia but hey, for linear algebra it's the shit.

It would be interesting to compare using Eigen's BCDSVD solv er vs handing it off to julia's library to compute the SVD and handing it back.

If you use LAPACK there’s a chance that some Fortran code will be run at some point. For example OpenBLAS does it behind the scenes of you give it a matrix to diagonalise.

I have been using armadillo now vs Eigen for a couple of years, which can be compiled against OpenBLAS and LAPACK, where matrix ops are handled by fortran routines.

http://arma.sourceforge.net/download.html

I've never used Eigen, could you explain specifically what lazy evaluation implies for this expression? Is it to do with the diagonal() call constraining the amount of work to be done?

Yeah, you got it. The diagonal() matrix is not actually created since you really don't need all the elements in the final computation.

Reply on news.ycombinator.com