4 ms·

> But in the end, it's FORTRAN all the way down. Even in Julia. That's not true. None of the Julia differential equation solver stack is calling into Fortran a

> But in the end, it's FORTRAN all the way down. Even in Julia. That's not true. None of the Julia differential equation solver stack is calling into Fortran anymore. We have our own BLAS tools that outperform OpenBLAS and MKL in the instances we use it for (mostly LU-factorization) and those are all written in pure Julia. See https://github.com/YingboMa/RecursiveFactorization.jl https://github.com/YingboMa/RecursiveFactorization.jl, https://github.com/JuliaSIMD/TriangularSolve.jl https://github.com/JuliaSIMD/TriangularSolve.jl, and https://github.com/JuliaLinearAlgebra/Octavian.jl https://github.com/JuliaLinearAlgebra/Octavian.jl. And this is one part of the DiffEq performance story. The performance of this of course is all validated on https://github.com/SciML/SciMLBenchmarks.jl https://github.com/SciML/SciMLBenchmarks.jl

gnufx 5y ago
Is the recursive comparison against OpenBLAS using ReLAPACK or the normal(?) one? I don't remember the story on building with that or know the different algorithms. Does Julia not use OpenBLAS at all now?
- ChrisRackauckas 5y ago
  Julia defaults to OpenBLAS but libblastrampoline makes it so that `using MKL` flips it to MKL on the fly. See the JuliaCon video for more details on that (https://www.youtube.com/watch?v=t6hptekOR7s https://www.youtube.com/watch?v=t6hptekOR7s). The recursive comparison is against OpenBLAS/LAPACK and MKL, see this PR for some (older) details: https://github.com/YingboMa/RecursiveFactorization.jl/pull/28 https://github.com/YingboMa/RecursiveFactorization.jl/pull/2... . What it really comes down to in the end is that OpenBLAS is rather bad, and MKL is optimized for Intel CPUs but not for AMD CPUs, so when the best CPUs are now all AMD CPUs, having a new set of BLAS tools and mixing that with recursive LAPACK tools is either as good or better on most modern systems. Then we see this in practice even when we build BLAS into Sundials for 1,000 ODE chemical reaction networks (https://benchmarks.sciml.ai/html/Bio/BCR.html https://benchmarks.sciml.ai/html/Bio/BCR.html).
  - gnufx 5y ago
    The question was whether the comparison was against OpenBLAS built with recursive LAPACK routines (not the default, I've now checked). I don't know whether that would use the same algorithm. I don't see great attraction in "flipping" to MKL. (In general I switch BLASes by dynamic linking if necessary.) I'm puzzled by the value of "rather bad" for OpenBLAS, at least for BLAS [1], but it seemed OK against MKL for LAPACK-based R benchmarks I measured on sandybridge. That said, I don't understand why people avoid AMD's BLAS/LAPACK (BLIS/libflame with enhancements intended for merging) on AMD systems. In this context: BLIS studiously avoids Fortran for some reason, even for test the Fortran interface. I have nothing against Julia, incidentally; I have the original Dylan book. 1. https://git.sr.ht/~fx/blis/tree/performance/docs/Performance.md https://git.sr.ht/~fx/blis/tree/performance/docs/Performance...
    - ChrisRackauckas 5y ago
      It doesn't look like you're measuring factorization performance? OpenBLAS matrix-matrix multiplication is fine, it just falls apart when going to things like Cholesky and LU. Also it falls apart on smaller problems, like <100x100 cases. > not the default, I've now checked Whatever the Julia default build is doing, so probably not the recursive LAPACK routines then if that's how it's being built. If there's a better default that's worth an issue. > That said, I don't understand why people avoid AMD's BLAS/LAPACK There just isn't a BLIS wrapper into Julia right now, and it's much easier to just write new BLAS tools than to build wrappers IMO. It makes it very easy to customize to nonstandard Julia number types too. But I do think that BLIS is a great project and I would like to see it replace OpenBLAS as the default. There's been some discussion to make it as easy as MKL (https://github.com/JuliaLinearAlgebra/BLIS.jl/issues/3 https://github.com/JuliaLinearAlgebra/BLIS.jl/issues/3).
      - gnufx 5y ago
        OK, no, I measured level 3 BLAS, and it sounded like a blanket dismissal of OpenBLAS. I'm surprised at it falling apart at all, but I haven't considered its LAPACK much and I value understanding from measurement. The set of R LA benchmarks I ran on sandybridge in response to theories about the Microsoft R distribution did include Cholesky. (I think the operations are a thin wrapper over LAPACK ones.) In that case OB was actually better than MKL for Cholesky with the versions I had then, so I'm puzzled. Probably I should revisit that sometime. OB seems to me a better general default than BLIS, though its threading implementation may not be robust. BLIS is even broken on POWER9, and I don't remember whether the dynamic micro-architecture selection ever added for other than x86. For really small matrices on x86, I'd front BLAS with libxsmm, but I don't remember where it starts winning, or know the threshold for BLIS' non-packed kernels. Anyway, thanks for the info. I should be embarrassed never to have made time to look seriously at Julia.

gnufx 5y ago

Is the recursive comparison against OpenBLAS using ReLAPACK or the normal(?) one? I don't remember the story on building with that or know the different algorithms. Does Julia not use OpenBLAS at all now?

ChrisRackauckas 5y ago
Julia defaults to OpenBLAS but libblastrampoline makes it so that `using MKL` flips it to MKL on the fly. See the JuliaCon video for more details on that (https://www.youtube.com/watch?v=t6hptekOR7s https://www.youtube.com/watch?v=t6hptekOR7s). The recursive comparison is against OpenBLAS/LAPACK and MKL, see this PR for some (older) details: https://github.com/YingboMa/RecursiveFactorization.jl/pull/28 https://github.com/YingboMa/RecursiveFactorization.jl/pull/2... . What it really comes down to in the end is that OpenBLAS is rather bad, and MKL is optimized for Intel CPUs but not for AMD CPUs, so when the best CPUs are now all AMD CPUs, having a new set of BLAS tools and mixing that with recursive LAPACK tools is either as good or better on most modern systems. Then we see this in practice even when we build BLAS into Sundials for 1,000 ODE chemical reaction networks (https://benchmarks.sciml.ai/html/Bio/BCR.html https://benchmarks.sciml.ai/html/Bio/BCR.html).
- gnufx 5y ago
  The question was whether the comparison was against OpenBLAS built with recursive LAPACK routines (not the default, I've now checked). I don't know whether that would use the same algorithm. I don't see great attraction in "flipping" to MKL. (In general I switch BLASes by dynamic linking if necessary.) I'm puzzled by the value of "rather bad" for OpenBLAS, at least for BLAS [1], but it seemed OK against MKL for LAPACK-based R benchmarks I measured on sandybridge. That said, I don't understand why people avoid AMD's BLAS/LAPACK (BLIS/libflame with enhancements intended for merging) on AMD systems. In this context: BLIS studiously avoids Fortran for some reason, even for test the Fortran interface. I have nothing against Julia, incidentally; I have the original Dylan book. 1. https://git.sr.ht/~fx/blis/tree/performance/docs/Performance.md https://git.sr.ht/~fx/blis/tree/performance/docs/Performance...
  - ChrisRackauckas 5y ago
    It doesn't look like you're measuring factorization performance? OpenBLAS matrix-matrix multiplication is fine, it just falls apart when going to things like Cholesky and LU. Also it falls apart on smaller problems, like <100x100 cases. > not the default, I've now checked Whatever the Julia default build is doing, so probably not the recursive LAPACK routines then if that's how it's being built. If there's a better default that's worth an issue. > That said, I don't understand why people avoid AMD's BLAS/LAPACK There just isn't a BLIS wrapper into Julia right now, and it's much easier to just write new BLAS tools than to build wrappers IMO. It makes it very easy to customize to nonstandard Julia number types too. But I do think that BLIS is a great project and I would like to see it replace OpenBLAS as the default. There's been some discussion to make it as easy as MKL (https://github.com/JuliaLinearAlgebra/BLIS.jl/issues/3 https://github.com/JuliaLinearAlgebra/BLIS.jl/issues/3).
    - gnufx 5y ago
      OK, no, I measured level 3 BLAS, and it sounded like a blanket dismissal of OpenBLAS. I'm surprised at it falling apart at all, but I haven't considered its LAPACK much and I value understanding from measurement. The set of R LA benchmarks I ran on sandybridge in response to theories about the Microsoft R distribution did include Cholesky. (I think the operations are a thin wrapper over LAPACK ones.) In that case OB was actually better than MKL for Cholesky with the versions I had then, so I'm puzzled. Probably I should revisit that sometime. OB seems to me a better general default than BLIS, though its threading implementation may not be robust. BLIS is even broken on POWER9, and I don't remember whether the dynamic micro-architecture selection ever added for other than x86. For really small matrices on x86, I'd front BLAS with libxsmm, but I don't remember where it starts winning, or know the threshold for BLIS' non-packed kernels. Anyway, thanks for the info. I should be embarrassed never to have made time to look seriously at Julia.