25 ms·

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 w

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.

kergonath 6y ago
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.
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
adamnemecek 6y ago
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.
- waynecochran 6y ago
  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.
mhh__ 6y ago
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?
- waynecochran 6y ago
  Yeah, you got it. The diagonal() matrix is not actually created since you really don't need all the elements in the final computation.
cinquemb 6y ago
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 http://arma.sourceforge.net/download.html
antoinealb 6y ago
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 https://eigen.tuxfamily.org/dox/TopicUsingBlasLapack.html

kergonath 6y ago

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.

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

adamnemecek 6y ago

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.

waynecochran 6y ago
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.

mhh__ 6y ago

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?

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

cinquemb 6y ago

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 http://arma.sourceforge.net/download.html

antoinealb 6y ago

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 https://eigen.tuxfamily.org/dox/TopicUsingBlasLapack.html